一種新穎的亞像素級SAR圖像水陸分割方法

毋 琳，牛世林，郭拯危,3,4，李 寧，3,4

(1.河南大學計算機與信息工程學院，河南開封 475004；2.河南大學環境與規劃學院，河南開封 475004；3.河南省智能技術與應用工程技術研究中心，河南開封 475004；4.河南省大數據分析與處理重點實驗室，河南開封 475004)

0 引言

地表水是地球水資源的重要組成部分，在地球的水文和生物化學循環中具有舉足輕重的作用。利用遙感圖像實現精確的水陸分割是開展基于水體各項研究的工作基礎，有著重要的現實意義。合成孔徑雷達(Synthetic Aperture Radar, SAR)是一種主動微波遙感技術，可穿透云層、降雨、霧霾等實現全天時、全天候對地觀測，且具有響應速度快、覆蓋范圍廣、周期性觀測等特點，目前已廣泛應用于大尺度地表調查與研究中[1]。

SAR圖像水陸分割研究發展至今，經典的閾值、聚類、主動輪廓模型(Active Contour Model，ACM)等方法，以及新興的深度學習算法都有較多的研究與應用[2-6]。這些方法雖然能夠有效地實現SAR圖像水陸分割，但分割尺度多停留在像素級別。對于星載SAR系統而言，由于受最小天線孔徑積的限制，常見工作模式的圖像空間分辨率通常在幾米至幾十米量級[7]。例如Sentinel-1衛星TOPS模式下的空間分辨率為20 m，國產高分三號(GF-3)衛星標準條帶模式下的空間分辨率為8 m。這些中低分辨率的圖像中，水域邊界上會存在大量既包含陸地又包含水域的混合象元，傳統算法無法在混合象元內確定真實水域輪廓位置，只能將整個象元劃分為陸地或者水域，極大影響了水陸分割結果的精度。同時，由于SAR圖像固有的相干斑噪聲在水域邊界附近分布雜亂，會破壞真實的水域邊界，這也是導致水陸分割精度難以提高的原因之一。而SAR圖像水陸分割中幾個像素的誤差對應于實際地物往往是數十米甚至上百米的差異。因此，設計一套性能良好且可實現亞像素級水陸分割的方法，對于更加有效地利用中低分辨率SAR圖像具有重要意義。

針對上述問題，本文提出了一種新穎的亞像素級SAR圖像水陸分割方法，該方法以改進非局域相干斑濾波算法與模糊C均值(Fuzzy C-Mean,FCM)聚類算法完成像素級水陸粗分割，在此基礎上根據水域特征，消除FCM聚類粗分割結果中的噪聲、陰影等干擾，利用雙三次樣條插值算法對獲得的像素級水域輪廓鄰域進行上采樣，并在該鄰域內應用幾何主動輪廓(Geometric Active Contour, GAC)模型獲取亞像素級水陸分割位置點，最終通過分段曲線擬合方法獲得完整的亞像素級水域輪廓。為了驗證本文所提方法的有效性，應用國產GF-3衛星的多模式SAR圖像進行驗證實驗，結果表明，與傳統方法相比，本文所提方法平均像素偏移精度提升了一個數量級，具有良好的水陸分割效果。

1 亞像素級水陸分割方法

本文提出的亞像素級SAR圖像水陸分割方法流程如圖1所示。該方法的第一步，使用改進的非局域濾波算法抑制圖像中的相干斑噪聲；第二步，應用FCM算法對水陸邊界進行像素級的粗分割；第三步，結合雙三次樣條插值和GAC模型，在像素級水域輪廓的基礎上實現亞像素級水陸分割，通過以上步驟可以獲取到更精確的亞像素級水域輪廓。

圖1 本文方法流程圖

1.1 非局域相干斑濾波算法

相干斑是一種SAR圖像特有的乘性噪聲，當電磁波照射粗糙表面時，電磁波到達的相位不同，導致回波產生干涉現象，使得SAR圖像中分布著大量由干涉引起的斑點噪聲，這些噪聲會使SAR圖像處理的難度大大增加。近年來，在抑制SAR圖像相干斑噪聲方面，非局域濾波算法取得了較好的效果。

非局域濾波算法中采用像素相關性反映圖像相鄰位置像素值的相關程度。文獻[8]證明了比值距離(Ratio Distance, RD)是一種魯棒性強的特征，可以用來描述SAR圖像中像素之間的相關性。本節將介紹一種基于RD的非局域SAR圖像濾波算法。

SAR圖像中每個像素的測量值Y，真實值X以及噪聲值Z三者的關系可以被描述為Y=X×Z。對于一景L視SAR幅度圖像，測量值通常用Goodman模型中的幅度模型描述，該模型又被稱為均方根伽馬(Nagakami-Rayleigh)分布[8]，如式(1)所示。

A≥0

(1)

式中，A為每個像素的觀測值，σ為待估計的真實強度值，則真實的幅度值X為σ的均方根。

非局域濾波算法中，定義以每個像素為中心的小正方形窗為一個相似塊。假設兩個相似塊的真實幅度值相同，觀測值的RD將從YNi/YNj轉變為ZNi/ZNj。將新的RD帶入到塊相似性測量方程S(i,j)以及兩個相似塊的聯合概率密度函數p(i,j)可得式(2)和式(3)。

(2)

(3)

式中：j是非局域窗中的像素，非局域窗是以像素i為中心、邊長為r1的正方形窗；G為標準高斯核函數；ZNk是相似塊Nk的觀測值之和，Nk是以像素k為中心、邊長為r2的正方形相似塊；p(ZNi/ZNj)是ZNi/ZNj的概率密度函數。

濾波后的SAR圖像I可由新定義的S(i,j)計算得到，如式(4)所示。

(4)

1.2 基于FCM算法的像素級水陸分割方法

FCM是一種無監督的聚類算法。一幅圖像可描述為集合I={i1,i2,…,iN}，傳統FCM算法通過設置c個不同的聚類中心，計算集合I中每個像素與各聚類中心的隸屬度，完成對圖像的分類。FCM算法執行時，需遍歷圖像中每一個像素計算隸屬度，當處理大尺度SAR圖像時，非常耗時。針對該問題，文獻[4]提出了一種基于灰度分布的FCM算法，采用如式(5)所示的目標函數Jm，將與聚類中心隸屬度的計算由像素轉變為灰度級，在灰度級上完成對圖像的分類，從而大大提高了算法的運行效率。

(5)

(6)

(7)

本文用于提取像素級水陸分割的FCM算法執行步驟如下：

(1)設置聚類中心數c=4，并隨機設置聚類中心數值；設置模糊指數m=2；設置迭代次數m1=15。

(2)根據式(6)計算樣本隸屬度，即SAR圖像中各灰度級到各聚類中心的隸屬度。

(3)根據式(7)更新各聚類中心。

(4)判斷迭代是否完成，未完成則循環執行步驟(2),步驟(3)。

(5)迭代完成后，若‖Ii-vj‖2<‖Ii-vk‖2,(j,k)∈[1,c]，則將灰度值等于i的像素劃分為第j類。

1.3 亞像素級水陸分割方法

本節提出了一種新穎的亞像素級水陸分割方法，該方案以像素級水陸分割結果為基礎，利用雙三次樣條插值算法提高其周邊區域的清晰度，并結合GAC模型求解梯度最大位置來確定亞像素級水陸邊界位置，方法流程如圖1中步驟3所示。

該方法是在以像素級水域邊界線的像素點為中心，半徑為r3的鄰域內進行水陸精細分割的，其有效性建立在兩個假設之上：第一，鄰域內需同時包含水域與陸地，如圖2(a)所示；第二，當雷達后向散射系數梯度值最大時，鄰域內的陸地與水域產生分界[9]，如圖2(b)所示。

(a)沿像素級水域輪廓建立SPW (b)SPW內SAR圖像 (c)SPW內每個像素灰度值 (d)經過插值后的SPW以及提取到的亞像素級水陸分割線

亞像素級水陸分割過程包括以下5個步驟：

(1)自然水域周邊地形一般較為復雜，陡峭山脈等地物會因起伏地形未受到電磁波照射，沒有回波，會產生較多的陰影現象，在SAR圖像中呈現為黑色區域。而水面對電磁波產生表面散射，后向散射系數較低，平靜水域在SAR圖像中也同樣呈現為與陰影相似的黑色區域。為實現對二者的區分，需進一步提取感興趣區域(Region of Interest, ROI)。在粗分割獲取的像素級水域中，記最大連通區域面積為Smax，設置面積閾值Ta=μ·Smax。當連通區域面積大于Ta時，被記為ROI，否則丟棄。同時，將提取到的ROI區域邊界信息記為像素級水域輪廓。本文根據經驗將μ設置為0.2。

(2)以步長k遍歷像素級水域輪廓，以訪問到的每一個像素為中心，建立邊長為r3的正方形亞像素處理窗(Subpixel Processing Window，SPW)，記錄該SPW中心坐標為(x1,y1)，計算得到第一個像素坐標為(x2,y2)=(x1-r3,y1-r3)，如圖2(a)中虛線框所示。在一個SPW中，以第一個像素為坐標原點，建立坐標系，則該SPW中所包含的像素級水域輪廓點可通過多項式擬合，表示為曲線方程：fw l(x,y)。方程階數取決于SPW的大小，SPW越大，包含的像素級水域輪廓點越多，多項式階數越高。通常，多項式階數過高可能會導致最終結果不具備魯棒性，相反小的SPW則會對高階多項式要求較小。

(3)利用雙三次樣條插值方法，對SPW內圖像進行λ倍的上采樣，增加局部像素數目，減小后續在亞像素尺度下調整水陸邊界時對最終邊界的影響，如圖2(c)與圖2(d)所示。

圖3 完整的亞像素級水陸分割線獲取示意圖

1.4 基于GAC模型的亞像素級水陸分割線獲取方法

GAC模型是基于Snake模型和水平集函數(Level-Set Function，LSF)發展而來的，它將能量最小化問題轉化為在黎曼空間中求最小曲線路徑的問題[5]。GAC模型對應的LSF可以表示為

(8)

式中，?φ為初始LSF，為梯度算子，g(x)為基于梯度信息的邊緣停止函數，α為控制輪廓膨脹或收縮的氣球力參數。GAC算法的執行步驟如下：

(1)按如式(9)所示方法初始化LSF。

(9)

式中，ψ為SPW內的圖像，ψ0為SPW中的水域區域，σψ0為由SPW中擬合的fwl函數獲得的曲線，ρ為常數參數，本文設置ρ=1。

(2)設置氣球力參數α，并根據式(8)演化LSF。本文令α=5。

(3)使用高斯濾波器調整LSF，如φ=φ×Gτ。τ為高斯函數的標準差，范圍為0.8到1.5。本文設置高斯核的邊長h=5，標準差τ=1。

(4)設置最大演化次數m2=30，若未達到演化次數，返回步驟(2)。

2 實驗結果與分析

2.1 研究區域及數據

實驗數據采用國產GF-3號衛星在丹江口水庫主庫區觀測得到的兩景SAR圖像，范圍如圖4中黑色框所示，圖像詳細信息如表1所示。國產GF-3號衛星是一顆C波段全極化SAR衛星，適用于多種工作場景，包括海洋監測、災害監測等[10]。

表1 實驗使用的GF-3號SAR圖像詳細參數

圖4 丹江口水庫位置

2.2 實驗結果

除本文上述步驟中已給出的實驗參數外，本文方法其他的實驗參數如下：(1)非局域相干斑濾波步驟中，非局域窗邊長r1=15，相似塊邊長r2=7，視數L=1。(2)亞像素水陸分割線提取方案中，步長k=2，SPW邊長r3=7，上采樣倍數λ=5。

圖5展示了非局域濾波算法對SAR圖像相干斑噪聲抑制的結果。兩景原始SAR圖像如圖5(a)和圖5(f)所示。為了更好地展示細節，選取區域1至區域4放大，原始圖像如圖5(b)、(d)、(g)和(i)所示，濾波后圖像如圖5(c)、(e)、(h)和(j)所示。由圖示結果可知，該濾波方法具有的良好平滑性能和邊緣保持效果，陸地與水域上的斑點噪聲被平滑，陸地呈均勻的灰色，水域呈均勻的黑色，水域邊界也得到了保持和加強，整體上提高了SAR圖像的均勻性，為后續處理創造了良好的條件。

圖5 相干斑濾波結果

圖6展示了像素級水陸分割線提取的結果和亞像素級水陸分割線提取的結果。為了更為清晰地觀察，選擇區域A到D放大，分割結果如圖6(b)、(c)、(e)和(f)所示。紅線代表亞像素級水域輪廓線，藍線代表像素級水域輪廓線。兩景SAR圖像的標稱分辨率分別為5 m和10 m，因此兩景SAR圖像中的一個像素分別代表25 m2和100 m2。結果的目視分析可知，亞像素級水域輪廓(紅線)比像素級(藍線)更加平滑，更加接近真實水陸邊界。因此，在水域SAR遙感監測中，亞像素級水陸分割具有更好的性能與應用前景。

圖6 亞像素級水陸分割線提取結果

2.3 對比實驗與定量分析

為了驗證本文所提方法的有效性，選取目前主流，且精度較高的馬爾可夫隨機場(Markov Random Field，MRF)方法[3]、FCM方法[4]和ACM方法[5]這三種SAR圖像水陸分割方法開展對比實驗，分割結果如圖7所示。同時，選擇虛警率(False Alarm Rate)、準確率(ACC)以及平均輪廓偏移距離(Average Contour Offset Distance, ACOD)等量化指標分別評估整體分割效果與局部輪廓質量[3]。參考水域為人工追蹤法[2]獲取的水域輪廓。評價指標計算方法如式(10)、式(11)和式(12)所示。

圖7 不同方法水域提取結果

(1)虛警率

(10)

(2)準確率

(11)

(3)平均輪廓偏移距離

(12)

式中，AE為算法檢測的水域面積，AR為人工檢測的參考水域面積，AE&R為AE與AR的交集面積，N是提取水域輪廓上的像素數，d(EinRBi;R)是像素i與實際水域輪廓之間的最短距離。

目視觀察各方法的實驗結果，本文所提方法與MRF和FCM方法相比，對于水域的提取更為精確，受SAR圖像陰影、相干斑噪聲以及水面復雜散射影響較少。ACM算法在大尺度SAR圖像中演化較慢，當達到最大演化次數時，仍未檢測出水庫邊緣的細小分支。該結果也可由如表2所示的定量分析結果驗證，本文方法對兩組實驗數據的提取準確率均大于99%，虛警率均小于0.5%，并且輪廓局部質量優異，輪廓平均偏移距離均小于0.2個像素，取得了更高精度的結果，優于其他方法。

表2 定量分析結果

3 結束語

本文針對中低分辨率SAR圖像精確水陸分割的難題，提出了一種新穎的亞像素級SAR圖像水陸分割方法。該方法利用改進的非局域濾波算法對SAR圖像的相干斑噪聲進行抑制；采用FCM算法完成像素級水陸粗分割；并在粗分割基礎上，應用基于雙三次樣條插值算法與GAC模型的方法完成亞像素級的水陸精分割。通過對丹江口水庫區域的兩組GF-3號衛星多模式SAR圖像的實驗驗證，與目前常見SAR圖像水陸分割算法相比，本文方法取得了更高的準確率以及更低的虛警率，平均像素偏移均小于0.2個像素，最小達到0.13個像素，精度提升了一個數量級。同時，該方法具有較強的魯棒性，可獲取更為精確的水域輪廓，對利用星載SAR圖像檢測水庫、湖泊等水域變化有較好的工程使用價值。