基于非負矩陣分解和極限學習機分類的合成孔徑雷達圖像變化檢測

王 平， 王宜懷， 余文森， 吳 薇

(1.武夷學院數學與計算機學院，武夷山 354300；2.認知計算與智能信息處理福建省高校重點實驗室，武夷山 354300； 3.蘇州大學計算機科學與技術學院，蘇州 215006； 4.軟件新技術與產業化協同創新中心，蘇州 215006)

合成孔徑雷達(synthetic aperture radar, SAR)圖像的變化檢測，目標是生成一個變化圖，表述兩時相或多時相圖在不同的標定時間之間的變化情況，其廣泛應用在環境監測、生態監測、災害評估、發展研究等領域[1]。SAR圖像以微波成像，其強度值與地面目標的后向散射成正比，易產生不易濾波處理的斑點噪聲，進而對變化檢測的精度有較大影響。從SAR圖像變化檢測技術流程來看，主要有圖像預處理、變化信息發現、區域提取或變化類型確定等幾個步驟，其中，變化信息發現和區域提取是變化檢測的研究熱點，目前從差異圖(difference image, DI)角度進行變化檢測是一條重要的途徑，廣泛應用的辦法是先構建兩時相SAR圖的DI，然后再對DI分析處理。

DI構建中，傳統的方法有差值法、比值法(R)[2]、改進比值法(IR)[3]、對數比值法(LR)[4]等，但這些方法均沒有利用像素點周圍的鄰域空間變化信息來減少斑點噪聲的影響，進而造成精度較低。近年來提出的組合差異法的思路是將各種方法的優點結合起來，特別是尋求噪聲抑制和圖像變化細節保留之間的均衡，諸如鄰域均值比法(MR)[5]、鄰域比法(NR)[6]等。

在DI分析中，傳統方法是直接將DI的像素點分為變化類和不變化類，Lu等[7]、Liu等[8]針對DI均采用k為2的K-means聚類；Krinidis等[9]利用改進的FLICM聚類方法對DI進行分析；王建明等[10]采用基于自適應距離的模糊C均值算法(FCM算法)；Li等[4]采用基于遺傳算法改進的FCM算法對DI進行處理；Li等[11]設計了識別DI特征空間的層次級聯FCM算法來得到變化和未變化的信息。但這些方法僅僅直接對DI像素進行聚類成兩類處理，依賴于DI構建方法的噪聲抑制，僅僅是粗粒度的分類，故而精度很難進一步提升。

近年來，無監督的增量學習模式被引入到了DI分析中，其具體步驟為：①通過閾值法或聚類法生成初始DI，其中包括高概率的變化類和不變化類，以及可視為散斑噪聲污染的不確定類像素。②從初始DI中尋求高概率的變化類和不變化類作為先驗知識(或偽標簽)，訓練算法所選擇的學習系統。③根據完成的學習系統對多時相SAR圖像進行解釋并生成最終變化檢測二值圖。Xiao等[5]將DI根據自組織映射網絡(SOM映射)為3類，即變化類、變化類和不確定類，再將深度神經網絡用于在網絡上訓練類似噪聲的訓練集，以減少變化類中的殘余噪聲并獲得最終變化圖；Mu等[12]在利用FCM將DI像素樣本聚類成3類的基礎上，利用加速遺傳算法進一步對不確定類的像素樣本進行分類而獲得變化檢測圖；Gao等[13]設計了層次FCM算法來分析DI得到樣本標簽，再將樣本結合原始圖像鄰域塊進行SVD網絡分類，進而得到變化檢測圖，對DI的噪聲殘留有一定的抑制效果；崔斌等[14]在層次FCM聚類的基礎上，利用卷積神經網絡對中間類進一步進行分類處理，但該方法對中間類處理算法復雜，耗時較長；Li等[15]提出一種顯著性檢測生成DI，并在DI的顯著性區域上提取鄰域塊特征進行聚類，并分為3類標簽，即變化類、不變化類和不確定類，選取前兩者數量均衡地構成樣本標簽，結合原始輸入兩時相SAR圖像的對應位置鄰域塊構成特征向量形成樣本集，訓練SDPCAnet，實現不確定類的最終分類，進而得到變化檢測二值圖。但這些方法在生成偽標簽時未考慮DI的鄰域空間信息，忽略了其噪聲殘留會影響標簽的正確性，進而影響檢測精度。

基于上述分析，提出在偽標簽生成時考慮DI中殘留噪聲的影響，提取DI的鄰域空間信息進行DI分析，采用非負矩陣分解(nonnegative matrix factorization, NMF)算法實現DI鄰域特征向量降維提取，并層次聚類生成較為穩定的偽標簽，再將偽標簽中的高概率的變化類和不變化類形成樣本標簽，結合輸入的兩時相SAR圖對應位置的鄰域塊組合并連接形成的特征向量樣本，訓練極限學習機(extreme learning machine, ELM)分類網絡，實現將不確定類像素進一步分類成變化類和不變化類，進而得到變化檢測二值圖，解決DI殘留噪聲影響的數據判別問題，提升檢測方法的魯棒性。

1 方法框架

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

圖1 方法的總體框架Fig.1 The overall framework of the method

方法的具體步驟如下。

(1)基于INR的DI構建。

(2)對DI采用NMF方法進行鄰域h×h特征降維提取，并進行層次FCM聚類，得到聚類差異圖，該矩陣元素值為0、0.5和1，代表樣本的偽標簽，其中，0和1分別代表DI對應像素點屬于高概率的不變化類和變化類，而0.5代表可能噪聲殘留影響的中間類。

(3)根據聚類結果的變化類和不變化類像素的偽標簽，且以這些像素點位置為中心，提取兩時相SAR圖的k×k鄰域塊并組合在一起作為鄰域塊特征向量樣本，構成訓練集，訓練ELM分類網絡學習其非線性變化關系，同時，將兩時相圖的所有k×k鄰域塊特征向量作為測試集。

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

2 相關要點原理

2.1 DI構建

DI構建采用INR方法，其基于NR思想，在IR的基礎上引入異質性測度θ，實現融合中心像素點強度和局部鄰域均值信息的線性組合，其公式為

(1)

式(1)中：DINR表示采用INR方法得到的DI；DINR(i,j)表示DINR中位置處于(i,j)位置的DI的像素強度值，DINR的維數為M×N；u1(x)、u2(x)定義為SAR圖像I1和I2的以(i,j)位置為中心鄰域空間x(大小為r×r)的強度均值；θ=σ(x)/μ(x)是鄰域空間x的異質性測度，σ(x)是鄰域x的標準差，μ(x)是鄰域x的均值；常數C解決分母可能趨近于0問題，常根據C=(kL)2取值，L為SAR圖像灰度級別，如8位灰度強度則取值為255，k一般取值為0.003。分析式(1)可知，DINR(i,j)越大，則兩時相SAR圖間的相似度越大，對應的像素點也越可能沒有變化，反之，DINR(i,j)越小則越可能發生了變化。該方法理論上可以得到較好的差異圖DI，且能實現斑點噪聲抑制和細節保留之間的折中。

另外一個問題是異質性測度值歸一化問題，式(1)中θ的理想值為[0,1]，但是當σ(x)>μ(x)時，可能出現θ>1情況，而造成DINR(i,j)取值可能為負，無法理想地表達變化關系，則異質性測度值進行統一歸一化的計算式可定義為

(2)

式(2)中：θ1,(i,j)、θ2,(i,j)分別為兩時相SAR圖以(i,j)位置為中心的鄰域x的異質性測度值；θmax為兩時相SAR圖的所有像素位置為中心的鄰域x的異質性測度的最大值；θ′1,(i,j)、θ′2,(i,j)分別是兩時相SAR圖以(i,j)位置為中心的鄰域x的歸一化異質性測度值。

2.2 DI鄰域特征的NMF降維提取并聚類

傳統的DI分析法中，大多直接對DI像素樣本直接進行聚類分析，所得的聚類結果或直接作為變化檢測二值圖結果，或作為進一步分析處理的偽標簽，但由于其沒有考慮DI中殘留噪聲的影響，精度較難進一步提升。考慮DI的鄰域特征信息，引入非負矩陣分解NMF的方法進行鄰域塊特征降維提取，再根據提取的特征進行FCM聚類分析，并按特征值聚類成3類進而得到聚類差異圖，即為變化類Ωc、中間類Ωi和不變化類Ωu。其中，變化類Ωc和不變化類Ωu分別為高概率的變化像素和不變化像素數，而中間類Ωi呈現變化與否的不確定性，其可能是噪聲殘留的影響。進一步，由于單一的FCM聚類可能會造成中間類Ωi的比重過大，使得有用信息的丟失過大，影響最終變化檢測的精度，因而采用層次FCM聚類的方法，具體方法的步驟如下。

輸入：差異圖矩陣DI。

(1)DI鄰域特征信息提取，以M×N個像素為中心，分別提取其鄰域大小為h×h的鄰域塊，并將其轉換為向量形式，最后將所有的特征向量按序組合成特征向量矩陣，其維數為2h×MN，如圖2所示。

圖2 DI特征提取示意圖Fig.2 Schematic diagram of DI feature extraction

(2)基于NMF特征向量矩陣降維，廣泛應用于圖像處理中的局部特征提取方法NMF對圖像變化細節保留方面具有很好的魯棒性，其迭代穩定后用兩個低維的非負矩陣W和H的乘積逼近原特征向量矩陣，所得的系數矩陣H保留了原始特征向量矩陣CM的本質特征的同時，可實現后續處理中占用更少的存儲空間和更少的計算資源，且能更好應用于聚類分析，能有效地抑制DI中殘留噪聲的影響。進一步，為了加速NMF的收斂速度，使用非負雙奇異值分解(non-negative double singular value decomposition, NNDSVD)方法進行分解初始化，且所分解的H矩陣的秩設置為2，進而得到低維的特征向量矩陣。

(3)根據所得到的2×MN維特征向量矩陣執行第1輪FCM聚類，將其向量聚類為上述3類：

(3)

(4)

(5)設置

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

經過層次FCM聚類后，得到了一張標有 {Ωc,Ωi,Ωu}的聚類差異圖，可令屬于Ωc的像素變化概率較大且令其值為0標簽，屬于Ωu的像素不變化概率較大且令其值為1標簽，這兩種像素可以構造下一小節ELM網絡的訓練樣本，而屬于Ωi的標簽值為0.5的像素為不確定變化與否的中間類，其進一步的標簽確定將由ELM網絡根據兩輸入時相圖的鄰域空間特征向量分類而完成，且由于特征向量來源于DI的鄰域空間，所得的偽標簽能有效抑制DI中殘留噪聲的影響。

2.3 ELM二分類

根據2.2節層次FCM聚類后的結果，且為了進一步解決Ωi的類別問題，將Ωc和Ωu作為兩輸入時相圖的鄰域特征樣本標簽構造ELM的訓練樣本集，訓練ELM網絡，再將所有鄰域特征樣本構造成ELM網絡的測試樣本集，最終可得到變化檢測圖，具體方法的步驟如下。

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

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

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

3 實驗及結果分析

3.1 數據集描述和實驗設置

實驗數據集使用4組真實且已配準預處理的SAR圖像數據，分別是Farmland_C地區數據、Ottawa地區數據、San_Francisco地區數據和Coastline地區數據的兩時相圖及手動真值圖，如圖4～圖7所示。

圖3 兩輸入時相圖鄰域特征向量的生成方法Fig.3 Method for generating neighborhood eigenvectors of two-input phase diagram

圖4 Farmland_C地區數據Fig.4 Farmland_C area data

圖5 Ottawa地區數據Fig.5 Ottawa area data

圖6 San_Francisco地區數據Fig.6 San_Francisco area data

圖7 Coastline地區數據Fig.7 Coastline area data

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

(1)將本文所提方法與方法1和方法2對比，其中：方法1不考慮DI中殘留噪聲的影響，將DI像素直接FCM聚類成變化類和不變化類像素，進而得到變化檢測二值圖；方法2是采用增量學習的模式，但與本文方法不同的地方是，雖然考慮了DI中殘留噪聲會影響檢測精度，但沒有在DI中利用鄰域空間信息來抑制噪聲進而獲得穩定的偽標簽，而是將DI像素直接進行層次FCM聚類成變化類、不變化類和不確定類像素而得到偽標簽。為不失公平性，DI構建均采用本文所提的INR方法，方法2的層次FCM聚類和ELM鄰域塊分類參數與本文相同。對比目的是驗證本文所提的考慮DI噪聲殘留，并利用NMF降維提取DI鄰域空間特征以提升偽標簽穩定性方法的有效性。

(2)將本文方法與其他幾種典型的增量學習模式的檢測方法進行對比，包括方法3的SOMDNCD方法[5]、方法4的S-aGA方法[12]、方法5的SDPCANet方法[15]。方法3采用在DI上進行鄰域空間中值濾波后進行增量學習模式的變化檢測方法；方法4在DI基礎上利用顯著性檢測模型提取顯著區域，以此實現對噪聲殘留區域進行過濾，并對顯著區域進行增量學習模式的變化檢測；方法5也是在DI基礎上利用顯著性檢測模型提取顯著區域，以此實現對噪聲殘留區域進行過濾，并在顯著區域上利用鄰域空間信息進行特征聚類生成偽標簽。為不失公平性，DI的構建均采用本文所提的INR方法，其他參數如其文獻所述。對比目的是驗證本文所提的增量學習模式變化檢測的方法與其他增量學習模式的檢測方法的優劣。

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

OE=FN+FP

(5)

PCC=(N-OE)/N

(6)

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

(7)

式中：TP為確定變化的像素且在變化檢測中被正確檢測為變化的像素數；TN為確定未變化的像素且在變化檢測中被正確檢測為未變化的像素數；則系數值計算中，Pr(a)和Pr(e)分別為

Pr(a)=(TN+TP)/N

(8)

Pr(e)=[(TP+FN)(TP+FP)+(TN+FN)(TN+FP)]/N2

(9)

KAPPA系數值越趨近于1，表明變化檢測的效果越好，根據其一致性等級劃分，一般為當其取值為0.81～1時，所得到的變化檢測結果跟手動真值圖更接近且趨于一致。

3.2 實驗結果及分析

(1)Farmland_C地區數據的變化檢測結果如圖8所示，定量分析的結果如表1所示。該數據集變化區域的像素值相對較小，且兩時相圖之間不變化區域之間的像素強度對比較大，使得構建的DI變化區域與不變化區域的對比度較低，因而對檢測算法的要求較高。從圖8中可以看出，圖8(a)的主變化區域的邊緣漏檢情況較多，且由于噪聲殘留干擾的誤判也比較多；圖8(b)可看出，增量學習模式檢測方法對檢測效果有明顯提升，但DI中噪聲殘留的干擾使得誤判也比較多，這主要是由于處理過程中DI中的噪聲殘留影響了層次聚類生成偽標簽的穩定性；圖8(c)對DI進行中值濾波抑制噪聲后，進行增量學習模式地變化檢測，從其二值圖可看出，變化區域的檢測效果較好，但中值濾波的方法對噪聲殘留區域的抑制有限；圖8(d)對DI進行顯著性區域提取后再針對該區域進行增量學習模式的變化檢測，可看出其可以有效地提升變化檢測精度；圖8(e) 對DI進行顯著性區域提取后再針對該區域提取鄰域空間信息進行特征聚類生成偽標簽，相對于圖8(d)，可看出其可以進一步提升變化檢測精度，說明對DI鄰域空間信息進行處理是有意義的，但由于其只針對顯著區域進行處理，對其他不變化區域的有用信息有大的丟失，在不變化區域反而增強了殘留噪聲的影響；圖8(f)是本文所提方法所得的二值圖，可以看出其取得了最好的變化效果，雖然右下角有局部變化信息丟失，但從輪廓細節上表現最好。從表1的定量分析數據可以看出，屬于增量學習模式的檢測方法漏檢數均較低，這些方法的PCC均在98%以上，且KAPPA值也在0.81以上，說明增量學習模式的檢測方法有助于提升變化檢測精度，且根據一致性等級劃分，這些檢測的結果都是有效的。

在這些方法里，本文方法考慮了DI的鄰域空間信息來抑制DI中殘留噪聲的影響，取得了最好的變化檢測數據，錯檢數降到103，噪聲干擾殘留也最少，其KAPPA與其他對比方法相比高2個百分點以上，這些數據對比也說明本文方法在DI分析中考慮DI的鄰域空間信息是有效的，由其生成的偽標簽是最穩定的，有助于提升變化檢測效果。

(2)Ottawa地區數據的變化檢測結果如圖9所示，定量分析的結果如表2所示。該數據集的兩時相圖成像斑點噪聲影響較小，且圖像不變化區域像素強度值變化也小，故而對算法的要求不高，同樣，從圖9整體可以看出，增量學習模式的檢測方法均取得不錯的效果，說明了這些檢測方法的有效性。進一步，圖9(f)本文方法和圖9(e)方法5、圖9(d)方法4從二值圖上看更接近于手動真值圖，尤其是上半部分的變化細節保留最好。從表2的數據對比看，屬于增量學習模式的檢測方法的KAPPA值也均大于0.86以上，也說明了這些檢測方法的有效性，其中，方法3由于對DI采用中值濾波后層次聚類得到偽標簽，該法所得的偽標簽的穩定性低于方法4、方法5和本文方法，其對DI殘留噪聲的抑制有限，故而其漏檢數較大，是造成KAPPA低于方法4、方法5和本文方法的原因，這點從該方法所得的二值圖的上部即可看出，本文方法和方法5的PCC值達到98%以上，KAPPA也達到0.93以上，這兩種方法均對DI有提取鄰域空間特征進行噪聲抑制，這說明提取DI鄰域空間特征是有意義的。

圖9 Ottawa地區變化檢測二值圖對比Fig.9 Comparison of binary maps for change detection in Ottawa area

另外，本文方法比方法2的變化檢測結果PCC高1.5個百分點，KAPPA也高6個百分點以上，這也說明對DI有提取鄰域空間特征進行噪聲抑制而得到增量學習的偽標簽更穩定，能有效抑制DI中殘留噪聲的影響，有助于提升變化檢測的精度。

表2 Ottawa地區變化檢測二值圖數據對比

(3)San_Francisco地區數據的變化檢測結果如圖10所示，定量分析的結果如表3所示。該數據集的噪聲干擾與Farmland_C地區類似，也對變化檢測算法的要求較高。圖10(a)方法1的直接聚類得變化二值圖，漏檢情況較多；圖10(b)方法2的效果明顯比方法1有了較大提升，但同樣DI噪聲殘留影響的誤判也比較多；從圖10(c)方法3開始到圖10(f)本文方法的變化二值圖可看出，有效地

圖10 San_Francisco地區變化檢測二值圖對比Fig.10 Comparison of binary maps for change detection in San_Francisco area

表3 San_Francisco地區變化檢測二值圖數據對比

抑制DI中殘留噪聲的影響能更進一步提升檢測效果，且其中圖10(f)本文方法取得了最好的檢測效果。從表3數據對比看，屬于增量學習模式的方法2、方法3、方法4、方法5和本文方法均取得不錯效果，其PCC均大于98%，KAPPA也達0.83以上，從一致性等級劃分看，也都是有意義的，其中，表現最好的是本文方法，PCC達99%，KAPPA達0.925 6，高于方法4、方法5近3個百分點，說明DI鄰域空間特征的提取對增量學習模式中的偽標簽的生成具有較強的穩定性，對DI噪聲殘留有較好的抑制作用。

(4)Coastline地區數據的變化檢測結果如圖11所示，定量分析的結果如表4所示。該數據集兩時相圖的兩時相圖成像斑點噪聲影響較大，且之間的像素強度值對比也比較大，這為變化檢測的算法也有較高的要求。同樣可以看出：圖11(a)方法1誤判情況和漏判情況均較多；圖11(b)方法2開始到圖11(f)本文方法的各種增量學習模式的檢測方法中，檢測效果均遠遠高于方法1，說明了增量學習模式的有效性；圖11(c)方法3到圖11(f)本文方法效果又比方法2有了明顯提升，這說明對DI中殘留噪聲的抑制是有意義的，能增強增量學習模式中偽標簽的穩定性。這些方法里圖11(e)方法5和圖11(f)

圖11 Coastline地區變化檢測二值圖對比Fig.11 Comparison of binary maps for change detection in Coastline area

表4 Coastline地區變化檢測二值圖數據對比

本文方法均采用了DI鄰域空間特征提取來抑制噪聲，取得了最好的檢測效果。從表4的數據對比也可以看出，本文方法的PCC達99.81%，KAPPA接近0.91，在增量學習模式的方法中取得了最好的效果。

由上述4組真實的SAR數據集圖像的變化檢測可以看出，增量學習模式的變化檢測方法，有助于提升變化檢測的精度，但DI中殘留噪聲會影響增量學習模式的偽標簽的穩定性，有效地抑制這些殘留噪聲將會進一步提升增量學習模式的變化檢測方法的精度。這4組真實數據中，本文所提的對DI采用鄰域空間特征降維提取來抑制噪聲殘留影響偽標簽穩定性的方法，有3組數據均相對于其他增量學習方法取得了最好的效果，Ottawa地區數據雖然檢測的效果略低于方法5，但也遠遠優于其他增量學習方法，這也可以看出本文方法針對DI中噪聲殘留較多的數據集效果更優。

4 結論

針對兩時相的SAR圖像變化檢測問題，提出了一種ELM網絡分類的增量學習模式的變化檢測方法，并為了提升偽標簽的穩定性，該方法在對DI分析處理中，以提取DI的鄰域空間特征信息來抑制其噪聲殘留的影響，進而可使得層次聚類得到的偽標簽具有高穩定性。實驗表明，當變化檢測的輸入兩時相SAR圖的噪聲干擾較為嚴重時，該方法能有效抑制噪聲，使得最終的變化檢測二值圖更趨近于手動真值圖，能使變化細節保留較好，變化檢測精度能得到有效提升，說明本文所提方法的有效性。下一步，將繼續對DI中殘留噪聲抑制進行分析處理。