基于地基雷達圖像的無監(jiān)督變化檢測

黃平平 任慧芳 譚維賢 段盈宏 徐 偉 劉 方

(內(nèi)蒙古工業(yè)大學信息工程學院 呼和浩特 010051)

(內(nèi)蒙古自治區(qū)雷達技術(shù)與應(yīng)用重點實驗室 呼和浩特 010051)

1 引言

我國地質(zhì)構(gòu)造復(fù)雜、地形地貌起伏變化大，極易發(fā)生滑坡、崩塌、泥石流等地質(zhì)災(zāi)害。此外，人類的生產(chǎn)、生活活動也很容易誘發(fā)滑坡災(zāi)害的發(fā)生，大量巖土邊坡的失穩(wěn)或滑動屢見不鮮，給工程建設(shè)和人民生命財產(chǎn)帶來了嚴重的危害和損失。災(zāi)害發(fā)生區(qū)域的變化檢測分析對安全生產(chǎn)有著重要意義。從遙感數(shù)據(jù)中獲取災(zāi)害信息是一種十分重要的災(zāi)害研究手段，而快速有效的變化檢測方法對于掌握災(zāi)害情況尤為重要。地基雷達技術(shù)具有全天時、全天候、大范圍、遠距離、非接觸等特點，在地質(zhì)災(zāi)害監(jiān)測、預(yù)警方面具有巨大優(yōu)勢。目前，地基雷達已廣泛應(yīng)用于滑坡、崩塌等地質(zhì)災(zāi)害的監(jiān)測，利用干涉測量原理地基雷達可以監(jiān)測到目標區(qū)域發(fā)生的微小形變量，根據(jù)累計形變量、形變速度、形變加速度等對地質(zhì)災(zāi)害進行實時監(jiān)測。然而，在長期施工建設(shè)監(jiān)測中，部分區(qū)域的變化并不會引起滑坡等災(zāi)害，但受人為因素、地質(zhì)因素、氣象因素等影響，導(dǎo)致監(jiān)測區(qū)域的雷達圖像失相干較為嚴重，給長期定量化監(jiān)測帶來了較大的難度。因此，迫切需要在定量監(jiān)測的基礎(chǔ)上，進一步開展變化檢測方面的應(yīng)用，為長期全面了解監(jiān)測區(qū)域的動態(tài)變化提供有效信息。針對上述問題，本文將研究一種基于地基雷達圖像的變化檢測方法。

一般來說，圖像變化檢測技術(shù)可以分為兩大類：監(jiān)督式和無監(jiān)督式。監(jiān)督變化檢測方法，需要先驗信息；而無監(jiān)督變化檢測方法直接自動比較多時相的遙感圖像，不需要任何額外信息，降低人為誤差的影響，符合實際應(yīng)用中先驗變化信息缺失的現(xiàn)實情況，因此得到了廣泛地研究[1]。雷達圖像中的無監(jiān)督變化檢測一般包括3個步驟：(1)圖像預(yù)處理；(2)獲取差分圖像；(3)分析差異圖像及后處理[2,3]。

近年來雷達圖像變化檢測的研究主要集中在機載和星載合成孔徑雷達(Synthetic Aperture Radar,SAR)圖像上，基于星載SAR圖像中的無監(jiān)督變化檢測，龔茂國等人[4]提出了一種基于鄰域的比值法來產(chǎn)生差異圖，雖然該算法優(yōu)于產(chǎn)生差異圖的傳統(tǒng)方法，如差值法和比值法，但是其檢測結(jié)果受噪聲干擾嚴重；為了提高變化檢測結(jié)果的準確率，使用均值比圖像和對數(shù)比圖像的互補信息來生成差異圖，取得了較好的檢測結(jié)果[5]。Sumaiya等人[6]采用了一種經(jīng)濟且簡單的對數(shù)均值閾值進行星載SAR圖像變化檢測。楊祥立等人[7]運用D-S證據(jù)理論(Dempster-Shafer evidence theory,D-S)融合高分辨率SAR影像的相干/非相干差異特征進行變化檢測，采用美國舊金山地區(qū)港口附近的TerraSAR-X數(shù)據(jù)切片進行實驗驗證。此外，在圖像變化檢測中一個重要的步驟是圖像的分割聚類，Elguebaly等人[8]提出了一種基于廣義高斯混合模型(Generalized Gaussian Mixture models,GGM)的高效無監(jiān)督圖像分割和變化檢測算法，論證了貝葉斯GGM在圖像分割中的應(yīng)用潛力。邵寧遠等人[9]提出了一種面向變化檢測的SAR圖像超像素協(xié)同分割算法，解決了變化檢測中存在的多時相圖像邊界和空間對應(yīng)關(guān)系不一致的問題。Carincotte等人[10]提出了一種新的模糊隱馬爾可夫鏈(Hidden Markov ChainS,HMCS)，從而解決了模糊變化檢測的統(tǒng)計方法。而最著名的模糊聚類算法是模糊C均值聚類(Fuzzy C-Means,FCM)分割算法[11]，它主要的缺點是在成像應(yīng)用中缺乏對空間環(huán)境信息的整合，導(dǎo)致其對噪聲和其它人為因素較為敏感。而FCM的變形，如魯棒FCM(Robust FCM,RFCM)[12]和空間條件FCM(the Spatial Conditional FCM,SCFCM)[13]在星載SAR圖像聚類中都得到了較好的效果。

星載SAR適用于大范圍對地成像，其成像范圍可達數(shù)千平方公里，基于星載SAR圖像變化檢測適用于監(jiān)測區(qū)域內(nèi)發(fā)生的大面積、特征明顯的變化。地基雷達通常針對特定的目標區(qū)域成像，多數(shù)情況下發(fā)生變化的特征不夠明顯，將星載SAR圖像變化檢測技術(shù)應(yīng)用到地基雷達圖像變化檢測中，會造成差異圖群分性比較差和聚類分割結(jié)果存在較多噪聲點等缺陷。本文提出了一種無監(jiān)督地基雷達圖像變化檢測方法，本方法通過對相干系數(shù)圖和均值對數(shù)比值圖進行非下采樣輪廓波變換(NonsubSampled Contourlet Transform,NSCT)和局部能量法得到合成差異圖，得到的合成差異圖具有較好的群分性。然后采用主成分分析法(Principal Component Analysis,PCA)提取特征向量，增強抗噪性。最后使用改進的模糊C均值聚類(FCM)對特征向量進行聚類得到最終的變化檢測結(jié)果。文章安排如下：本文第2部分提出了一種地基雷達圖像變化檢測方法；第3部分利用地基雷達實測數(shù)據(jù)進行結(jié)果分析并驗證該方法的有效性；第4部分對全文做出總結(jié)。

2 地基雷達圖像變化檢測算法

本文所提基于地基雷達圖像的無監(jiān)督變化檢測方法主要分為4個步驟：

(1) 初始差異圖獲取。先取兩幅已配準的地基雷達圖像X1與X2，通過相干系數(shù)算子以及均值對數(shù)比算子生成相干系數(shù)差異圖和均值對數(shù)比差異圖；

(2) NSCT變換。對得到的兩幅初始差異圖，進行NSCT變換并融合低頻系數(shù)，然后使用融合的低頻系數(shù)分別與初始差異圖的高頻系數(shù)進行NSCT反變換得到兩幅融合結(jié)果圖，對這兩幅融合結(jié)果圖進行像素能量合成得到最終的合成差異圖；

(3) PCA提取特征向量。利用PCA算法提取出合成差異圖的像素特征；

(4) 改進的FCM聚類。最后通過一種加相似性懲法項的模糊C-均值聚類算法(FCM)聚類，獲取變化檢測結(jié)果。

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

圖1 地基雷達圖像無監(jiān)督變化檢測流程圖Fig.1 Flow chart for unsupervised change detection based on ground-based radar image

2.1 初始差異圖獲取

由于地基雷達圖像自身存在乘性斑點噪聲，對數(shù)比值法能夠?qū)⒊诵园唿c噪聲轉(zhuǎn)化為加性噪聲，便于有效去除斑點噪聲。此外，對數(shù)運算作用于兩幅地基雷達圖中能夠?qū)ο袼氐幕叶戎档谋戎底龇蔷€性拉伸，可以在一定程度上增強變化和非變化區(qū)域的對比。為了充分考慮像素點周圍的特征，采用均值對數(shù)比獲取的差異圖為Xd，計算公式為

其中，m1(i,j)與m2(i,j)分別為圖像X1和X2在像素點(i,j)處的3×3鄰域均值。

雷達圖像之間的相干系數(shù)用來衡量影像之間的相似程度，它提供了檢測區(qū)域重要的變化信息。采用相干系數(shù)圖，可以得到變化區(qū)域非常明顯和背景信息較為平滑的差異圖[14,15]。圖像X1和X2的相干系數(shù)定義如式(2)所示

式中，*表示共軛復(fù)數(shù)，X1和X2為經(jīng)過配準濾波的地基雷達圖像，相干系數(shù)Xm(i,j)的值在0～1之間，0表示失相干，1表示完全相干。

采用對數(shù)比值圖和相干系數(shù)圖相結(jié)合的方法，綜合利用各自的優(yōu)勢，能夠得到群分性較好的融合差異圖，為后續(xù)的圖像聚類分割提供重要依據(jù)。

2.2 NSCT變換

NSCT是基于非下采樣金字塔濾波器組(Non Subsamlped Pyramid Filter Banks,NSPFB)和非下采樣定向濾波器組(Non Subsampled Directional Filter Banks,NSDFB)進行的，每個子帶圖像與原始圖像在形狀和大小上都是完全一樣的，因此，NSCT具有多分辨率、局部化、方向性、各向異性和位移不變性等特點[16]。該方法能較好地捕捉圖像的幾何細節(jié)，保持目標信息和邊緣信息[17]。

基于NSCT的圖像融合算法通常將圖像進行NSCT分解，針對不同頻率域的特點選擇不同的融合規(guī)則。對數(shù)比值圖Xd和相干系數(shù)圖Xm的融合結(jié)果圖為Fd，具體實現(xiàn)步驟如下：

(1) NSCT變換。輪廓波變換使用“9-7”型拉普拉斯金字塔濾波器，使用“pkva”作為方向組濾波器，每個金字塔級方向濾波器組分解層數(shù)為4。經(jīng)NSCT變換得到差異圖Xd和Xm對應(yīng)的低頻系數(shù)Ld和Lm，高頻系數(shù)Hd和Hm，高、低頻系數(shù)反映了高、低頻圖像信息。

(2) 加權(quán)平均法。由于圖像的低頻信息通常反映圖像的結(jié)構(gòu)信息[18]，因此在融合過程中，為了凸顯變化區(qū)域，針對低頻系數(shù)采取加權(quán)平均的方式，對差異圖Xd和Xm經(jīng)NSCT變換得到的低頻系數(shù)Ld和Lm進行融合，得到融合后的低頻系數(shù)Lf，計算過程為[19]

其中，α(α ≥0)為加權(quán)平均的融合系數(shù)，本文實驗中α=0.3。

(3) NSCT反變換。用融合的低頻系數(shù)Lf分別與差異圖Xd和Xm的高頻系數(shù)進行NSCT反變換，得到融合結(jié)果圖Fd和Fm[19]。

因為地基雷達圖像的相干系數(shù)圖的背景相對平坦，可以采用局部能量法抑制背景信息[5]，得到最終的合成差異圖Fh為

其中，Ed(i,j)和Em(i,j)分別表示差異圖Fd和Fm的局部能量，Li,j表示以像素(i,j)為中心，大小為3×3的鄰域，F(xiàn)d(q)和Fm(q)表示鄰域內(nèi)第q個像素值[10]。

2.3 PCA提取特征向量

主成分分析(PCA)，也稱Karhunen-Leve變換(簡稱K-L變換)，是一種線性變換。它是一種均方誤差最小的最佳正交變換，是在統(tǒng)計特征基礎(chǔ)上的多維線性變換。本文采用PCA方法對融合差異圖進行特征提取，由于采用分塊的思想，使得該方法具有一定的抗噪能力[20]。

將合成差異圖分為n個大小為3×3的相鄰的小塊，小塊之間互不重疊。每個小塊可以看作一個矩陣，然后將每個小塊矩陣轉(zhuǎn)換為列向量Pt，其中t=1,2,…,n，并計算全部列向量的均值向量

將所有列向量重構(gòu)為32×n的矩陣S，計算S的協(xié)方差矩陣C

其中，T 表示矩陣轉(zhuǎn)置。協(xié)方差矩陣C是一個大小為32×32的矩陣。對協(xié)方差矩陣進行特征值分解，求出特征值和特征向量，按特征值從大到小的順序排列，選出對應(yīng)的特征向量，這些特征向量形成一組正交基，取特征向量的前r個列向量構(gòu)成矩陣A,A是一個大小為32×r的矩陣。取合成差異圖Fh中每個像素所在的3×3鄰域小塊，將每個小塊轉(zhuǎn)化為列向量Yη,η=1,2,…,k,k為合成差異圖Fh中像素的總個數(shù)。將所有的列向量映射到矩陣A上，即

其中，Yη為每個像素代表的特征向量，AT為A的轉(zhuǎn)置矩陣，這樣每個像素都用一個R維的特征向量來表示，[Y1Y2…Yk]構(gòu)成特征向量空間Q，Q即為一個r×k矩陣。

2.4 改進的FCM聚類

2.4.1 FCM算法

FCM是一種無監(jiān)督的模糊聚類方法，傳統(tǒng)FCM[11]目標函數(shù)的形式為

其中，uηi是像素點Yη在i類的隸屬度矩陣，c為分類數(shù)，m為模糊權(quán)重(一般為2),η為圖像的像素點個數(shù)，Ci為聚類中心。

uηi表示Yη隸屬于第i類的概率，介于0和1之間。因此，約束必須滿足

考慮到式(10)中的約束，通過拉格朗日數(shù)乘算子計算當式(9)中的目標函數(shù)達到最小值時，uηi為

聚類中心Ci為

當用于圖像聚類時，F(xiàn)CM沒有考慮到相鄰圖像像素之間的空間關(guān)系[21]。采用變形的FCM為魯棒模糊C-均值(RFCM)[12]的目標函數(shù)采取了如式(13)所示形式

其中，Nη是像素點Yη的空間鄰域集合，Pi是除i類外的所有類的集合，β是控制附加空間懲罰項效果的正偏差。

增加的空間懲罰項鼓勵像素點在同一類中擁有有較高的聚類資格。最小化式(13)得到隸屬度函數(shù)為

Ci的計算過程為式(12)，因為懲罰項不依賴于它。

FCM的另一個變型也在其目標函數(shù)中包含了空間關(guān)系，它是在目標函數(shù)中引入的空間條件FCM(也稱SCFCM)[13]。具體來說，它的目標函數(shù)采用如式(15)的形式

最小化式(15)得到隸屬度函數(shù)為

Ci的計算過程為式(12)。

2.4.2 改進的FCM算法

為了彌補FCM和其變形RFCM和SCFCM分割率低以及對于地基雷達圖像噪聲敏感等問題。本文對FCM的目標函數(shù)進行改進，在原始的FCM的目標函數(shù)中加入相似性懲罰項。目的是使圖像中相類似的像素點以更大概率聚類到相同的類，提高聚類效果并抑制噪聲點。改進后的聚類算法的目標函數(shù)的表達式為

其中，Yη為輸入每個像素的特征向量，分類數(shù)c=2(變化類和非變化類)，m模糊權(quán)重(一般為2),β控制懲罰項的積極因子，本文中β=0.2615。Nη為Yη及其周圍鄰域的一個集合，l為此集合中的一個元素。

改進的目標函數(shù)在保留傳統(tǒng)FCM目標函數(shù)的基礎(chǔ)上，添加了相似性懲罰項。需要強調(diào)的是，在式(17)的鄰域中，鄰域是以相似性而不是以空間定義的，這樣可以找到每個輸入圖像的最相似鄰域，以鼓勵相似的像素以更大的概率聚類到同一類中。通過拉格朗日數(shù)乘算子計算當式(17)中的目標函數(shù)達到最小值時，uηi為

由于相似度懲罰項不包含Ci，聚類中心采用式(12)中給出的FCM算法的相同形式。改進的FCM聚類算法計算過程如表1所示。

通過比較每個像素的隸屬度uηi的大小，uη1和uη2分別代表第η個像素對應(yīng)兩類隸屬度，采用如式(19)所示的判決條件對uη1和uη2進行二值化處理

R為最終變化檢測結(jié)果圖。

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

我國西南某省發(fā)生地震之后山體嚴重垮塌，形成堰塞體，為確保堰塞體整治工程安全、順利進行，在現(xiàn)場部署地基雷達LSA對堰塞體治理施工過程進行監(jiān)測。經(jīng)現(xiàn)場調(diào)研，地基雷達LSA布置在距堰塞體300 m正對面，雷達的監(jiān)測范圍完全可以有效覆蓋堰塞體。表2為地基雷達LSA系統(tǒng)參數(shù)，圖2所示為監(jiān)測現(xiàn)場光學圖。監(jiān)測期間，施工地出現(xiàn)持續(xù)性強降水，導(dǎo)致邊坡發(fā)生滑坡。如圖3(a)為滑坡前雷達圖，數(shù)據(jù)獲取時間為2019年7月15日，圖3(b)為滑坡后雷達圖，數(shù)據(jù)獲取時間為2019年7月19日。截取滑坡前后數(shù)據(jù)切片對進行實驗，如圖3(c)和圖3(d)所示。

表1 改進的FCM算法計算步驟Tab.1 Improved FCM algorithm calculation steps

圖4 初始差異圖Fig.4 Initial difference images

采用均值對數(shù)比值算子和相干系數(shù)算子對圖3(c)和圖3(d)的實驗數(shù)據(jù)進行處理，得到均值對數(shù)比值圖和相干系數(shù)圖，如圖4所示。將相干系數(shù)圖和對數(shù)比值圖經(jīng)過NSCT處理融合(方法1)產(chǎn)生的合成差異圖Fh與通過鄰域均值比和對數(shù)比值圖經(jīng)過NSCT處理融合(方法2)產(chǎn)生的差異圖進行對比，如圖5(a)和圖5(b)所示，其中colorbar代表圖像像素值域，方法2為星載SAR圖像變化檢測生成差異圖最常采用的方法[19]。將兩種方法產(chǎn)生的差異圖顯示在3維圖中，3維圖如圖5(c)和圖5(d)所示，X,Y軸代表圖像像素的位置信息，Z軸代表差異度值越大發(fā)生變化的可能性越大，可以明顯發(fā)現(xiàn)通過相干系數(shù)圖和對數(shù)比值圖經(jīng)過NSCT處理融合的差異圖群分性最好，這就代表了本文產(chǎn)生差異圖的方法更容易進行聚類分割，進而區(qū)分變化區(qū)域與非變化區(qū)域。

對兩時相的地基雷達圖像進行變化檢測，本文方法產(chǎn)生的合成差異圖采用改進的FCM算法與傳統(tǒng)的FCM算法的變形(包括RFCM和SCFCM)進行聚類結(jié)果對比，由于采用FCM算法進行聚類，無法得出聚類結(jié)果，在本文中不作分析。其結(jié)果如圖6所示，圖6(a)為RFCM聚類結(jié)果，圖6(b)為SCFCM聚類結(jié)果，圖6(c)為采用本文改進的FCM聚類算法得到的變化檢測結(jié)果，即式(19)的變化檢測結(jié)果R，其中藍色的像素值為0，表示變化區(qū)域；黃色的像素值為1，表示未變化區(qū)域。

圖6(a)中可以明顯地看到，RFCM算法聚類結(jié)果含有較多噪聲點，無法聚類出結(jié)果。圖6(b)中SCFCM算法去除較多噪聲點，但依然有少量噪聲存在，如圖6(b)中紅色矩形區(qū)域。而圖6(c)中本文算法很好地濾除了周圍的噪聲點，得到的變化檢測結(jié)果更加準確。通過實地考察，證實實驗數(shù)據(jù)對應(yīng)區(qū)域發(fā)生明顯的滑坡，圖7所示為監(jiān)測現(xiàn)場光學圖，圖7(b)中黃色勾選區(qū)域發(fā)生了明顯的滑坡。此外，圖6(c)和圖6(b)相比，橢圓紅框中檢測到的變化區(qū)域存在較大差異，本文方法檢測到該部分并沒有發(fā)生明顯的變化區(qū)域，在滑坡現(xiàn)場實地考察過程中發(fā)現(xiàn)滑坡前后此處為地質(zhì)條件較為穩(wěn)定的巖石，并沒有隨著滑坡的發(fā)生有所變化，只有兩塊巖石中間的部分小區(qū)域發(fā)生浮土下滑，如圖7(b)中最小的黃色框區(qū)域，由此可確定檢測到的變化結(jié)果與實地考察結(jié)果相符。地基雷達圖像變化檢測結(jié)果與地形3維數(shù)據(jù)融合結(jié)果如圖8所示，根據(jù)融合圖，可以準確地定位變化區(qū)域，如圖8中紅色框區(qū)域，通過融合圖，進一步驗證了本文方法的有效性，在實際應(yīng)用中，可以為安全生產(chǎn)提供有效監(jiān)測數(shù)據(jù)。

圖5 差異圖對比Fig.5 Difference images comparison

圖6 聚類結(jié)果對比Fig.6 Clustering result comparison

圖7 監(jiān)測現(xiàn)場光學圖Fig.7 Optical image of monitoring site

4 結(jié)束語

圖8 變化檢測結(jié)果與地形3維數(shù)據(jù)融合圖Fig.8 Effect of fusion of change detection results and UAV tilt photography

對于地基雷達圖像采用均值比以及對數(shù)比都無法得到群分性較好的差異圖，本文提出使用相干系數(shù)與均值對數(shù)比得到群分性較好的差異圖，使得地基雷達圖像在變化檢測過程中更容易進行聚類分割。此外，本文對傳統(tǒng)FCM聚類算法進行改進，使圖像中相類似的像素點以更大概率聚類到相同的類，以提高聚類效果并抑制噪聲點。通過將所提出的方法應(yīng)用于地基雷達實測數(shù)據(jù)，并將改進的FCM與FCM的另外兩種變形RFCM和SCFCM聚類算法進行聚類結(jié)果對比，結(jié)果表明，本文所提方法具有更好的聚類效果，并且變化檢測結(jié)果在保留變化區(qū)域的同時噪聲點明顯減少。未來的工作中，擬進一步探索地基雷達圖像的無監(jiān)督變化檢測方法，以得到更為精確的變化檢測結(jié)果，為安全生產(chǎn)、排險、救援提供數(shù)據(jù)參考。