鄭建麗,曹建軍,劉 洋
(1.中國水產(chǎn)科學(xué)研究院漁業(yè)機(jī)械儀器研究所,上海 200092;2.遼寧師范大學(xué)計算機(jī)與信息技術(shù)學(xué)院,遼寧 大連 116081)
船載雷達(dá)是放置在船上用來觀測船舶周圍各種水面目標(biāo)和即時氣象環(huán)境的重要觀測工具。在船舶航行中遇到霧、霾、雨等現(xiàn)象導(dǎo)致能見距離受到影響,無法憑目視發(fā)現(xiàn)目標(biāo)的條件下,雷達(dá)表現(xiàn)出受環(huán)境影響小、觀測距離遠(yuǎn)的優(yōu)勢。目前,常用的航海雷達(dá)有SPERRY-MK 系 列、JRC-JMA 系 列 和FURUNO-FAR系列等。根據(jù)國際電工委員會IEC62388—2013 標(biāo)準(zhǔn),船載雷達(dá)的顯示圖像應(yīng)能夠進(jìn)行連續(xù)且平滑的雷達(dá)圖像更新[1]。但目前主流的雷達(dá)終端只能滿足部分觀測需求,還無法達(dá)到連續(xù)平滑的顯示效果,難以觀測海上近地大氣層一些快速生成和移動的強(qiáng)對流氣象目標(biāo)。從觀測對象的角度分析,海洋環(huán)境水汽熱量充足,下墊面平滑,氣象目標(biāo)形成和移動速度較快。約70%的龍卷風(fēng)從出現(xiàn)到消失的生命周期為1~10 min,海上冷風(fēng)的移動速度可達(dá)70 km/h 以上[2]。這些快速生成和變化的氣象目標(biāo)對雷達(dá)數(shù)據(jù)處理提出了高清顯示和更快速度的要求。從觀測平臺的角度分析,船載雷達(dá)因平臺穩(wěn)定性差和能量系統(tǒng)相對封閉等原因,探測距離一般為40~80 km。相較于陸基雷達(dá)通常200~300 km 的探測距離,船載雷達(dá)的目標(biāo)探測距離更短,因而有效應(yīng)對時間更少。
陸基雷達(dá)處理的數(shù)據(jù)量通常遠(yuǎn)遠(yuǎn)高于船載雷達(dá),但傳統(tǒng)陸基雷達(dá)對于提升數(shù)據(jù)處理速度的解決方案不適用于船載雷達(dá)。通常陸基雷達(dá)產(chǎn)品生成一般分為前端和后端兩部分,前端計算機(jī)用來控制雷達(dá)及其機(jī)構(gòu)并收集反射信號,后端計算機(jī)用來將雷達(dá)信號反演成圖像產(chǎn)品。在陸基雷達(dá)的反演過程中后端計算機(jī)通常為計算能力較強(qiáng)的超級計算機(jī),通過改善硬件的方式提高反演速度。但由于小型船舶制造成本和能量環(huán)境等限制因素,應(yīng)用改進(jìn)硬件來提高反演速度的方法并不適用于目前的漁船和游艇等小型船舶的情況。
目前,國內(nèi)海洋漁船一般采用基于(8~12 GHz)波段的小型雷達(dá),通過天線發(fā)射電磁波到海面,電磁波與雷達(dá)波長相當(dāng)?shù)拿?xì)波產(chǎn)生Bragg 散射,后向散射回波后被天線接收,形成“海雜波”,接收終端把連續(xù)接收的回波形成圖像序列,而波長較長的重力波通過對毛細(xì)波的調(diào)制作用表現(xiàn)在海雜波圖像上,從中提取風(fēng)、浪、流等海洋環(huán)境信息[3-5]。本文分析傳統(tǒng)雷達(dá)圖像反演過程和成像步驟,利用對雷達(dá)圖層封閉輪廓的重建,提出基于輪廓重建的船載雷達(dá)MAP-MRF 快速反演方法,解決漁船雷達(dá)圖像反演過程中大量重復(fù)計算缺失數(shù)據(jù)的重建插值過程。通過統(tǒng)計力學(xué)模型對MRF 中格點劃分進(jìn)行精確建模,利用MAP-MRF 方法對基于輪廓重建的雷達(dá)圖層邊界進(jìn)行估計。將整幅圖像的插值重建轉(zhuǎn)化為對圖層輪廓的精確分割,實現(xiàn)雷達(dá)圖像的快速反演,從而提高顯示終端的刷新率,實現(xiàn)對海上極端氣象目標(biāo)的及時發(fā)現(xiàn)和連續(xù)追蹤,滿足漁民對天氣氣象的變化情況特別是壞氣象實時了解的需求。在證明基于輪廓反演方法有效性的同時,驗證應(yīng)用基于統(tǒng)計力學(xué)模型的MAP-MRF 算法對雷達(dá)圖像進(jìn)行相位估計的可行性,為進(jìn)一步從概率模型方向處理雷達(dá)圖像提供新的分析思路和處理方法。主要內(nèi)容分為3 個部分:一是通過對漁船雷達(dá)數(shù)據(jù)和反演過程的分析,明確反演方法改進(jìn)重點和改進(jìn)內(nèi)容;二是對PPI 模式下的雷達(dá)基數(shù)據(jù)圖像應(yīng)用MAP-MRF 算法重建封閉輪廓進(jìn)行精確重建,進(jìn)而反演成像;三是將試驗結(jié)果與較有代表性的傳統(tǒng)反演方法在成像質(zhì)量和運(yùn)行速度上進(jìn)行分析和比較。
雷達(dá)圖像的重建可分為基于時間維度和基于空間維度的兩種過程,二者分別對應(yīng)雷達(dá)反演過程中時間分辨率(temporal resolution)和空間分辨率(spatial resolution)。本文主要討論單幀雷達(dá)圖像的空間超分辨率重建方法。
雷達(dá)的掃描方式是建立在球坐標(biāo)系中,但其展示和應(yīng)用通常在笛卡爾坐標(biāo)系下,因此雷達(dá)圖像的反演處理中包含了兩個坐標(biāo)系下的數(shù)據(jù)坐標(biāo)轉(zhuǎn)換。當(dāng)坐標(biāo)系變換后,雷達(dá)的輻射結(jié)構(gòu)會造成近雷達(dá)端數(shù)據(jù)采樣密集混疊,遠(yuǎn)雷達(dá)端數(shù)據(jù)采樣稀疏缺失的情況。因此,雷達(dá)通常會對坐標(biāo)系變換后的雷達(dá)圖像進(jìn)行重建或插值操作來改善以上問題。將雷達(dá)圖像重建算法分為雷達(dá)格點分析法和通用插值算法兩類。
雷達(dá)格點分析法是在雷達(dá)和氣象領(lǐng)域常用的數(shù)據(jù)插值和分析方法,如Barnes 分析法和Cressman 分析法等[6-9]。Barnes 分析法能夠把分布不均勻的數(shù)據(jù)映射到規(guī)則網(wǎng)格坐標(biāo)系的指數(shù)權(quán)重,在映射到笛卡爾坐標(biāo)系時通常應(yīng)用均一Barnes 插值,在插值到球坐標(biāo)系時一般應(yīng)用自適應(yīng)Barnes 插值。重建雷達(dá)圖像具有平滑、失真較小的優(yōu)點,此外Barnes 分析法還可用于基于站點序列的插值和糾錯。Cressman 分析法用于從基數(shù)據(jù)向格點空間的映射,通過加權(quán)影響半徑內(nèi)基數(shù)據(jù)對格點值進(jìn)行計算,采用迭代影響半徑的方式優(yōu)化精度依據(jù)影響半徑。Cressman 計算速度快且重建準(zhǔn)確率高,但當(dāng)映射空間密度高于基數(shù)據(jù)密度時算法穩(wěn)定性較弱,邊緣附近易出現(xiàn)極值點。
通用插值算法因應(yīng)用范圍廣,在性能和比較分析方面有比較深入的研究[10-14]。其中最常用的算法包括最鄰近插值法、雙線性插值法、基于小波的插值法,以及近年來被廣泛研究的基于機(jī)器學(xué)習(xí)理論的插值法等。其中最鄰近插值法是最簡單也是最快速的一種方法,計算復(fù)雜度最低,對于噪聲敏感。雙線性插值法是在實際應(yīng)用中被廣泛采用的算法,時間復(fù)雜度較低,峰值信噪比、均方差及主觀評價等較最鄰近插值法有明顯改善。基于小波的插值法是利用小波分解后的高頻信息預(yù)測高分辨率圖像,重構(gòu)高分辨率圖像的方法,是應(yīng)用變換域插值的一種主要方法。機(jī)器學(xué)習(xí)的圖像超分辨率重建法屬于基于樣例的重建,其中比較有代表性的如DONG W 等[15]基于稀疏編碼重建,以及SHI W 等[16]對卷積網(wǎng)絡(luò)的深度訓(xùn)練圖像重建等,對于普通圖像深度學(xué)習(xí)方法具有較好的重建效果,但在氣象雷達(dá)領(lǐng)域面臨著探測目標(biāo)通常不具有穩(wěn)定的結(jié)構(gòu)特征,重建內(nèi)容缺少對應(yīng)訓(xùn)練數(shù)據(jù)的問題。
在圖像重建的研究中,基于邊緣信息的重建是比較有代表性的方向。圖像邊緣表達(dá)了更加豐富的圖像信息,通過解決重建過程中的邊緣模糊效應(yīng)和鋸齒狀不連續(xù)等問題能夠有效改善整幅重建效果[17-18]。因此有效利用圖像邊緣信息是有效改善重建效果、提升重建效率的一個有效方向。受此啟發(fā),提出精確分割閾值圖層輪廓實現(xiàn)雷達(dá)圖像重建的方法。同時,漁船雷達(dá)圖像是灰度圖像,圖像中的艦船目標(biāo)和陸地區(qū)域灰度級較高,海面區(qū)域灰度值較低,目標(biāo)與目標(biāo)之間不會重疊,對此圖像重建方法具有一定的自適應(yīng)性。
雷達(dá)圖像描述的輪廓邊緣既對應(yīng)物理特征的等高線,又描述氣象目標(biāo)內(nèi)部降水粒子的尺度和密度分布。這些粒子懸浮在大氣中,其分布和運(yùn)動特征可以用大氣流體運(yùn)動模型描述。
此類模型主要是指基于納維—斯托克斯(Navier-Stokes,NS)方程的各類多相流模型,模擬方法以傳統(tǒng)計算流體力學(xué)方法為主,光滑粒子方法、格子玻爾茲曼方法(Lattice Boltzmann Method,LBM)應(yīng)用比較廣泛[19]。
2004 年,Dankert 從雷達(dá)圖像序列中反演海面風(fēng)場,研究了以光流運(yùn)動估計技術(shù)為基礎(chǔ)的風(fēng)場反演方法。此方法不需要校正雷達(dá)系統(tǒng)的相位,具有較強(qiáng)的實用價值。然而,由于反演風(fēng)場的光流方法對圖像質(zhì)量的要求較高,以及風(fēng)場反演中光流方程的求解方法不夠完善等原因,反演效果仍不夠理想,有待進(jìn)一步改進(jìn)[20]。而格子玻爾茲曼是基于正四邊形或正六邊形的均勻?qū)ΨQ網(wǎng)絡(luò),但因流體壓縮性和黏性等影響,在計算精度和計算效率上存在不足。
近年來,新的算法在精度和效率上有所提高但求解過程變得復(fù)雜,本文應(yīng)用Ising 模型對MRF-MAP 計算框架中的勢函數(shù)進(jìn)行了替換[21]。首先LEE T D 等[22]在20 世紀(jì)50 年代證明了Ising 模型和格子玻爾茲曼模型在力學(xué)統(tǒng)計領(lǐng)域的等價性,其次Ising 模型更適應(yīng)本文的二值分割應(yīng)用,更方便計算。
參照雷達(dá)壓縮格式完成雷達(dá)基數(shù)據(jù)sweep{SKs}的讀取之后可以得到基數(shù)據(jù)原始圖像Sk,基數(shù)據(jù)原始圖像為一張m×n大 小的灰度圖,m為雷達(dá)有效探測距離內(nèi)的徑向分辨率,n為Sk掃描一周包含的雷達(dá)射線數(shù),即Sk的雷達(dá)方向角分辨率。對基數(shù)據(jù)圖像Sk進(jìn)行PPI 模式的快速反演處理,提出反演過程如圖1 所示。
圖1 基于輪廓重建的快速反演流程Fig.1 Flowchart of rapid retrieval based on contour reconstruction
2.1.1 圖像去噪雷達(dá)接收機(jī)是一種處理回波信號的設(shè)備,可以對信號進(jìn)行變換、放大等操作,從而提取接收的微小回波信號[23]。而在雷達(dá)接收機(jī)接收信號的同時,也會收到一些對設(shè)備靈敏度產(chǎn)生負(fù)面影響的噪聲,這些噪聲影響了雷達(dá)圖像的質(zhì)量,對準(zhǔn)確發(fā)現(xiàn)和識別氣象目標(biāo)產(chǎn)生不良影響,因此雷達(dá)反演過程中通常需要去噪。在過濾噪聲的同時也會過濾掉一些圖像的細(xì)節(jié)特征,使圖像紋理變得平滑,這對后續(xù)輪廓重建的計算量會產(chǎn)生直接影響,因此以去噪平滑并保持明顯輪廓結(jié)構(gòu)為標(biāo)準(zhǔn),對比了均值濾波、小波閾值濾波和自適應(yīng)均值濾波等方法,并選取了自適應(yīng)均值濾波作為去噪處理方法[24]。
2.1.2 閾值化提取輪廓
在雷達(dá)圖像顯示中,通常將連續(xù)變化的單通道雷達(dá)灰度圖像閾值轉(zhuǎn)化成偽彩圖顯示,這在各閾值圖層之間形成了完整的封閉輪廓。將0~70 dBZ 級雷達(dá)反射圖像閾值化為4-bit 的彩色圖層Sk{Th0...Th15}, 并提取圖層輪廓。
2.1.3 坐標(biāo)變換
在PPI 模式坐標(biāo)變換中,基雷達(dá)圖像在球坐標(biāo)下轉(zhuǎn)換為笛卡爾坐標(biāo),只對各圖層輪廓線像素進(jìn)行計算,坐標(biāo)變換后的輪廓線將由連續(xù)變?yōu)椴贿B續(xù)。
圖層輪廓線在坐標(biāo)變換后從連續(xù)線變?yōu)殡x散點,為重建封閉輪廓,以輪廓線相鄰像素點為對角,將對角矩形范圍內(nèi)像素點進(jìn)行雙線性插值填充,形成隨機(jī)場。嘗試應(yīng)用MRF 分割算法來找到隨機(jī)場中兩個相鄰圖層的精確邊緣。以下闡述基于Ising 統(tǒng)計力學(xué)模型應(yīng)用MRF-MAP 算法對隨機(jī)場進(jìn)行估計。
2.2.1 定義MRF
在隨機(jī)場中的輪廓線是一組可以用馬爾可夫性描述的隨機(jī)過程,因此用MRF 方法對輪廓線兩側(cè)的格點分布進(jìn)行估計。定義S為一個二維格點位置集,s為這組位置集中的一個格點,s∈S。設(shè)x為一個隨機(jī)變量,x=0~255,將圖像X定義為一個隨機(jī)場,x∈X,xs為格點s對應(yīng)的像素值,δs記作與格點s相鄰的其他格點,并且相鄰關(guān)系符合對稱性,既s∈δγ?γ ∈δs。定義子團(tuán)(clique)C為一組具有這樣相鄰關(guān)系的格點。那么隨
在統(tǒng)計熱力學(xué)中,一般認(rèn)為能量越低系統(tǒng)越穩(wěn)定,出現(xiàn)概率越高。
2.2.2 MRF-MAP 計算框架
在完成圖層隨機(jī)場定義后,應(yīng)用MRF-MAP 框架對圖像中真實的x進(jìn)行估計。假設(shè)y為圖像觀測值,那么MAP 的求解為
在Ising 模型中粒子間的相互作用關(guān)系由他們的自旋方向決定,自旋產(chǎn)生了“↑”和“↓”的兩種狀態(tài),可以描述圖像的二值像素。通過Ising 模型的物理概念可知自旋方向相反的粒子對越多,系統(tǒng)的能量越高,越不穩(wěn)定,處于這種狀態(tài)的概率越低。模型中系統(tǒng)能量用自旋方向相反的粒子對相交邊界的邊長描述[27]。
2.2.3 Metropolis 算法
Metropolis 算法為一種固定溫度的模擬退火算法,是通過賦予搜索過程一種時變且最終趨于零的概率突跳性,從而可有效避免陷入局部極小并最終趨于全局最優(yōu)串行結(jié)構(gòu)的優(yōu)化算法[28]。應(yīng)用Metropolis 算法在圖像中計算能量函數(shù)的全局最優(yōu)值,具體步驟如下。
在完成雷達(dá)圖層輪廓線的重建后,對輪廓線內(nèi)區(qū)域進(jìn)行漫水填充賦值。將整個圖像設(shè)置為閾值顏色,然后對背景進(jìn)行填充。因為漫水填充法需要對封閉輪廓內(nèi)種子像素點標(biāo)記區(qū)域進(jìn)行填充,而雷達(dá)閾值圖層的前景多為不連續(xù)區(qū)域,所以采用背景填充。在完成每張雷達(dá)圖層后,按閾值從小到大順序?qū)⒗走_(dá)圖層集成為雷達(dá)反演圖。
硬件環(huán)境為Intel Core2 E7400 CPU 和8Gb 內(nèi)存,軟件環(huán)境為Windows7 操作系統(tǒng),VS2012 編譯和OpenCV 庫。采用美國CASA(Collaborative Adaptive Sensing of the Atmosphere)項目的X 波段多普勒雷達(dá)產(chǎn)品為數(shù)據(jù)[29]。雷達(dá)探測半徑40~60 km,具有可對近地空間觀測、安裝運(yùn)行成本低、可陸基布置也可車船安裝的優(yōu)點。雷達(dá)數(shù)據(jù)壓縮格式為“uf”,數(shù)據(jù)采集地點為某兩個雷達(dá)站(Lat:34.812 9,Lon:—97.931 3;Lat:35.031 2,Lon:—97.956 7),采 集 時 間 為2009年5 月14 日,掃描模式為仰角2°~4°的PPI 模式。雷達(dá)基數(shù)據(jù)文件包含CZ 和VE 兩個字段,即反射率字段和風(fēng)速字段,采用反射率字段為例進(jìn)行試驗。在反射率字段中,雷達(dá)的方向角分辨率為426,徑向分辨率為996,雷達(dá)基數(shù)據(jù)矩陣成像如圖2 所示。
圖2 雷達(dá)基數(shù)據(jù)文件成像Fig.2 Radar base data file imaging
首先對圖像進(jìn)行去噪處理,去噪處理不僅能提高圖像質(zhì)量,還能平滑圖像輪廓,去掉在基數(shù)據(jù)圖像中不能形成封閉輪廓的像素點,減少后續(xù)處理的計算量。分析3 種常用去噪結(jié)果,并應(yīng)用自適應(yīng)均值法,去噪效果對比如圖3 所示。
圖3 雷達(dá)基數(shù)據(jù)圖像去噪比較Fig.3 Comparison of radar cardinality according to image denoising
對處理后的基數(shù)據(jù)圖像進(jìn)行了4-bit 圖層的閾值化并提取了圖層輪廓。
在向PPI 模式的坐標(biāo)變換中,不同于其他方法,只對輪廓線像素位置進(jìn)行計算。在變換后為了將離散的輪廓線像素點連接起來,以兩個相鄰像素點為對角的矩形區(qū)域進(jìn)行了雙線性插值,使圖層重新成為封閉區(qū)域。
對連續(xù)的矩形區(qū)域應(yīng)用MAP-MRF 算法進(jìn)行精確輪廓線的計算,保證了雷達(dá)遠(yuǎn)端的輪廓線精度。對封閉輪廓區(qū)域內(nèi)進(jìn)行漫水填充如圖4 所示,并將圖層進(jìn)行合并,形成完整的PPI 模式雷達(dá)掃描圖像。
圖4 雷達(dá)圖層反演結(jié)果Fig.4 Retrieval results of radar layers
將快速雷達(dá)圖像反演結(jié)果在圖像質(zhì)量和處理時間上與其他雷達(dá)圖像重建方法進(jìn)行對比,對比方法包括最鄰近插值成像、Barnes 分析成像及雙線性插值成像,結(jié)果如圖5 所示。
圖5 雷達(dá)圖像反演結(jié)果Fig.5 Radar image retrieval results
在成像質(zhì)量的對比分析上,通常依據(jù)人的主觀感受及一些客觀評價指標(biāo),如峰值信噪比(PSNR)和均方誤差(MSE)等。但因為快速反演方法的圖像區(qū)域大部分為填充區(qū)域,所以為了客觀評價反演效果,采用直接驗證法,以單束雷達(dá)基數(shù)據(jù)作為真值進(jìn)行分析。在對比分析中,采樣了兩個雷達(dá)站掃描空域的交集區(qū)域,以閾值化后的一個雷達(dá)站基數(shù)據(jù)為真值,對另一個雷達(dá)站雷達(dá)波束上的400 個像素進(jìn)行了采樣統(tǒng)計。分析不同重建方法處理后的圖像中像素能夠被正確分類的比例,結(jié)果如表1 所示[30-32]。在重建方法對比中,準(zhǔn)確率最高的是Barnes 分析法,本文方法優(yōu)于最鄰近插值法,在延單一雷達(dá)照射波束上的準(zhǔn)確率統(tǒng)計與其他兩種方法差距不大。
表1 數(shù)據(jù)重建準(zhǔn)確率對比Tab.1 Data reconstruction accuracy comparison
在算法效率對比中,為了便于分析,對雷達(dá)基數(shù)據(jù)圖像進(jìn)行像素級的插值重建處理?;鶖?shù)據(jù)圖像大小為996 像素×426 像素,重建圖像為2 000 像素×2 000像素,重建圖像像素數(shù)量約為原始圖像9 倍。對10 組基數(shù)據(jù)文件進(jìn)行了處理,測試樣本數(shù)據(jù)閾值化后圖層為8~13 層。在未對代碼進(jìn)行優(yōu)化的情況下測得運(yùn)行時間如圖6 所示。
從運(yùn)行時間來看,基于輪廓重建的MAP-MRF 方法能大大縮短雷達(dá)圖像的反演時間。但由圖6 可知,與其他方法不同,本文方法測試樣本的運(yùn)行時間方差較大。主要影響因素是雷達(dá)基數(shù)據(jù)的范圍和圖層中輪廓線劃分的隨機(jī)場大小。通過試驗得出劃分閾值后的圖層數(shù)量對算法處理時間影響較大。輪廓線像素數(shù)量影響劃分隨機(jī)場的大小,進(jìn)而圖像需要的迭代時長不同。
圖6 文中方法與傳統(tǒng)方法運(yùn)行時間對比Fig.6 Processing time of proposed method compared with traditional methods
海洋捕撈業(yè)一直有“起風(fēng)魚”之說,即在壞天氣系統(tǒng)如臺風(fēng)侵入前,魚群往往大而集中并浮表層,因為上述天氣系統(tǒng)將帶來大風(fēng),漁民抓住這一時機(jī)進(jìn)行“迎風(fēng)頭”生產(chǎn),魚獲量常常較大。但是此時作業(yè)常常會伴隨大風(fēng)大浪,重大漁船海損事故幾乎都是在此氣象條件下發(fā)生的,因此準(zhǔn)確掌握氣象變化顯得尤為重要。針對漁船類小型船舶上雷達(dá)探測顯示時間間隔長,應(yīng)對海上極端氣象目標(biāo)時間短的問題,提出了基于輪廓重建的MAP-MRF 快速雷達(dá)圖像反演方法。應(yīng)用插值過程針對輪廓線計算量較少的原理,闡述了Ising 模型對大氣流體運(yùn)動建模的可行性,引入該模型對輪廓線的分布進(jìn)行描述,替換了傳統(tǒng)的格子氣模型對介觀尺度大氣流體運(yùn)動的描述。通過MAP-MRF 算法進(jìn)行輪廓線內(nèi)外區(qū)域的分割,以及CASA 試驗場真實數(shù)據(jù),驗證了該方法的插值準(zhǔn)確性。與氣象和圖像插值領(lǐng)域主要方法進(jìn)行對比,論證了針對輪廓進(jìn)行快速反演插值的方法具有更快的處理速度,有利于漁船雷達(dá)系統(tǒng)達(dá)到更高的刷新率,實現(xiàn)對海上極端氣象目標(biāo)的及時發(fā)現(xiàn)和連續(xù)追蹤。通過試驗分析,該方法對于雷達(dá)圖像的反演速度有較大提升,優(yōu)化效果受雷達(dá)探測對象的空間分布情況影響較大。