基于輪廓重建的漁船雷達MAP-MRF 快速反演研究

鄭建麗，曹建軍，劉 洋

（1.中國水產科學研究院漁業機械儀器研究所，上海 200092；2.遼寧師范大學計算機與信息技術學院，遼寧 大連 116081）

0 引 言

船載雷達是放置在船上用來觀測船舶周圍各種水面目標和即時氣象環境的重要觀測工具。在船舶航行中遇到霧、霾、雨等現象導致能見距離受到影響，無法憑目視發現目標的條件下，雷達表現出受環境影響小、觀測距離遠的優勢。目前，常用的航海雷達有SPERRY-MK 系 列、JRC-JMA 系 列 和FURUNO-FAR系列等。根據國際電工委員會IEC62388—2013 標準，船載雷達的顯示圖像應能夠進行連續且平滑的雷達圖像更新[1]。但目前主流的雷達終端只能滿足部分觀測需求，還無法達到連續平滑的顯示效果，難以觀測海上近地大氣層一些快速生成和移動的強對流氣象目標。從觀測對象的角度分析，海洋環境水汽熱量充足，下墊面平滑，氣象目標形成和移動速度較快。約70%的龍卷風從出現到消失的生命周期為1～10 min，海上冷風的移動速度可達70 km/h 以上[2]。這些快速生成和變化的氣象目標對雷達數據處理提出了高清顯示和更快速度的要求。從觀測平臺的角度分析，船載雷達因平臺穩定性差和能量系統相對封閉等原因，探測距離一般為40～80 km。相較于陸基雷達通常200～300 km 的探測距離，船載雷達的目標探測距離更短，因而有效應對時間更少。

陸基雷達處理的數據量通常遠遠高于船載雷達，但傳統陸基雷達對于提升數據處理速度的解決方案不適用于船載雷達。通常陸基雷達產品生成一般分為前端和后端兩部分，前端計算機用來控制雷達及其機構并收集反射信號，后端計算機用來將雷達信號反演成圖像產品。在陸基雷達的反演過程中后端計算機通常為計算能力較強的超級計算機，通過改善硬件的方式提高反演速度。但由于小型船舶制造成本和能量環境等限制因素，應用改進硬件來提高反演速度的方法并不適用于目前的漁船和游艇等小型船舶的情況。

目前，國內海洋漁船一般采用基于（8～12 GHz）波段的小型雷達，通過天線發射電磁波到海面，電磁波與雷達波長相當的毛細波產生Bragg 散射，后向散射回波后被天線接收，形成“海雜波”，接收終端把連續接收的回波形成圖像序列，而波長較長的重力波通過對毛細波的調制作用表現在海雜波圖像上，從中提取風、浪、流等海洋環境信息[3-5]。本文分析傳統雷達圖像反演過程和成像步驟，利用對雷達圖層封閉輪廓的重建，提出基于輪廓重建的船載雷達MAP-MRF 快速反演方法，解決漁船雷達圖像反演過程中大量重復計算缺失數據的重建插值過程。通過統計力學模型對MRF 中格點劃分進行精確建模，利用MAP-MRF 方法對基于輪廓重建的雷達圖層邊界進行估計。將整幅圖像的插值重建轉化為對圖層輪廓的精確分割，實現雷達圖像的快速反演，從而提高顯示終端的刷新率，實現對海上極端氣象目標的及時發現和連續追蹤，滿足漁民對天氣氣象的變化情況特別是壞氣象實時了解的需求。在證明基于輪廓反演方法有效性的同時，驗證應用基于統計力學模型的MAP-MRF 算法對雷達圖像進行相位估計的可行性，為進一步從概率模型方向處理雷達圖像提供新的分析思路和處理方法。主要內容分為3 個部分：一是通過對漁船雷達數據和反演過程的分析，明確反演方法改進重點和改進內容；二是對PPI 模式下的雷達基數據圖像應用MAP-MRF 算法重建封閉輪廓進行精確重建，進而反演成像；三是將試驗結果與較有代表性的傳統反演方法在成像質量和運行速度上進行分析和比較。

1 研究內容

1.1 雷達圖像反演

雷達圖像的重建可分為基于時間維度和基于空間維度的兩種過程，二者分別對應雷達反演過程中時間分辨率（temporal resolution）和空間分辨率（spatial resolution）。本文主要討論單幀雷達圖像的空間超分辨率重建方法。

雷達的掃描方式是建立在球坐標系中，但其展示和應用通常在笛卡爾坐標系下，因此雷達圖像的反演處理中包含了兩個坐標系下的數據坐標轉換。當坐標系變換后，雷達的輻射結構會造成近雷達端數據采樣密集混疊，遠雷達端數據采樣稀疏缺失的情況。因此，雷達通常會對坐標系變換后的雷達圖像進行重建或插值操作來改善以上問題。將雷達圖像重建算法分為雷達格點分析法和通用插值算法兩類。

雷達格點分析法是在雷達和氣象領域常用的數據插值和分析方法，如Barnes 分析法和Cressman 分析法等[6-9]。Barnes 分析法能夠把分布不均勻的數據映射到規則網格坐標系的指數權重，在映射到笛卡爾坐標系時通常應用均一Barnes 插值，在插值到球坐標系時一般應用自適應Barnes 插值。重建雷達圖像具有平滑、失真較小的優點，此外Barnes 分析法還可用于基于站點序列的插值和糾錯。Cressman 分析法用于從基數據向格點空間的映射，通過加權影響半徑內基數據對格點值進行計算，采用迭代影響半徑的方式優化精度依據影響半徑。Cressman 計算速度快且重建準確率高，但當映射空間密度高于基數據密度時算法穩定性較弱，邊緣附近易出現極值點。

通用插值算法因應用范圍廣，在性能和比較分析方面有比較深入的研究[10-14]。其中最常用的算法包括最鄰近插值法、雙線性插值法、基于小波的插值法，以及近年來被廣泛研究的基于機器學習理論的插值法等。其中最鄰近插值法是最簡單也是最快速的一種方法，計算復雜度最低，對于噪聲敏感。雙線性插值法是在實際應用中被廣泛采用的算法，時間復雜度較低，峰值信噪比、均方差及主觀評價等較最鄰近插值法有明顯改善?；谛〔ǖ牟逯捣ㄊ抢眯〔ǚ纸夂蟮母哳l信息預測高分辨率圖像，重構高分辨率圖像的方法，是應用變換域插值的一種主要方法。機器學習的圖像超分辨率重建法屬于基于樣例的重建，其中比較有代表性的如DONG W 等[15]基于稀疏編碼重建，以及SHI W 等[16]對卷積網絡的深度訓練圖像重建等，對于普通圖像深度學習方法具有較好的重建效果，但在氣象雷達領域面臨著探測目標通常不具有穩定的結構特征，重建內容缺少對應訓練數據的問題。

在圖像重建的研究中，基于邊緣信息的重建是比較有代表性的方向。圖像邊緣表達了更加豐富的圖像信息，通過解決重建過程中的邊緣模糊效應和鋸齒狀不連續等問題能夠有效改善整幅重建效果[17-18]。因此有效利用圖像邊緣信息是有效改善重建效果、提升重建效率的一個有效方向。受此啟發，提出精確分割閾值圖層輪廓實現雷達圖像重建的方法。同時，漁船雷達圖像是灰度圖像，圖像中的艦船目標和陸地區域灰度級較高，海面區域灰度值較低，目標與目標之間不會重疊，對此圖像重建方法具有一定的自適應性。

1.2 空氣運動描述模型

雷達圖像描述的輪廓邊緣既對應物理特征的等高線，又描述氣象目標內部降水粒子的尺度和密度分布。這些粒子懸浮在大氣中，其分布和運動特征可以用大氣流體運動模型描述。

此類模型主要是指基于納維—斯托克斯（Navier-Stokes，NS）方程的各類多相流模型，模擬方法以傳統計算流體力學方法為主，光滑粒子方法、格子玻爾茲曼方法（Lattice Boltzmann Method，LBM）應用比較廣泛[19]。

2004 年，Dankert 從雷達圖像序列中反演海面風場，研究了以光流運動估計技術為基礎的風場反演方法。此方法不需要校正雷達系統的相位，具有較強的實用價值。然而，由于反演風場的光流方法對圖像質量的要求較高，以及風場反演中光流方程的求解方法不夠完善等原因，反演效果仍不夠理想，有待進一步改進[20]。而格子玻爾茲曼是基于正四邊形或正六邊形的均勻對稱網絡，但因流體壓縮性和黏性等影響，在計算精度和計算效率上存在不足。

近年來，新的算法在精度和效率上有所提高但求解過程變得復雜，本文應用Ising 模型對MRF-MAP 計算框架中的勢函數進行了替換[21]。首先LEE T D 等[22]在20 世紀50 年代證明了Ising 模型和格子玻爾茲曼模型在力學統計領域的等價性，其次Ising 模型更適應本文的二值分割應用，更方便計算。

2 輪廓重建快速反演算法

參照雷達壓縮格式完成雷達基數據sweep{SKs}的讀取之后可以得到基數據原始圖像Sk，基數據原始圖像為一張m×n大 小的灰度圖，m為雷達有效探測距離內的徑向分辨率，n為Sk掃描一周包含的雷達射線數，即Sk的雷達方向角分辨率。對基數據圖像Sk進行PPI 模式的快速反演處理，提出反演過程如圖1 所示。

圖1 基于輪廓重建的快速反演流程Fig.1 Flowchart of rapid retrieval based on contour reconstruction

2.1 雷達圖像預處理

2.1.1 圖像去噪雷達接收機是一種處理回波信號的設備，可以對信號進行變換、放大等操作，從而提取接收的微小回波信號[23]。而在雷達接收機接收信號的同時，也會收到一些對設備靈敏度產生負面影響的噪聲，這些噪聲影響了雷達圖像的質量，對準確發現和識別氣象目標產生不良影響，因此雷達反演過程中通常需要去噪。在過濾噪聲的同時也會過濾掉一些圖像的細節特征，使圖像紋理變得平滑，這對后續輪廓重建的計算量會產生直接影響，因此以去噪平滑并保持明顯輪廓結構為標準，對比了均值濾波、小波閾值濾波和自適應均值濾波等方法，并選取了自適應均值濾波作為去噪處理方法[24]。

2.1.2 閾值化提取輪廓

在雷達圖像顯示中，通常將連續變化的單通道雷達灰度圖像閾值轉化成偽彩圖顯示，這在各閾值圖層之間形成了完整的封閉輪廓。將0～70 dBZ 級雷達反射圖像閾值化為4-bit 的彩色圖層Sk{Th0...Th15}, 并提取圖層輪廓。

2.1.3 坐標變換

在PPI 模式坐標變換中，基雷達圖像在球坐標下轉換為笛卡爾坐標，只對各圖層輪廓線像素進行計算，坐標變換后的輪廓線將由連續變為不連續。

2.2 基于MAP-MRF 的封閉輪廓重建

圖層輪廓線在坐標變換后從連續線變為離散點，為重建封閉輪廓，以輪廓線相鄰像素點為對角，將對角矩形范圍內像素點進行雙線性插值填充，形成隨機場。嘗試應用MRF 分割算法來找到隨機場中兩個相鄰圖層的精確邊緣。以下闡述基于Ising 統計力學模型應用MRF-MAP 算法對隨機場進行估計。

2.2.1 定義MRF

在隨機場中的輪廓線是一組可以用馬爾可夫性描述的隨機過程，因此用MRF 方法對輪廓線兩側的格點分布進行估計。定義S為一個二維格點位置集，s為這組位置集中的一個格點，s∈S。設x為一個隨機變量，x=0～255，將圖像X定義為一個隨機場，x∈X，xs為格點s對應的像素值，δs記作與格點s相鄰的其他格點，并且相鄰關系符合對稱性，既s∈δγ?γ ∈δs。定義子團（clique）C為一組具有這樣相鄰關系的格點。那么隨

在統計熱力學中，一般認為能量越低系統越穩定，出現概率越高。

2.2.2 MRF-MAP 計算框架

在完成圖層隨機場定義后，應用MRF-MAP 框架對圖像中真實的x進行估計。假設y為圖像觀測值，那么MAP 的求解為

在Ising 模型中粒子間的相互作用關系由他們的自旋方向決定，自旋產生了“↑”和“↓”的兩種狀態，可以描述圖像的二值像素。通過Ising 模型的物理概念可知自旋方向相反的粒子對越多，系統的能量越高，越不穩定，處于這種狀態的概率越低。模型中系統能量用自旋方向相反的粒子對相交邊界的邊長描述[27]。

2.2.3 Metropolis 算法

Metropolis 算法為一種固定溫度的模擬退火算法，是通過賦予搜索過程一種時變且最終趨于零的概率突跳性，從而可有效避免陷入局部極小并最終趨于全局最優串行結構的優化算法[28]。應用Metropolis 算法在圖像中計算能量函數的全局最優值，具體步驟如下。

2.3 閾值圖層集成

在完成雷達圖層輪廓線的重建后，對輪廓線內區域進行漫水填充賦值。將整個圖像設置為閾值顏色，然后對背景進行填充。因為漫水填充法需要對封閉輪廓內種子像素點標記區域進行填充，而雷達閾值圖層的前景多為不連續區域，所以采用背景填充。在完成每張雷達圖層后，按閾值從小到大順序將雷達圖層集成為雷達反演圖。

3 反演試驗和分析

3.1 試驗數據和仿真環境

硬件環境為Intel Core2 E7400 CPU 和8Gb 內存，軟件環境為Windows7 操作系統，VS2012 編譯和OpenCV 庫。采用美國CASA（Collaborative Adaptive Sensing of the Atmosphere）項目的X 波段多普勒雷達產品為數據[29]。雷達探測半徑40～60 km，具有可對近地空間觀測、安裝運行成本低、可陸基布置也可車船安裝的優點。雷達數據壓縮格式為“uf”，數據采集地點為某兩個雷達站（Lat：34.812 9，Lon：—97.931 3；Lat：35.031 2，Lon：—97.956 7），采 集 時 間 為2009年5 月14 日，掃描模式為仰角2°～4°的PPI 模式。雷達基數據文件包含CZ 和VE 兩個字段，即反射率字段和風速字段，采用反射率字段為例進行試驗。在反射率字段中，雷達的方向角分辨率為426，徑向分辨率為996，雷達基數據矩陣成像如圖2 所示。

圖2 雷達基數據文件成像Fig.2 Radar base data file imaging

3.2 雷達反演過程

首先對圖像進行去噪處理，去噪處理不僅能提高圖像質量，還能平滑圖像輪廓，去掉在基數據圖像中不能形成封閉輪廓的像素點，減少后續處理的計算量。分析3 種常用去噪結果，并應用自適應均值法，去噪效果對比如圖3 所示。

圖3 雷達基數據圖像去噪比較Fig.3 Comparison of radar cardinality according to image denoising

對處理后的基數據圖像進行了4-bit 圖層的閾值化并提取了圖層輪廓。

在向PPI 模式的坐標變換中，不同于其他方法，只對輪廓線像素位置進行計算。在變換后為了將離散的輪廓線像素點連接起來，以兩個相鄰像素點為對角的矩形區域進行了雙線性插值，使圖層重新成為封閉區域。

對連續的矩形區域應用MAP-MRF 算法進行精確輪廓線的計算，保證了雷達遠端的輪廓線精度。對封閉輪廓區域內進行漫水填充如圖4 所示，并將圖層進行合并，形成完整的PPI 模式雷達掃描圖像。

圖4 雷達圖層反演結果Fig.4 Retrieval results of radar layers

3.3 試驗結果對比分析

將快速雷達圖像反演結果在圖像質量和處理時間上與其他雷達圖像重建方法進行對比，對比方法包括最鄰近插值成像、Barnes 分析成像及雙線性插值成像，結果如圖5 所示。

圖5 雷達圖像反演結果Fig.5 Radar image retrieval results

在成像質量的對比分析上，通常依據人的主觀感受及一些客觀評價指標，如峰值信噪比（PSNR）和均方誤差（MSE）等。但因為快速反演方法的圖像區域大部分為填充區域，所以為了客觀評價反演效果，采用直接驗證法，以單束雷達基數據作為真值進行分析。在對比分析中，采樣了兩個雷達站掃描空域的交集區域，以閾值化后的一個雷達站基數據為真值，對另一個雷達站雷達波束上的400 個像素進行了采樣統計。分析不同重建方法處理后的圖像中像素能夠被正確分類的比例，結果如表1 所示[30-32]。在重建方法對比中，準確率最高的是Barnes 分析法，本文方法優于最鄰近插值法，在延單一雷達照射波束上的準確率統計與其他兩種方法差距不大。

表1 數據重建準確率對比Tab.1 Data reconstruction accuracy comparison

在算法效率對比中，為了便于分析，對雷達基數據圖像進行像素級的插值重建處理?；鶖祿D像大小為996 像素×426 像素，重建圖像為2 000 像素×2 000像素，重建圖像像素數量約為原始圖像9 倍。對10 組基數據文件進行了處理，測試樣本數據閾值化后圖層為8～13 層。在未對代碼進行優化的情況下測得運行時間如圖6 所示。

從運行時間來看，基于輪廓重建的MAP-MRF 方法能大大縮短雷達圖像的反演時間。但由圖6 可知，與其他方法不同，本文方法測試樣本的運行時間方差較大。主要影響因素是雷達基數據的范圍和圖層中輪廓線劃分的隨機場大小。通過試驗得出劃分閾值后的圖層數量對算法處理時間影響較大。輪廓線像素數量影響劃分隨機場的大小，進而圖像需要的迭代時長不同。

圖6 文中方法與傳統方法運行時間對比Fig.6 Processing time of proposed method compared with traditional methods

4 結束語

海洋捕撈業一直有“起風魚”之說，即在壞天氣系統如臺風侵入前，魚群往往大而集中并浮表層，因為上述天氣系統將帶來大風，漁民抓住這一時機進行“迎風頭”生產，魚獲量常常較大。但是此時作業常常會伴隨大風大浪，重大漁船海損事故幾乎都是在此氣象條件下發生的，因此準確掌握氣象變化顯得尤為重要。針對漁船類小型船舶上雷達探測顯示時間間隔長，應對海上極端氣象目標時間短的問題，提出了基于輪廓重建的MAP-MRF 快速雷達圖像反演方法。應用插值過程針對輪廓線計算量較少的原理，闡述了Ising 模型對大氣流體運動建模的可行性，引入該模型對輪廓線的分布進行描述，替換了傳統的格子氣模型對介觀尺度大氣流體運動的描述。通過MAP-MRF 算法進行輪廓線內外區域的分割，以及CASA 試驗場真實數據，驗證了該方法的插值準確性。與氣象和圖像插值領域主要方法進行對比，論證了針對輪廓進行快速反演插值的方法具有更快的處理速度，有利于漁船雷達系統達到更高的刷新率，實現對海上極端氣象目標的及時發現和連續追蹤。通過試驗分析，該方法對于雷達圖像的反演速度有較大提升，優化效果受雷達探測對象的空間分布情況影響較大。