王 平 王宜懷 劉長勇 彭 濤

1(武夷學院數(shù)學與計算機學院 福建 武夷山 354300)2(蘇州大學計算機科學與技術學院 江蘇 蘇州 215006)3(認知計算與智能信息處理福建省高校重點實驗室 福建 武夷山 354300)

0 引 言

合成孔徑雷達(SAR)圖像的變化檢測，目標是生成一個變化圖，表述兩時相或多時相圖在不同標定時間之間的變化情況，其廣泛應用在環(huán)境監(jiān)測、生態(tài)監(jiān)測、城市發(fā)展研究、農林監(jiān)測、自然災害評估等領域[1]。SAR圖像以微波成像，圖像強度值與地面目標的后向散射成正比，易產生不易濾波處理的斑點噪聲，進而對變化檢測精度造成較大影響。如何在噪聲抑制和變化細節(jié)保留之間進行均衡，提升變化檢測精度，是處理過程中需要解決的重要問題。

從SAR圖像變化檢測技術流程來看，主要有圖像預處理、變化信息發(fā)現(xiàn)、區(qū)域提取或變化類型確定等幾個步驟。其中，信息發(fā)現(xiàn)和區(qū)域提取是變化檢測的研究熱點，且方法眾多。從差異圖(DI)角度進行變化檢測處理是一種常用的途徑，其辦法是先構建兩時相SAR圖的DI，然后再對DI分析處理。

DI構建方法中，傳統(tǒng)的方法有差值法、比值法(R)[2]、改進比值法(IR)[3]、對數(shù)比值法(LR)[4]、均值法等，但這些方法均沒有利用像素點周圍的鄰域空間變化信息來減少斑點噪聲的影響，進而造成精度較低。近年來提出的組合差異法的思路是將各種方法的優(yōu)點結合起來，特別是尋求噪聲抑制和圖像變化細節(jié)保留之間的均衡[5]。MR方法在處理中引入像素點的鄰域空間均值信息，能抑制斑點噪聲，圖像平滑，但未考慮區(qū)域異質性，使得變化邊緣等異質性較大的地方檢測不穩(wěn)定[6]。Liu等[7]利用MR方法和差值法進行組合獲取DI，具有較好的效果。而NR方法采用自適應斑點噪聲濾波算法常用的局部鄰域異質性測度來控制鄰域空間信息對中心像素點變化檢測的影響，其本質是將IR方法和空間鄰域均值信息通過異質性測度來融合，但DI像素樣本分布不利于DI分析[8]。

在DI分析中，常用方法是將DI的像素點分為變化類和不變化類。文獻[5，7]針對DI采用k為2的K-means聚類，但對斑點噪聲影響的魯棒性較差，誤檢率較大。王建明等[9]采用基于自適應距離的FCM算法，將DI聚類出變化和不變化的區(qū)域類別，精度有所提高。文獻[4]采用基于遺傳算法改進的FCM算法對DI進行處理，識別變化類和不變化類，進而得到變化檢測圖。文獻[10]在DI上基于多粒度級聯(lián)森林和多尺度融化的方法進行處理，進而得到變化檢測圖。文獻[11]設計了識別DI特征空間的層次級聯(lián)FCM算法來得到變化和未變化的信息。但上述方法的樣本分類依據(jù)僅僅來源于DI，未考慮DI中殘留噪聲的影響，僅僅是粗粒度的分類。文獻[12]在利用FCM將DI像素樣本聚類成三類的基礎上，利用加速遺傳算法進一步對不確定變化與否的像素樣本進行分類進而獲得變化檢測圖。文獻[13]設計了層次FCM算法來分析DI得到樣本標簽，然后再將樣本結合原始圖像鄰域塊進行SVD網(wǎng)絡分類，進而得到變化檢測圖，對DI的噪聲殘留有一定的抑制效果。崔斌等[14]在層次FCM聚類的基礎上，利用卷積神經(jīng)網(wǎng)絡對中間類進一步進行分類處理，但該方法對中間類處理算法復雜，耗時較長。

基于此，本文依托IR方法的高效性和NR方法抗噪性強的特性，提出了新的鄰域比方法INR，使得所構建的DI樣本分布更利于分析。而在DI分析中，常用的直接聚類得到的檢測結果忽略了DI中殘留噪聲的影響。本文利用層次聚類獲得樣本標簽、輸入時相圖鄰域塊特征向量作為樣本進行ELM網(wǎng)絡分類訓練而生成變化檢測二值圖，依托ELM網(wǎng)絡簡單有效的非線性映射能力，能解決聚類分析方法中有用信息丟失的問題，尤其是殘留噪聲影響的數(shù)據(jù)判別，實現(xiàn)對中間類細粒度處理，使得檢測效果明顯。

1 方法框架

假設有提前配準預處理的兩時相SAR圖像I1和I2，則本文所提方法的總體框架如圖1所示。

圖1 本文所提方法的總體框架

本文所提方法的具體步驟如下：

(1) 基于INR的DI構建。

(2) 對DI進行層次FCM聚類，并得到聚類差異圖，該矩陣元素值為0、0.5和1。其中，0和1分別代表DI對應像素點屬于不變化類和變化類，0.5代表可能噪聲殘留影響的中間類。

(3) 根據(jù)FCM聚類結果的變化類和不變化類像素，且以這些像素點為中心，提取兩時相圖的k×k圖像塊并組合在一起作為鄰域塊特征向量，依托ELM網(wǎng)絡簡單有效的非線性映射能力，構造ELM網(wǎng)絡的訓練樣本集，訓練ELM網(wǎng)絡，并將兩時相圖的所有k×k鄰域塊特征向量作為測試樣本集。

(4) 根據(jù)訓練好的ELM網(wǎng)絡，對測試樣本集進行分類，可實現(xiàn)將FCM聚類結果值為0.5的中間類像素進一步確定為變化類和不變化類，最終得到與該兩時相SAR圖關聯(lián)的變化檢測二值圖。

2 相關要點原理

2.1 差異圖矩陣(DI)計算

SAR圖像的配準預處理雖然實現(xiàn)了圖像斑點噪聲濾波和對齊等處理，但其殘留的乘性斑點噪聲仍然會嚴重影響變化檢測的精度，因而，好的DI生成方法能提高變化檢測的精度。假設In(i,j)表示SAR圖像In中處于(i,j)位置上的像素強度值，i∈[1,M]，j∈[1,N]，其中，M、N分別對應時相SAR圖像In的維度，則R方法可以表示為：

(1)

式中：DR表示采用R方法得到的DI；DR(i,j)表示DR中處于(i,j)位置的DI的像素強度值。經(jīng)過分析，R方法只能檢測諸如像素點強度值增大或減小變化類型中的一種，基于R方法改進得到IR方法則可以實現(xiàn)各種變化類型的檢測。IR方法可以表示為：

(2)

則DIR(i,j)值越大，代表該像素點對應的原始輸入兩時相SAR圖的像素未發(fā)生變化，反之，則發(fā)生了變化。實際應用表明，像素比值法對校正和輻射測量的誤差都具有較好的魯棒性，且融合SAR圖像的鄰域強度均值信息生成的DI可以減少乘性或加性斑點噪聲的影響，這也是提出MR、NR方法的依據(jù)。MR方法直接根據(jù)IR方法來改進，其表示為：

(3)

式中：u1(x)、u2(x)分別表示SAR圖像I1和I2的以(i,j)位置為中心鄰域空間x(大小為r×r)的強度均值。DMR(i,j)值越大，代表該像素點對應的原始輸入兩時相SAR圖的像素未發(fā)生變化，反之，則發(fā)生了變化。分析式(3)可知，MR方法不僅可以檢測各種類型的變化，且在一定程度上能抑制斑點噪聲的影響，能提高兩時相SAR圖像變化檢測精度。NR方法正是結合了IR和MR方法優(yōu)點，采用異質性測度值θ來作為權重比，其表示為：

(4)

式中：Ωx表示兩幅SAR圖像以(i,j)位置為中心的鄰域空間x的集合；I1(x)、I2(x)定義為取以(i,j)位置為中心的鄰域空間x塊矩陣像素強度值；θ=σ(x)/μ(x)是鄰域Ωx的異質性測度，σ(x)是鄰域Ωx的標準差，μ(x)是鄰域Ωx的均值。較小的θ對應一個同質區(qū)域，則式(4)的后一部分所占比重大，反之，較大的θ對應一個異質區(qū)域，則式(4)的前一部分所占比重大。

基于NR思想，本文在IR的基礎上引入異質性測度θ，實現(xiàn)融合中心像素點強度和局部鄰域均值信息的線性組合，形成本文提出的INR方法，其表示為：

DINR(i,j)=

(5)

式中：DINR表示采用INR方法得到的DI；DINR(i,j)表示DINR中位置處于(i,j)位置的DI的像素強度值；u1(x)、u2(x)定義如前所述；θ1和θ2分別是SAR圖像I1和I2是鄰域x的異質性測度；常數(shù)C解決分母可能趨近于0問題，常根據(jù)C=(kL)2取值，L為SAR圖像灰度級別，如8位灰度強度則取值為255，k一般取值為0.003。

分析式(5)可知，DINR(i,j)值越大，則兩時相SAR圖間的相似度越大，對應的像素點也越可能沒有變化，反之，DINR(i,j)值越小則越可能發(fā)生了變化。這可解釋為：若有一個以像素位置(i,j)為中心的同質區(qū)域，此時θ取值小，表明此區(qū)域在前后時相圖像中是未變化的，則式中鄰域x的均值信息u(x)占比重較大，而可能受斑點噪聲干擾的中心像素點強度值I(i,j)占比重小，使得比值法的min/max模式能很好地識別出該像素點及其鄰域空間的不變性；反之，如果是異質區(qū)域，則θ取值大，在式(5)中，中心像素點強度值I(i,j)占比重大，此時min/max模式又更能細致反映各種變化類型的像素點變化情況，實現(xiàn)了變化細節(jié)保留。

另外一個問題是異質性測度值歸一化，式(5)中θ1和θ2的理想值為[0,1]，但是當σ(x)>μ(x)時，可能出現(xiàn)該值大于1情況，而造成DINR(i,j)取值可能為負，無法理想地表達變化關系，可以使用下式對異質性測度值進行統(tǒng)一歸一化：

(6)

2.2 層次FCM聚類

根據(jù)所提的INR方法進行計算得到的DI，通過FCM聚類方法[15]，將其元素分3類，即變化類Ωc、中間類Ωi和不變化類Ωu。其中，變化類Ωc和不變化類Ωu為變化概率和不變化概率高的像素點，而中間類Ωi呈現(xiàn)變化與否的不確定性，其可能是噪聲殘留的影響。進一步，由于單一的FCM聚類可能會造成中間類Ωi的比重過大，使得有用信息的丟失過大，影響最終變化檢測的精度，因而采用層次FCM聚類的方法，具體方法的步驟如下：

輸入：差異圖矩陣DI

(7)

(8)

(5) 轉至第(3)步，直到t=m后退出。

輸出：層次FCM聚類結果得到表示為帶有標簽{Ωc,Ωi,Ωu}的DI像素集合。

經(jīng)過層次FCM聚類后，我們得到了一張標有{Ωc,Ωi,Ωu}的變化圖。屬于Ωc的像素變化概率較大且令其值為0標簽，屬于Ωu的像素不變化概率較大且令其值為1標簽，這兩種像素可以構造ELM網(wǎng)絡的訓練樣本。而屬于Ωi的標簽值為0.5的像素為不確定變化與否的中間類，其進一步的標簽確定將由ELM網(wǎng)絡根據(jù)空間鄰域塊特征向量分類而完成，這樣就解決了可能的DI殘留噪聲影響。

2.3 ELM二分類

根據(jù)2.2節(jié)層次FCM聚類后的結果，且為了進一步解決Ωi的類別問題，將Ωc和Ωu作為兩輸入時相圖的鄰域塊特征向量的樣本標簽構造分類器的訓練樣本集，訓練ELM網(wǎng)絡[16]，并將所有鄰域塊特征樣本構造成ELM網(wǎng)絡的測試樣本集，最終可得到變化檢測圖。具體方法的步驟如下：

(1) 生成以Ωc和Ωu類對應像素為中心的兩輸入時相圖的鄰域塊特征，并構造標簽為Ωc和Ωu的訓練樣本集。其中，鄰域塊特征向量構造方法是以對應位置像素為中心，取大小為k×k的像素塊構成，將兩輸入時相圖的鄰域塊特征按如圖2所示(圖中k取3)形式轉換并連接為行向量形式，作為ELM網(wǎng)絡的樣本特征向量，再加上中心像素點的類別標簽，即構成ELM網(wǎng)絡的訓練樣本集。同理，按同樣的方法按序選取所有的像素為中心，取大小為k×k的像素塊構造圖像鄰域塊特征，即構成ELM網(wǎng)絡的測試樣本集。

圖2 兩輸入時相圖鄰域特征向量的生成方法

(2) 基于訓練樣本集對ELM網(wǎng)絡進行訓練。

(3) 基于測試樣本集，利用訓練好的ELM進行二分類，最終得到變化檢測圖。

3 實 驗

3.1 數(shù)據(jù)集描述和實驗設置

實驗數(shù)據(jù)集使用三組真實且已配準預處理的SAR圖像數(shù)據(jù)，分別是Bern地區(qū)數(shù)據(jù)、Ottawa地區(qū)數(shù)據(jù)和Sardinia地區(qū)的兩時相圖及手動真值圖，如圖3-圖5所示。

圖3 Bern地區(qū)數(shù)據(jù)

圖4 Ottawa地區(qū)數(shù)據(jù)

圖5 Sardinia地區(qū)數(shù)據(jù)

為了驗證本文算法的有效性，具體實驗設置如下：

(1) 四種DI構建方法的對比：本文所提的INR方法與IR方法、MR方法和NR方法，通過其所生成的DI及直方圖進行分析。

(2) 變化檢測二值圖對比：將四種方法生成的DI按本文設計的DI分析方法取得變化檢測二值圖，并進行性能的對比。

(3) 本文方法與其他典型的幾種變化檢測方法對比，以評估本文方法的有效性。

本文所提算法的參數(shù)設置為鄰域空間x均值矩陣 大小r取值為3，層次FCM聚類的第一輪聚類數(shù)為3，第二輪聚類數(shù)m為7，聚類預設參數(shù)σ1取值為1.12，σ2取值為1.10，鄰域塊特征矩陣大小k取值為5，ELM分類器輸入層節(jié)點數(shù)為50，隱含層節(jié)點數(shù)為10，輸出層節(jié)點數(shù)為2。數(shù)據(jù)對比采用較為廣泛使用的定量分析方法，假設圖像總像素個數(shù)為N，則定量數(shù)據(jù)包括漏檢像素數(shù)FN、誤檢像素數(shù)FP、總錯誤數(shù)OE、正確率PCC和KAPPA系數(shù)。

OE=FN+FP

(9)

PCC=(N-OE)/N

(10)

KAPPA=(Pr(a)-Pr(e))/(1-Pr(e))

(11)

式中：TP為確定變化的像素且在變化檢測中被正確檢測為變化的像素數(shù)；TN為確定未變化的像素且在變化檢測中被正確檢測為未變化的像素總數(shù)；Pr(e)=((TP+FN)×(TP+FP)+(TN+FN)×(TN+FP))/N2，Pr(a)=(TN+TP)/N。KAPPA系數(shù)值越趨近于1，表明變化檢測的效果越好，根據(jù)其一致性等級劃分，一般為當其取值為0.81～1時，得到的變化檢測結果跟手動真值圖更接近且趨于一致。

3.2 實驗結果及分析

(1) 四種DI構建方法的對比：本文INR方法和IR方法、MR方法和NR方法對兩組數(shù)據(jù)處理得到的DI及其直方圖分別如圖6-圖8所示。

(a) IR方法及其直方圖

(b) MR方法及其直方圖

(c) NR方法及其直方圖

(d) 本文INR方法及其直方圖圖6 Bern地區(qū)數(shù)據(jù)DI對比

(a) IR方法及其直方圖

(b) MR方法及其直方圖

(c) NR方法及其直方圖

(d) 本文INR方法及其直方圖圖7 Ottawa地區(qū)數(shù)據(jù)DI對比

(a) IR方法及其直方圖

(b) MR方法及其直方圖

(c) NR方法及其直方圖

(d) 本文INR方法及其直方圖圖8 Sardinia地區(qū)數(shù)據(jù)DI對比

差異圖DI反映兩時相SAR圖像之間的變化情況，尤其需突出變化區(qū)域和輪廓，且還應有高的噪聲抑制魯棒性。在其直方圖中，將會呈現(xiàn)兩個波峰的數(shù)據(jù)分布，其中未發(fā)生變化像素數(shù)的波峰應高、集中且更靠近零灰度值附近，而發(fā)生變化像素數(shù)對應的波峰應該更靠近1灰度值附近。從得到的DI及其直方圖可以看出，如圖6(a)、圖7(a)和圖8(a)所示，由于IR方法按單像素比的min/max求得DI，未考慮鄰域空間信息和區(qū)域異質性，對圖像殘留斑點噪聲的抑制不夠，細微的強度值變化也會在DI中反映出來，所得DI的直方圖毛刺較多，且直方圖分布也表明樣本數(shù)據(jù)峰值分布不合理，尤其是噪聲殘留較多的圖6(a)，沒有明顯地體現(xiàn)變化像素和未變化像素的區(qū)分，勢必影響后續(xù)的DI分析。對于MR方法，其采用鄰域均值比的min/max求得DI，能較好地抑制斑點噪聲的影響，如圖6(b)、圖7(b)和圖8(b)所示，得到的DI平滑度好，但忽略了區(qū)域異質性，對變化輪廓細節(jié)信息的丟失較為嚴重，因而邊緣表現(xiàn)不好，所對應的直方圖顯示其數(shù)據(jù)峰值分布也不太好。對于NR方法，其結合了鄰域空間信息和區(qū)域異質性測度，融合了IR方法和MR方法的優(yōu)點，去噪和細節(jié)保留也較好，但從其所生成的DI直方圖可以看出，如圖6(c)、圖7(c)和圖8(c)所示，波峰數(shù)據(jù)的分布仍然存在有不太好的現(xiàn)象，對進一步的DI分析有影響。而本文所提的INR方法中引入異質性測度來加權考慮中心像素灰度值和鄰域灰度均值信息，并以IR方法的min/max方式求得DI，從圖6(d)、圖7(d)和圖8(d)可以看出，圖像平滑性增強的基礎上仍然突出保留了變化輪廓細節(jié)，其直方圖顯示該法數(shù)據(jù)峰值分布最好，毛刺少，這將更利于后續(xù)的DI分析處理求得變化檢測圖。

(2) 變化檢測圖對比：將上面的DI分別按本文的DI分析方法得到的變化檢測二值圖，如圖9-圖11所示，定量分析的數(shù)據(jù)對比如表1-表3所示。

圖9 Bern地區(qū)變化檢測二值圖對比

圖10 Ottawa地區(qū)變化檢測二值圖對比

圖11 Sardinia地區(qū)變化檢測二值圖對比

表1 Bern地區(qū)變化檢測二值圖數(shù)據(jù)對比

表2 Ottawa地區(qū)變化檢測二值圖數(shù)據(jù)對比

表3 Sardinia地區(qū)變化檢測二值圖數(shù)據(jù)對比

由圖9可以看出，(d)的變化檢測圖輪廓和細節(jié)保留最好，與真實的手動真值圖最接近，(c)和(b)次之，(a)雖然輪廓和細節(jié)保留也不錯，但圖噪聲干擾太多，說明IR方法噪聲殘留較多。從表1數(shù)據(jù)對比來看，IR方法誤檢像素數(shù)最多，是造成其正確率低的主要原因；MR方法的漏檢數(shù)較高，是因為其采用鄰域均值比來構建DI，增強了DI的平滑性，但變化輪廓的邊緣信息丟失增多；NR方法的DI取得了不錯效果，說明其引入的鄰域均值信息和區(qū)域異質性測度是有意義的；而本文的INR方法得到變化檢測圖各項數(shù)據(jù)最好，其PCC達到99.67%，KAPPA系數(shù)達到0.866 9，且比NR方法高出1個百分點。

由圖10可以看出，(d)變化細節(jié)保留最好，尤其是上半部分的變化細節(jié)更清晰，(c)次之，(a)和(b)的漏檢情況較多。從表2數(shù)據(jù)可以看出，本文INR方法得到變化檢測圖各項數(shù)據(jù)最好，漏檢數(shù)和誤檢數(shù)均較低，PCC達到97.03%，KAPPA系數(shù)達到0.879 6，且比NR方法高了3個多百分點。對比數(shù)據(jù)均表明本文所提的INR方法來生成DI更有利于DI分析。

根據(jù)Sardinia地區(qū)的手動真值圖可以知道其變化細節(jié)要求較高，由圖11可以看出，仍然取得了較好的變化細節(jié)保留，(c)和(b)次之。從表3數(shù)據(jù)可看出，雖然幾種方法的誤檢像素數(shù)均較高，但本文INR方法得到變化檢測圖各項數(shù)據(jù)最好，誤檢數(shù)相對較低，PCC達到97.65%，KAPPA系數(shù)達到0.812 0，且比NR方法高了2個多百分點。對比數(shù)據(jù)均也表明本文所提的INR方法來生成DI更有利于DI分析。

(3) 本文方法與其他典型的幾種變化檢測方法的對比：這些典型的SAR變化檢測方法有K-means方法、RFLICM方法、SIFT方法、S-PCANet方法和MRFFCM方法，且為實現(xiàn)對比公平性，這些方法的DI構建均用本文所提的方法，其他參數(shù)設置均按其文獻闡述取值。三組SAR圖像數(shù)據(jù)的結果對比分別如圖12-圖14所示，定量分析的數(shù)據(jù)對比分別如表4-表6所示。

圖12 Bern地區(qū)六種檢測方法的變化檢測圖對比

圖13 Ottawa地區(qū)六種檢測方法的變化檢測圖對比

圖14 Sardinia地區(qū)六種檢測方法的變化檢測圖對比

表4 Bern地區(qū)六種檢測方法數(shù)據(jù)對比

表5 Ottawa地區(qū)六種檢測方法數(shù)據(jù)對比

表6 Sardinia地區(qū)六種檢測方法數(shù)據(jù)對比

從上述三組SAR圖像數(shù)據(jù)經(jīng)六種方法處理得到的變化檢測二值圖對比可以看出，在這些方法中，本文方法和S-PCANet方法表現(xiàn)最穩(wěn)定，因為這兩種方法均對DI中殘留的噪聲進行了抑制，所得到變化檢測二值圖的輪廓和變化細節(jié)均較為清晰，更趨近于手動真值圖。從表4-表6的數(shù)據(jù)對比進一步看出，本文方法的PCC值均在97%以上，KAPPA值也均大于0.81，數(shù)值居于前兩位，尤其是Sardinia地區(qū)檢測細節(jié)要求高，本文方法也獲得了不錯的結果，根據(jù)其一致性等級劃分，所得到的變化檢測結果均是有效的，趨近于手動真值圖，這表明本文設計的DI分析方法是有效的，比其他典型的幾種變化檢測方法更優(yōu)。

4 結 語

本文針對兩時相的SAR圖像變化檢測問題，提出了一種新的鄰域比方法INR來構建DI并設計了一種新的DI組合分析方法。實驗表明，通過該INR方法構建的DI能更清晰地反映兩時相圖的變化性質，能有效抑制斑點噪聲的影響，在圖像平滑性增強的基礎上仍然能突出變化輪廓細節(jié)，且使得DI的像素樣本分布更加合理，更利于DI分析。通過本文設計的DI組合分析方法，利用層次聚類獲得樣本標簽、鄰域塊特征向量作為樣本進行ELM網(wǎng)絡分類，得到變化檢測二值圖。實驗結果再次表明INR方法有效的同時，數(shù)據(jù)對比的PCC值和KAPPA系數(shù)值也說明了本文所設計的DI分析方法的有效性，其能對DI中殘留的噪聲進行抑制，變化細節(jié)保留較好，使得變化檢測精度得以提升。下一步，將繼續(xù)對DI中殘留噪聲抑制進行分析處理，考慮新的特征提取方法。