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

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

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

0 引言

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

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

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

1 亞像素級水陸分割方法

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

圖1 本文方法流程圖

1.1 非局域相干斑濾波算法

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

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

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

A≥0

(1)

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

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

(2)

(3)

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

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

(4)

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

FCM是一種無監(jiān)督的聚類算法。一幅圖像可描述為集合I={i1,i2,…,iN}，傳統(tǒng)FCM算法通過設(shè)置c個不同的聚類中心，計算集合I中每個像素與各聚類中心的隸屬度，完成對圖像的分類。FCM算法執(zhí)行時，需遍歷圖像中每一個像素計算隸屬度，當處理大尺度SAR圖像時，非常耗時。針對該問題，文獻[4]提出了一種基于灰度分布的FCM算法，采用如式(5)所示的目標函數(shù)Jm，將與聚類中心隸屬度的計算由像素轉(zhuǎn)變?yōu)榛叶燃墸诨叶燃壣贤瓿蓪D像的分類，從而大大提高了算法的運行效率。

(5)

(6)

(7)

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

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

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

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

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

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

1.3 亞像素級水陸分割方法

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

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

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

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

(1)自然水域周邊地形一般較為復(fù)雜，陡峭山脈等地物會因起伏地形未受到電磁波照射，沒有回波，會產(chǎn)生較多的陰影現(xiàn)象，在SAR圖像中呈現(xiàn)為黑色區(qū)域。而水面對電磁波產(chǎn)生表面散射，后向散射系數(shù)較低，平靜水域在SAR圖像中也同樣呈現(xiàn)為與陰影相似的黑色區(qū)域。為實現(xiàn)對二者的區(qū)分，需進一步提取感興趣區(qū)域(Region of Interest, ROI)。在粗分割獲取的像素級水域中，記最大連通區(qū)域面積為Smax，設(shè)置面積閾值Ta=μ·Smax。當連通區(qū)域面積大于Ta時，被記為ROI，否則丟棄。同時，將提取到的ROI區(qū)域邊界信息記為像素級水域輪廓。本文根據(jù)經(jīng)驗將μ設(shè)置為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)。方程階數(shù)取決于SPW的大小，SPW越大，包含的像素級水域輪廓點越多，多項式階數(shù)越高。通常，多項式階數(shù)過高可能會導(dǎo)致最終結(jié)果不具備魯棒性，相反小的SPW則會對高階多項式要求較小。

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

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

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

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

(8)

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

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

(9)

式中，ψ為SPW內(nèi)的圖像，ψ0為SPW中的水域區(qū)域，σψ0為由SPW中擬合的fwl函數(shù)獲得的曲線，ρ為常數(shù)參數(shù)，本文設(shè)置ρ=1。

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

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

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

2 實驗結(jié)果與分析

2.1 研究區(qū)域及數(shù)據(jù)

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

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

圖4 丹江口水庫位置

2.2 實驗結(jié)果

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

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

圖5 相干斑濾波結(jié)果

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

圖6 亞像素級水陸分割線提取結(jié)果

2.3 對比實驗與定量分析

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

圖7 不同方法水域提取結(jié)果

(1)虛警率

(10)

(2)準確率

(11)

(3)平均輪廓偏移距離

(12)

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

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

表2 定量分析結(jié)果

3 結(jié)束語

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