基于多源遙感時空譜特征融合的滑坡災害檢測方法

陳善靜 向朝參 康 青 吳 韜 劉 凱 馮 亮 鄧 濤

1(陸軍勤務學院 重慶 401311) 2(重慶大學計算機學院 重慶 400044) 3(國防科技大學電子對抗學院 合肥 230037) 4(78092部隊 成都 610031)chengshanjing_11@163.com)

我國約70%的地域為山區，地質災害發生頻率高，每年因各種地質災害造成的損失數以億計，給國家和人民帶來巨大損失.我國西南地區多山地，地質、地貌穩定性差，在地震或暴雨等誘發因素作用下滑坡等地質災害可能呈大面積高發狀態.在大面積滑坡災害發生時快速獲取滑坡區域分布、數量、規模等災情信息對救援決策和防災減災都有著重要意義.當前針對滑坡的檢測識別與災情信息提取多以現地勘測為主，存在耗時長、工作量大、效費比低，提取災情信息不直觀、不精確、不全面，受人為因素和地理環境影響大，難以對災害地區實現快速大面積精確檢測與綜合評估等問題.

在航空遙感滑坡檢測識別方面，文獻[8]提出了一種基于航空遙感圖像的半主動滑坡危險區域制圖方法，該方法利用變化向量分析和閾值分割提取候選滑坡危險區域，以形態學方法消除誤差，最后以邊緣分級進化和區域分級進化對滑坡目標進行快速識別.該方法主要適用于高分辨率航空正射影像，而對多類型、多譜段、低分辨率遙感圖像其適用效果不夠理想.文獻[9]以2013年蘆山地震震區發生的滑坡為研究對象，通過構建區域地質特色的地震滑坡樣本庫，引入遷移學習機制實現基于無人機高分辨率遙感影像地震滑坡信息提取.該方法對地震滑坡識別精度較高，滿足地震滑坡災害環境宏觀調查、滑坡災害體監測等應急需求，但是在訓練樣本完整性、全面性和特征準確提取方面還有待進一步完善和優化.文獻[10]以機載激光雷達(light detection and ranging, LiDAR)提取的DEM數據和QuickBird衛星圖像為基礎，通過小波變換對2類數據進行融合，并根據目標光譜、紋理和局部空間位置特征對其進行面向對象的目標分類，最終實現了對熱帶城市地區滑坡的精確檢測，但該方法對LiDAR遙感數據的空間分辨率要求較高，而低分辨率遙感圖像融合處理效果不夠理想.

本文從大面積滑坡發生區域災情信息快速獲取的客觀需求出發，以災前災后多波段光學遙感圖像為基礎，在綜合分析滑坡區域災前災后時空譜特征變化規律的基礎上，提出了一種基于多源遙感時空譜特征融合的滑坡災害檢測方法，通過對滑坡發生前后遙感圖像進行光譜空間和尺度空間配準、多特征融合遙感影像數據構建、基于SVM(support vector machine)的地物目標識別、典型形狀特征分類，實現對滑坡目標的精確識別.

1 基于多維特征融合的遙感影像集構建

1.1 多源遙感影像光譜空間配準

Fig. 1 Results of spectral space registration for remote sensing images(remote sensing images of Sentinel-2)圖1 遙感圖像進行光譜空間配準后結果(Sentinel-2衛星遙感影像)

在空天遙感成像中各種光學成像設備的光譜輻射分辨率各不相同，獲取的遙感圖像定標和校正方法也存在一定差異，這就造成了多源光學遙感圖像在多個譜段上的輻射亮度分布呈現較大差異.即使同一衛星對相同場景多次成像其獲取的光學圖像也可能出現較為明顯的光譜特征變化.以Sentinel-2衛星遙感成像為例，其對同一區域滑坡前后遙感成像獲得的多光譜圖像就存在明顯的色差，如圖1(a)(b)所示.滑坡后遙感圖像各地物顏色與真實目標較為接近，但滑坡前的遙感圖像整體偏綠.同一地物目標在災前災后圖像的顏色或光譜特征差異將直接影響后期數據融合處理和目標檢測應用的準確性和可靠性，因此，本文在深入分析滑坡災前災后遙感圖像光譜特征的基礎上提出對同一場景多時相多源遙感圖像進行光譜空間分布配準.對于第p波段的圖像，其直方圖分布概率函數可表示為

(1)

其中，sk為圖像第k級灰度值，0≤sk≤1，k=0,1,…,L-1，np,k為第p波段圖像中灰度值為k的像素的個數，L為圖像灰度級，N為該波段圖像中像素的總數.直方圖提供了特定譜段遙感圖像中各種灰度分布情況，是一幅圖像所有灰度或亮度值的整體描述[11].將直方圖分布轉換為積累分布函數：

(2)

Fig. 2 The histogram distribution of unregistrated/ registrated image for landslide area圖2 配準前后滑坡區域的遙感圖像直方圖分布情況

本文在多源遙感影像光譜空間配準中以災后滑坡區域各波段遙感圖像直方圖分布函數為基準，對災前遙感圖像進行對應的直方圖相似變換.根據單映射規則從小到大依次找到能使式(3)最小的k和l，將災前圖像直方圖概率分布Qp,before(sk)對應到災后圖像直方圖概率分布Qp,after(sl).

minψ(k,l)=|tp,before(k)-tp,after(l)|,

(3)

其中，tp,before(k)和tp,after(l)為滑坡區域災前和災后p波段圖像的直方圖積累分布函數，k∈[0,Lbefore-1]，l∈[0,Lafter-1].Lbefore和Lafter為災前和災后p波段圖像的灰度級數.對災前遙感圖像進行變換的結果如圖1(c)所示，變換前災前災后遙感圖像中紅色分量的直方圖分布如圖2(a)(b)所示，通過對比可以發現，災前圖像中紅色分量像素灰度級整體偏低，與災后圖像存在紅色分量的直方圖分布差異明顯，造成2幅圖像出現顯著色差.通過直方圖相似變換后，災前圖像的紅色分量直方圖分布整體與災后圖像紅色分量直方圖分布趨于相似，如同圖2(c)所示，對比圖1(b)(c)可以發現，利用直方圖變換后獲得的滑坡發生前遙感圖像與滑坡后遙感圖像整體色差基本消除，相同地物目標在不同遙感圖像中顏色或光譜特征基本趨于一致.

1.2 多源遙感影像尺度空間配準

多源遙感平臺由于成像高度、角度和探測器空間分辨率等因素各不相同，其獲得的圖像空間尺度上也存在一定差異，需要開展災前災后多源遙感影像尺度空間配準.SIFT(scale invariant feature transform)算法在特征點提取方面具有縮放不變性、旋轉不變性、仿射不變性，還有一定的抗光照變化和抗視點變換性能[12-13]，因此本文利用SIFT算法對災前災后滑坡區域遙感圖像進行尺度空間配準，其基本步驟包括：尺度空間極值點提取、提取特征點量化表征和特征點匹配.

對于二維圖像I(x,y)的高斯尺度空間定義為L(x,y,t)，該尺度空間可由參數為t的高斯核G(x,y,t)與I(x,y)卷積得到：

(4)

其中

(5)

其中，x,y,t分別為L(x,y,t)的像素位置坐標和尺度參數.尺度參數t為尺度空間因子，它是高斯正態分布的方差，表征圖像被平滑的程度.為高效地在尺度空間內檢測出穩定的特征點，使用尺度空間中差分高斯(difference of Gaussian, DOG)極值作為判斷依據.DOG算子定義為：

D(x,y,t)=(G(x,y,kt)-G(x,y,t))*
I(x,y)=L(x,y,kt)-L(x,y,t).

(6)

利用多級金字塔模型檢測DOG局部極值.將每個像素跟同一尺度的周圍鄰域8個像素和相鄰尺度對應位置的9×2個像素總共26個像素點進行比較，尋找極值點并作為特征點保存.根據鄰域像素梯度分布特性，以梯度模m(x,y)和梯度方向θ(x,y)對高斯圖像中每個點L(x,y)進行量化表征.

m(x,y)={[L(x+1,y)-L(x-1,y)]2+
[L(x,y+1)-L(x,y-1)]2},

(7)

(8)

利用SIFT算法分別提取災前災后遙感影像對應特征點的m(x,y)與θ(x,y)，利用二維范數尋找最優匹配特征點，實現滑坡災前災后遙感影像的特征點匹配，圖1(b)(c)的特征點匹配結果如圖3所示：

Fig. 3 Registration result of characteristic points for remote sensing images before and after the landslide occurrence圖3 滑坡災前災后遙感影像特征點匹配結果

1.3 基于時空譜特征融合的遙感影像數據集構建

以經過光譜空間和尺度空間配準的災前災后遙感圖像構建新的多維特征融合影像數據集.將滑坡區域災前災后時序變化特征轉化為多光譜圖像的光譜特征，為后期時空譜特征融合應用奠定基礎.以Sentinel-2遙感圖像為例，在綜合分析滑坡地區背景地物光學特征、遙感成像有效波段和圖像分辨率的基礎上，選擇能覆蓋滑坡顏色和光譜特征的B(440～538 nm),G(537～582 nm),R(646～684 nm)和NR(760～908 nm)這4個典型特征波段災前災后圖像構建新的影像數據，如圖4所示.新構建的影像數據中，前4個波段為目標區域災前遙感圖像，后4個波段為目標區域災后遙感圖像.由于經過遙感光譜空間和尺度空間配準，圖像中相同位置地物目標光譜信息和災前災后圖像變化信息均整合為新構建的8波段影像數據的光譜特征信息.在新構建的數據中部分地物目標的光譜分布曲線如圖5所示.從圖5中可以看出，植物和建筑物目標災前災后圖像中未發生明顯變化，因此災前和災后B,G,R,NR這4個波段的光譜曲線空間分布基本相似，而新翻土地在災前災后圖像中發生了變化，災前和災后4個波段的光譜曲線空間分布出現差異，這種差異主要為目標區域的變化信息，通過對該信息的有效利用即可實現對特定地物目標的變化檢測和識別.

Fig. 4 The data of new eight bands images圖4 新構建的8波段遙感影像數據集

Fig. 5 The object spectrum of new eight bands images set圖5 新構建的8波段影像數據中部分地物光譜分布曲線

2 基于時空譜特征融合的滑坡檢測方法

2.1 基于SVM的滑坡地物目標檢測

SVM是一種基于VC(vapnick-chervonenkis)維理論和結構風險最小化原理的經典機器學習分類方法，在圖像處理和目標識別等方面都有著重要應用.本文將其用于新構建8波段遙感圖像中滑坡目標檢測.SVM求解最優超平面問題等價于求解如下方程[14]：

(9)

(10)

其中K(·)是滿足Mercer條件的核函數.

2.2 滑坡典型空間形狀特征分選

發生滑坡區域山體表面的植被會受到不同程度的破壞，新裸露出的土石在亮度、色調、空間形狀方面與周圍植被和人工地物有明顯差異.本文選取滑坡中較為常見的2種基本形狀——簸箕形和三角形，構建基礎形狀模型，如圖6所示：

Fig. 6 Fundamental shape models of landslides圖6 滑坡常見基礎形狀模型

通過軸向長寬比、面積參數和不變矩特征對遙感圖像滑坡區域進行分選，實現對滑坡目標的精確識別.

1) 軸向長寬比

通過對滑坡基礎形狀模型的分析，建立表征滑坡形狀特征的軸向長寬比(如圖7所示)：

(11)

Fig. 7 The axis length ratio of fundamental shape models of landslides圖7 滑坡常見基礎形狀模型的軸向長寬比

2) 面積參數

當大面積區域內發生滑坡災害時，小型滑坡數量較多，但通常造成危害不大；而特大型滑坡或者巨型滑坡數量通常較少，前期可能已有一定的監測數據.這2類滑坡在遙感檢測與識別中可適度弱化，而對于中、大型滑坡其危害性和未知性均較大，將成為遙感檢測與識別的重點.因此本文通過設置合理的滑坡區域面積參數Rsize，實現對典型滑坡的分選與識別.

3) 不變矩特征

考慮到滑坡在遙感圖像中可能存在旋轉和尺度的不確定性，本文利用不變矩特征作為檢測依據，對滑坡遙感檢測結果進行特征匹配識別.對于圖像I(x,y)的階幾何矩(標準矩)定義為[11]：

(12)

其中，p=0,1,…，W;q=0,1,…，H;W和H為圖像的寬高.p+q階中心矩定義為

(13)

(14)

利用二階和三階歸一化中心矩構造7個不變矩特征參數M1～M7如下：

(15)

提取滑坡檢測結果中的符合面積參數要求的連通區域作為疑似滑坡區域，分別計算基礎形狀模型與各疑似滑坡區域的不變矩特征參數，以不變矩的絕對差值作為疑似滑坡區域的選擇標準.

(16)

其中，Mi-possible為疑似滑坡區域不變矩特征參數M1～M7，i∈{1,2,…,7}，Mi-standard為基礎形狀模型不變矩特征參數M1～M7，i∈{1,2,…,7}.

2.3 多源遙感圖像時空譜特征融合檢測流程

本文提出的基于多源遙感時空譜特征融合的滑坡災害檢測方法以滑坡區域災前災后多光譜遙感圖像為基礎，經過多源遙感影像圖像光譜空間配準和尺度空間配準后，構建包含地物光譜和時變特征信息的多特征融合遙感影像數據集，在此基礎上利用SVM對滑坡區域進行檢測識別，獲得疑似滑坡區域，再結合滑坡空間形狀基礎特征模型，依據軸向長寬比、面積參數和不變矩等特征指標對滑坡區域進行精確分選識別，實現對大面積范圍內的滑坡災害快速檢測和精確識別.本文提出方法的基本流程如圖8所示:

Fig. 8 Basic flow chart of the proposed method in our paper圖8 本文滑坡遙感檢測方法流程圖

Fig. 9 Field investigation and information collection by UAV圖9 現地勘察和無人機信息采集

3 實驗結果與分析

為驗證本文提出方法的有效性，利用Sentinel-2獲得的重慶市沙坪壩大學城某地區的滑坡前后多光譜遙感圖像(B,G,R,NR這4個波段)進行了滑坡遙感檢測實驗.實驗將目標檢測和識別中常用的大津法(Otsu)、Bayes分類器、傳統的SVM、傳統的SVDD(support vector data description)、神經網絡(neural networks, NN)等方法[15-17]與本文方法的檢測識別效果進行了對比.算法運行平臺：Intel Core i5-4210 2.4 GHz CPU， 4 GB內存，MATLAB2010.通過現地勘察和無人機信息采集后確定訓練樣本和檢驗樣本區域(如圖9和圖10所示)，在遙感圖上標記訓練樣本分布如圖10所示，訓練用滑坡樣本25個、訓練用背景樣本112個，其中，植被樣本32個、建筑物樣本32個、水體樣本16個、訓練場樣本32個.檢測效果量化評估中所用到的檢測樣本如圖11所示，其中，滑坡樣本137個、背景地物樣本514個.各對比方法參數設置和使用的波段如表1所示.各種檢測方法對滑坡及背景地物的識別情況如圖12所示.對滑坡目標識別精度如表2所示.

Fig.10 Distribution of training samples for each detection method圖10 各檢測方法在所用到訓練樣本分布

Fig. 11 Distribution of test samples圖11 檢測樣本分布情況

Table 1 Used Bands and Parameter of Comparative Methods表1 各對比方法使用的波段和參數設置

Fig. 12 The detection result of comparative methods on Sentinel-2 remote sensing images圖12 Sentinel-2遙感圖像中各對比方法滑坡識別結果

Table 2 Accuracy of Each Method for Landslide Detection表2 各種方法滑坡目標識別精度表 %

結合圖12和表2可以看出，傳統的Otsu方法主要根據地物目標的亮度信息對目標進行分割和識別，難以區分建筑物或道路等亮度較高的人工地物目標，整體識別精度較差，虛警和漏警的像素均較多.Bayes分類器對滑坡識別精度較高，但虛警率也較高，檢測結果圖出現了大量的噪聲像素點.利用災后的R,G,B,NR這4個波段進行基于SVM的滑坡檢測，其正確識別率達67.88%，但漏警率較高，并且與滑坡顏色相近的許多訓練場也存在識別錯誤.SVDD方法整體識別效果較好，但虛警率較高，結合圖12(d)可以看出，該方法對建筑和道路等目標仍難以區分.NN方法能區分大部分建筑和道路目標，但對滑坡的正確識別率和漏警率均不理想.利用災前災后8波段圖像進行基于SVM的滑坡檢測方法整體識別效果較好，結合圖12(f)可以看出，部分新建工程開挖的地面由于其光譜特征和變化特征與滑坡較為相似，因此被錯誤地識別成了滑坡目標，所以其誤警率達10.70%.本文方法在基于SVM的檢測方法基礎上加入滑坡基礎形狀特征，實現了對檢測結果的精確分選和識別，最終的識別效果如圖12(g)所示，結合表2可知，本文方法有效去除了誤分類點，滑坡的正確識別率達到95%以上，同時誤警率下降到了4.47%.通過多種滑坡檢測和識別方法的對比可以看出，本文提出的方法整體識別精確優于其他多種方法.

4 檢測方法討論

本文提出的滑坡災害檢測與識別方法主要利用了遙感圖像中目標的光譜特征、變化特征和空間形狀特征.具體實現過程如下：

1) 在經光譜-尺度空間配準的災前和災后圖像基礎上構建多維特征融合的遙感影像數據集，將災前災后遙感圖像中光譜信息和變化信息融合后表征為新建遙感影像數據集的各種光譜曲線.

2) 利用SVM對各種模式的光譜曲線進行識別，提取疑似滑坡區域.此時SVM識別滑坡目標的主要依據為目標的光譜信息和變化信息，結合表2 中SVM(8 Bands)和SVM(4 Bands)兩欄數據可以看出，融合變換信息和光譜信息后，滑坡正確識別率提升了27.01%.

3) 結合滑坡典型空間形狀特征對SVM的識別結果進行分選和優化.該環節將滑坡的空間形狀信息融合到目標識別算法中，最終實現光譜特征、變化特征和形狀特征的融合檢測應用.結合表2中實驗數據可知，本文最終方法的滑坡正確識別率得到進一步提升，同時誤警率較未加形狀特征的SVM(8 Bands)下降了6.23%.

在滑坡災害遙感探測和災情提取方面，災前災后遙感圖像對比分析提取常常是比較有效的方法，但受衛星重返周期和氣象條件的影響，近實時的災后圖像和滑坡前短期內災前圖像往往不容易獲取，災前圖像和災后圖像較大的時間差可能引入較多的變化信息，影響本文方法對滑坡識別的精度.

5 結論與展望

本文從遙感圖像時序變化特征、光譜特征和空間形狀特征融合利用角度出發，提出一種基于多源遙感時空譜特征融合的滑坡災害檢測方法，在多源遙感圖像光譜空間配準和尺度空間配準的基礎上構建基于多維特征融合的遙感影像數據集，實現目標區域時序變化和光譜特征的有效融合，利用SVM對滑坡進行檢測與識別，結合空間形狀特征對檢測結果進行分選和優化，實現光學遙感影像時、空、譜多維特征信息的滑坡檢測應用，為滑坡災害的遙感檢測和識別提供了一種新方法.

下一步工作中我們將研究航空、航天多源遙感信息融合處理和變化特征提取方法，通過多源遙感信息融合和互補提升變化特征提取的有效性和準確性，進而實現對滑坡目標的快速精確探測.