李玲玉,張 毅

(中國科學(xué)院電子學(xué)研究所,北京 100190)

0 引 言

合成孔徑雷達(dá)(Synthetic Aperture Radar, SAR)是一種高分辨率主動式微波成像系統(tǒng)，不受復(fù)雜氣候等條件的限制，可全天時全天候地獲取影像，甚至可透過地表、植被、云層獲取其掩蓋的信息。且SAR具有成本低、覆蓋廣的特點，使其在水體變化檢測等領(lǐng)域具有獨特的優(yōu)勢[1]。利用SAR進(jìn)行水體變化的檢測可優(yōu)化水體資源配置，完善防洪減災(zāi)體系，完善湖泊等生態(tài)環(huán)境建設(shè)，進(jìn)行水體附近綜合管理。各種應(yīng)用的難點在于如何提高SAR圖像水體變化檢測精度[2]。

為精確得到SAR圖像的變化信息，研究人員對其進(jìn)行了大量的研究。變化檢測方法有直接比較法和分類后比較法。研究人員通常用直接比較法，首先利用變化信息構(gòu)建差異圖，然后針對差異圖進(jìn)行分析，區(qū)分出變化類與非變化類。

傳統(tǒng)的差異圖生成法主要包含如下3類：1)基于像素的算子，如Dekker[3]提出的比值法(LR)、Celik[4]提出的對數(shù)比值法(LR-D)等；2)基于鄰域的算子，如Inglada等[5]提出的均值比算子(MR)、Gong等[6]提出的鄰域比值法(NR)，以及Xiong等[7]提出的最大似然比算子(LLR)等；3)幾種方法的結(jié)合，因為一般單一的差異圖生成算子不能很好地反映出圖像之間的差異。如Hou等[8]提出的LR與高斯對數(shù)比結(jié)合的方法，Zheng等[9]提出的差值法與比值法結(jié)合的方法，Yan等[10]提出的LR與對數(shù)均值比結(jié)合的方法。上述方法均沒有結(jié)合單像素算子和鄰域算子的優(yōu)點。

差異圖分析法包括閾值法、主動輪廓法和聚類法等。閾值法能快速區(qū)分變化區(qū)域和非變化區(qū)域，然而當(dāng)變化類和不變類之間存在較強(qiáng)重疊或無法準(zhǔn)確建模時，使用閾值法提取結(jié)果較差。主動輪廓法可以在一定程度上降低相干斑噪聲的影響，但在初始輪廓線與真實輪廓線相距較遠(yuǎn)時誤差較大。聚類法，如Bezdek等[11]提出的模糊C均值(FCM)方法，它較依賴于人工給出的參數(shù)，沒有考慮鄰域信息；Krinidis等[12]提出的局部信息模糊C均值聚類(FLICM)方法，解決了上述問題，但仍容易被孤立噪聲點影響；Gong等[13]提出的基于改進(jìn)MRF能量函數(shù)的FCM(MRFFCM)算法，降低了散斑的影響，但不能保證噪聲與細(xì)節(jié)的平衡。而空間、紋理等信息可降低孤立噪聲的影響，引入更多的SAR圖像特征。

傳統(tǒng)方法精度有限，為了提高水體變化檢測精度，本文結(jié)合單像素算子和鄰域算子的優(yōu)點，將空間、紋理等信息與適應(yīng)調(diào)節(jié)抑制噪聲和細(xì)節(jié)平衡的模糊因子結(jié)合運用到算法中，提出一種基于新差異算子和紋理的SAR圖像變化檢測方法，該方法分為3步:

1)根據(jù)SAR圖像的特征，結(jié)合LR和LLR這2種算子，提出一種新的差異算子。

2)提出一種新的FLICM方法，將目標(biāo)分為3類。

3)根據(jù)差異圖像分割閾值將過渡區(qū)域再次分類，得到最終的變化檢測結(jié)果。

本文算法的具體流程如圖1所示。

圖1 本文的SAR圖像水體變化檢測算法流程

1 本文的SAR圖像水體變化檢測算法

1.1 差異圖像的生成

針對SAR圖像的噪聲特征，比值法的比率運算符被廣泛使用于SAR圖像的變化檢測中，它通常以對數(shù)刻度表示，可將SAR圖像的乘性噪聲變?yōu)榧有栽肼?，增?qiáng)了SAR圖像對輻射誤差的魯棒性。本文基于比值法的LR算子與LLR算子，分析其特性，提出將兩者結(jié)合的變化檢測新算子。

LR是應(yīng)用非常廣泛的比率檢測算子，可以很好地反映單像素差，但其忽視了鄰域信息，LR的計算公式為：

(1)

其中，X1、X2為變化前后的圖像。DLR最小值為1，代表不變的像素。而變化像素的值通常較大，因此采用對數(shù)運算符將其縮小。采用對數(shù)運算符后，最小值對應(yīng)為0，采用絕對值運算符把正負(fù)像素視為變化像素。因此LR變?yōu)椋?/p>

(2)

水體一般為鏡面散射，因此水體在SAR圖像中近似為均勻區(qū)域，伽馬分布能很好地擬合水體的SAR圖像。LLR算子假設(shè)SAR圖像服從伽瑪分布，并構(gòu)造了似然比檢驗[7]，采用3×3的矩形窗代表局部信息，既能反映鄰域差異，又具有較好的抗噪聲能力。仍然將對數(shù)運算符及絕對值運算符運用到簡化后的LLR算子中計算：

(3)

其中，Ω為圖像像素(m,n)在3×3窗口內(nèi)的像素集。上述的2個算子，不變像素的值在0附近形成了一個尖峰，變化的像素形成一條長尾。最后將差異圖像歸一化到[0,255]。本文利用1997年7月和8月加拿大渥太華的2幅RADARSAT SAR圖像作為數(shù)據(jù)集1來展示本文算法的優(yōu)勢，該圖像分辨率為12 m，大小為350×290，如圖2(a)和圖2(b)所示，圖2(a)為洪水后的SAR圖像，圖2(b)為洪水前的SAR圖像，參考差異圖像如圖2(c)所示，DLR、DLLR歸一化直方圖如圖3(a)和圖3(b)所示。

(a) 7月原始圖像

(b) 8月原始圖像

(c) 參考差異圖像

(a) DLR歸一化直方圖

(b) DLLR歸一化直方圖

本文通過將LR算子與LLR算子相乘構(gòu)成SAR圖像水體變化檢測新算子。在單像素差和領(lǐng)域差同時都很強(qiáng)時，新算子的值將會變大，在單像素差和領(lǐng)域差同時都很弱時，新算子的值將會變小。而在有噪聲造成單像素差弱但鄰域差強(qiáng)，或單像素差強(qiáng)但鄰域差弱時，新算子的值介于兩者之間。最終的差異圖生成新算子為：

DI=DLR×DLLR

(4)

上述DLR、DLLR、DI這3種算子生成的差異圖像如圖4(a)～圖4(c)所示。

(a) DLR差異圖

(b) DLLR差異圖

(c) DI差異圖

DI歸一化直方圖如圖5(a)所示。本節(jié)利用差異圖像的相鄰直方圖比值進(jìn)行閾值分割，得到初始分割圖像，如圖5(b)所示。

(a) DI歸一化直方圖

(b) 差異圖像DI的相鄰直方圖比值圖

差異矩陣直方圖中，峰值表示無變化區(qū)域，震蕩部分為變化區(qū)域，可以得到變化與非變化的檢測門限在其中間，具體值可用相鄰直方圖比值圖得出。差異圖像的相鄰直方圖比值計算公式為：

(5)

其中，hDI為差異圖像直方圖，i∈[i0,255)，i0為hDI最大值對應(yīng)的灰度。自動選擇比率小于1的點作為過渡點的第一個點為閾值,即：

k={i|η(i)1}

(6)

分析DLR、DLLR、DI算子的直方圖，所提出的DI的直方圖在不變部分顯示出更尖銳的峰值，在變化部分顯示出更長的尾部。該算子放大了變化區(qū)域與不變區(qū)域的特征，同時降低了噪聲效果。

分析差異圖像的相鄰直方圖比值圖，如圖6(a)所示，得到差異圖的初始分割。DLLR與DI的初始分割如圖6(b)、圖6(c)所示。其中，DLLR的閾值為文獻(xiàn)[14]中所提出的閾值。與參考圖對比，可以得到DLLR的PCC為0.9611，kappa系數(shù)為0.8392；DI的PCC為0.9746，kappa系數(shù)為0.9038。對差異圖進(jìn)行閾值分割，DI的變化檢測結(jié)果好于DLLR。

(a) 相鄰直方圖比值圖

(b) DLLR差異圖初始分割

(c) DI差異圖初始分割

1.2 基于紋理的差異圖獲取

現(xiàn)今，有許多紋理分析技術(shù)，如灰度共生矩陣(GLCM)、傅里葉功率譜、馬爾科夫隨機(jī)場(MRF)和Gabor濾波器[14～17]等。其中，Gabor濾波器在分類等問題上應(yīng)用廣泛，文獻(xiàn)[18]驗證了Gabor與其他紋理特征相比，可以在計算成本低的情況下達(dá)到同樣的性能。二維高斯濾波器公式如下：

(7)

利用尺度和方向的不同，將圖像與二維Gabor核進(jìn)行卷積，得到了DI(x，y)的N組Gabor響應(yīng)。本文選擇尺度和方向個數(shù)分別為4、6。Gabor濾波器及其對應(yīng)的響應(yīng)如圖7(a)和圖7(b)所示。

(a) Gabor濾波器

(b) 對應(yīng)的響應(yīng)

(c) 生成的紋理

將每組響應(yīng)變?yōu)橐涣校玫紾abor差分圖像矩陣為：

Gabor=(gbi1,gbi2,…,gbii,…,gbin)

(8)

其中，n=24，gbi為每組Gabor濾波器對應(yīng)的響應(yīng)。除紋理信息外，本節(jié)還加入了2個空間特征，它們?yōu)闅w一化后的圖像x坐標(biāo)和y坐標(biāo)。

因此，差分圖像的像素可以由這26維特征向量表示，為了減少冗余，采用PCA進(jìn)行降維，留下特征值超過最大特征值的5%的特征向量。最后，采用公式(9)將其變?yōu)橐痪S進(jìn)行分析[19]。

(9)

將其整形為原始圖像大小，得到的結(jié)果如圖7(c)所示。

1.3 改進(jìn)的FLICM方法

改進(jìn)的FLICM方法稱為結(jié)合紋理的FLICM方法(FLICM_texture)，該方法分為3個步驟：估計初始聚類中心、FLICM_texture的聚類和閾值分割的再次分類。

1.3.1 初始聚類中心的獲得

由于變化類與非變化類之間有一定的過渡區(qū)域，該區(qū)域有灰度級低的植被，或者受到了噪聲影響。因此本文將差異圖分為3類：變化區(qū)域、非變化區(qū)域和過渡區(qū)域。

由1.1節(jié)中提出的DI差異圖像的相鄰直方圖比值閾值法提取出核心變化區(qū)域，即：

(10)

其中，k為1.1節(jié)中對應(yīng)的閾值。經(jīng)過對比分析，該閾值分割的精度較高，因此可利用它得到3個初始聚類中心：

(11)

1.3.2 FLICM_texture方法

2010年，Stelios在FCM目標(biāo)函數(shù)中引入一個新的因子，提高了對噪聲和異常值的魯棒性，不需提前計算參數(shù)，不再需要人工選擇參數(shù)。該模糊因子計算公式如下：

(12)

其中，k為聚類中心，j為第i個像素局部窗口的像素集合，m為每個模糊隸屬度的加權(quán)指數(shù)。dij為像素i與j之間的空間歐氏距離，ukj為第j個像素在第k個聚類中心的隸屬度，vk為聚類k的中心?？梢钥闯觯蜃覩ki沒有使用任何控制圖像噪聲和圖像細(xì)節(jié)之間平衡的參數(shù)，且該因子包含了局部信息。

將局部空間和灰度信息應(yīng)用于目標(biāo)函數(shù)中，并將紋理信息添加到目標(biāo)函數(shù)中，得到最后的目標(biāo)函數(shù)如下式：

(13)

其中，yi、wk分別為紋理圖像中像素與聚類中心的值。使上式最小化的隸屬度及聚類中心表達(dá)式為：

(14)

(15)

FLICM_texture算法流程如下：

1)設(shè)置類別數(shù)c、模糊參數(shù)m和停止條件ε。

2)將1.3.1節(jié)得到的初始聚類中心設(shè)為該算法的初始聚類中心，計算紋理聚類中心，設(shè)置循環(huán)數(shù)b=0。

3)用式(14)計算隸屬度。

4)用式(11)修正聚類中心，計算新的紋理聚類中心。

5)當(dāng)max{u(b)-u(b+1)}<ε時停止，否則，令b=b+1，進(jìn)入步驟3。

當(dāng)算法收斂后，進(jìn)行去模糊化處理，將滿足下面條件的像素i劃分給第k類。

Ci=argk{max {uki}},k=1,2,…，c

(16)

由式(16)可以看出，該方法并未增加FLICM的時間復(fù)雜度。

生成紋理時，二維Gabor濾波器中的行數(shù)m的選取過程：利用FLICM_texture算法檢測圖像，當(dāng)m=9時，PCC為0.9539，kappa系數(shù)為0.8047，F(xiàn)A為0.0083，MA為0.2855。當(dāng)m=39時，PCC為0.9579，kappa系數(shù)為0.8238，F(xiàn)A為0.0083，MA為0.2599。m=39時，其PCC、kappa系數(shù)比m=9時高，F(xiàn)A、MA比m=9時低，顯然，m=39時的檢測效果好于m=9，因此令二維Gabor濾波器中的行數(shù)為39，其中PCC、kappa、FA、MA系數(shù)說明見第2章，生成的紋理圖像如圖8所示。

(a) m=9時紋理圖像

(b) m=39時紋理圖像圖8 紋理圖像

在一切參數(shù)相同的條件下，對FLICM和FLICM_texture將差異圖分為2類，得到的檢測結(jié)果如圖9(a)和圖9(b)所示。在本文的方法中，去除FLICM_texture的紋理及完整的本文方法得到的結(jié)果如圖9(c)和圖9(d)所示。

(a) FLICM結(jié)果 (b) FLICM_texture結(jié)果

(c) 不加紋理的本文方法 (d) 本文方法

對于圖9(a)和圖9(b)，與參考圖像對比，可以得到FLICM的PCC為0.9499，kappa系數(shù)為0.7843，MA為0.3152；FLICM_texture的PCC為0.9580，kappa系數(shù)為0.8238，MA為0.2599。PCC增長了約0.85%，kappa系數(shù)增長了約5.0%，MA下降了約17.5%。對于圖9(c)和圖9(d)，圖9(c)的PCC為0.9761，kappa系數(shù)為0.9067，MA為0.1226；圖9(d)的PCC為0.9800，kappa系數(shù)為0.9203，MA為0.0776。這驗證了FLICM_texture方法的有效性。

若只是使用FLICM_texture方法將差異圖等分為2類則精度較低，因此用FLICM_texture方法將其分為3類：水體、背景和過渡區(qū)域，使用閾值法將過渡區(qū)域再次劃分為2類。

1.3.3 閾值分類

由于變化類與非變化類之間有一定的過渡區(qū)域，該區(qū)域有灰度級低的植被，或者受到了噪聲影響。本文采用閾值法對其進(jìn)行二次分類。判斷公式為：

(17)

其中，Ωg為差分圖像過渡區(qū)域的像素集合，Dm為該區(qū)域的像素，k為1.1節(jié)中利用相鄰直方圖比值得到的閾值。

2 定量分析變化精度

將本文提出的方法得到的變化檢測結(jié)果與FCM、模糊局部C均值法(FLICM)、基于馬爾科夫隨機(jī)場的ICM算法(MRF&ICM)[20]、閾值法做對比，利用以下4個定量指標(biāo)分析算法的有效性：1)虛警(FA)表示錯誤檢測到的未變化像素的百分比;2)漏警(MA)表示未檢測到的已變化像素的百分比;3)kappa系數(shù)表示2幅圖像的一致性程度;4)正確分類率(PCC)表示百分比正確分類[10]。

表1 各種方法在加拿大渥太華上空數(shù)據(jù)集的變化檢測結(jié)果

1)數(shù)據(jù)集1為渥太華數(shù)據(jù)集，它分別是RADARSAT在1997年5月和8月在加拿大渥太華上空拍攝的SAR圖像的切片(290×350)，它們?yōu)樵?jīng)遭受洪水侵襲的地區(qū)。各種方法在渥太華數(shù)據(jù)集的變化檢測精度如表1所示，各種方法在渥太華數(shù)據(jù)集的變化檢測結(jié)果圖如圖10所示。

(a) FCM+DLR (b) FLICM+DLR

(c) MRF&ICM+DLR (d) FCM+DLLR

(e) FLICM+DLLR (f) 閾值+DLLR

(g) MRF&ICM+DLLR (h) DI+閾值

(i) 本文方法 (j) 參考圖

2)數(shù)據(jù)集2為伯爾尼數(shù)據(jù)集，它分別是ERS-2在1999年4月和5月在瑞士伯爾尼市上空拍攝的SAR圖像的切片(301×301)，如圖11(a)和圖11(b)所示，DI差異圖如圖11(c)所示。在這2個日期之間，部分城市、伯爾尼機(jī)場遭到了阿雷河的淹沒。各種方法在伯爾尼數(shù)據(jù)集的變化檢測精度如表2所示，各種方法在伯爾尼數(shù)據(jù)集的變化檢測結(jié)果如圖12所示。

(a) 4月原始圖像

(b) 5月原始圖像

(c) DI生成的差異圖圖11 伯爾尼原始圖像及生成的差異圖

表2 各種方法在伯爾尼數(shù)據(jù)集的變化檢測結(jié)果

(a) FCM+DLR (b) FLICM+DLR

(c) MRF&ICM+DLR (d) FCM+DLLR

(e) FLICM+DLLR (f) 閾值+DLLR

(g) MRF&ICM+DLLR (h) DI+閾值

(i) 本文方法 (j) 參考圖

3)數(shù)據(jù)集3為印度金奈數(shù)據(jù)集，它分別是哨兵在2018年8月4日和2019年7月30日于金奈上空拍攝的切片(1373×1927)，如圖13(a)和圖13(b)所示，DI差異圖如圖13(c)所示。其中金奈2019年遭受嚴(yán)重干旱，造成湖泊等水體大面積縮減[21]。各種方法在金奈數(shù)據(jù)集的變化檢測精度如表3所示，各種方法在金奈數(shù)據(jù)集的變化檢測結(jié)果如圖14所示。比較水體變化面積可以得到，金奈附近水體相比于去年減少了73.24%，其中金奈最大的湖Puzhal湖面積減少了近83.88%。

(a) 2018年8月原始圖像

(b) 2019年7月原始圖像

(c) DI差異圖

表3 各種方法在金奈數(shù)據(jù)集的變化檢測結(jié)果

(a) FCM+DLR (b) FLICM+DLR

(c) MRF&ICM+DLR (d) FCM+DLLR

(e) FLICM+DLLR (f) 閾值+DLLR

(g) MRF&ICM+DLLR (h) DI+閾值

(i) 本文方法 (j) 參考圖

綜上可以看出，由于差異圖中引入局部平均，差異圖像中一些細(xì)小的變化丟失了，但變化丟失相對較少。比較由DLLR和DI得到閾值分割后的結(jié)果的PCC與kappa系數(shù)，發(fā)現(xiàn)本文提出的新算子具有良好的性能。將本文方法與其他方法相比，本文方法提高了檢測率，降低了虛警率，提高了kappa系數(shù)，因此本文方法具有較高的魯棒性。

3 結(jié)束語

對SAR圖像水體變化進(jìn)行檢測具有完善防洪減災(zāi)體系、進(jìn)行水體附近綜合管理等作用。本文針對多時相SAR變化檢測中孤立噪聲點、需人工輸入部分參數(shù)、信息應(yīng)用不全等問題，提出了基于新差異算子和紋理的變化檢測算法。

將本文方法應(yīng)用于加拿大渥太華上空、瑞士伯爾尼市上空、印度金奈上空的SAR水體變化檢測，可以看出，本文方法在3個地區(qū)的PCC和kappa系數(shù)上均優(yōu)于FCM、FLICM、MRF&ICM、閾值法，因此本文方法提高了SAR圖像水體變化檢測精度。但是，本文提出的新差異算子對低后向散射區(qū)域較為敏感，如果SAR圖像中有大的低后向散射區(qū)域，虛警率可能會變大。下一步將針對低后向散射區(qū)域進(jìn)行處理，使其更具有通用性。