基于地基雷達圖像的無監督變化檢測

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

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

(內蒙古自治區雷達技術與應用重點實驗室 呼和浩特 010051)

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

2.1 初始差異圖獲取

由于地基雷達圖像自身存在乘性斑點噪聲，對數比值法能夠將乘性斑點噪聲轉化為加性噪聲，便于有效去除斑點噪聲。此外，對數運算作用于兩幅地基雷達圖中能夠對像素的灰度值的比值做非線性拉伸，可以在一定程度上增強變化和非變化區域的對比。為了充分考慮像素點周圍的特征，采用均值對數比獲取的差異圖為Xd，計算公式為

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

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

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

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

2.2 NSCT變換

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

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

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

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

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

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

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

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

2.3 PCA提取特征向量

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

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

將所有列向量重構為32×n的矩陣S，計算S的協方差矩陣C

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

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

2.4 改進的FCM聚類

2.4.1 FCM算法

FCM是一種無監督的模糊聚類方法，傳統FCM[11]目標函數的形式為

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

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

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

聚類中心Ci為

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

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

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

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

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

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

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

2.4.2 改進的FCM算法

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

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

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

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

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

R為最終變化檢測結果圖。

3 實驗結果與分析

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

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

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

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

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

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

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

圖6 聚類結果對比Fig.6 Clustering result comparison

圖7 監測現場光學圖Fig.7 Optical image of monitoring site

4 結束語

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

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