TOPSIS融合極化信息的復雜背景下滑坡檢測方法

高歐陽， 牛朝陽， 劉 偉， 張浩波， 鄒瑋琦， 張雅歌

(中國人民解放軍戰(zhàn)略支援部隊信息工程大學數(shù)據(jù)與目標工程學院, 河南鄭州 450001)

0 引 言

利用光學或合成孔徑雷達(Synthetic Aperture Radar, SAR)遙感衛(wèi)星影像檢測災后滑坡的位置和邊界，調(diào)查滑坡的空間分布特征，能夠為災害防治與災區(qū)重建提供重要信息[1]。滑坡的形成常常伴隨有惡劣的天氣情況，其中降雨是滑坡最主要的誘發(fā)因素之一[2]，陰雨天氣云層的遮擋會導致光學遙感衛(wèi)星探測手段的失效。SAR作為一種主動式對地觀測系統(tǒng)，通常工作在微波波段，波長較長，基本不受太陽光照、氣候條件等環(huán)境因素的影響，能夠進行全天時、全天候的不間斷觀測[3]。利用極化合成孔徑雷達(Polarimetric Synthetic Aperture Radar，PolSAR)圖像進行滑坡檢測已經(jīng)成為重要的研究方向。

PolSAR通過發(fā)射和接收極化狀態(tài)正交的電磁波，能夠獲取地物的全極化信息[4]，基于極化分解技術(shù)對地物散射機理進行建模與分析，能夠準確地理解地物的物理特性[5]。在極化分解的基礎上進行滑坡檢測主要有兩種思路[6]，一是利用多閾值法進行滑坡檢測，即對多項極化特征分別設置閾值，同時滿足所有閾值條件的區(qū)域被確定為滑坡區(qū)域，此類方法具有簡單易行的優(yōu)點，但是在實施僅限于特定環(huán)境并需要人工閾值的干預，普適性不強；二是利用監(jiān)督分類或非監(jiān)督分類的方法對災后PolSAR圖像進行分類以區(qū)分出屬于滑坡的區(qū)域。相比于多閾值法，后者得益于分類器的復雜結(jié)構(gòu)，可以提供更高的滑坡檢測準確率，因此分類法比多閾值法更常用于滑坡檢測。Niu等[6]在總結(jié)前人工作的基礎上，首次將變化檢測(Change Detection，CD)與層次分析法(Analytic Hierarchy Process，AHP)相結(jié)合實現(xiàn)復雜地貌背景下的滑坡檢測。該方法用CD檢測發(fā)生變化的區(qū)域，即預篩選滑坡區(qū)域，然后通過AHP剔除其他變化區(qū)域?qū)聶z測造成的干擾，但是從復雜地貌背景中檢測滑坡的準確率為73.25%，田地對滑坡檢測造成的虛警率為14.74%，還有較大的改善空間。

綜上，極化測量技術(shù)可以充分利用極化信息來區(qū)分不同散射特性的地貌，但是對于散射特性相似的地貌難以取得理想的效果。因此，本文在前人研究的基礎上，首先通過優(yōu)劣解距離法(Technique for Order Preference by Similarity to an Ideal Solution，TOPSIS)綜合多個極化特征參數(shù)，充分利用PolSAR圖像所包含的極化信息提取出與滑坡散射特征相似的區(qū)域，實現(xiàn)對疑似滑坡區(qū)域的初步篩選；然后通過相干圖對地形地貌變化非常敏感的特性來確定發(fā)生變化的區(qū)域，從而抑制疑似滑坡區(qū)域中未變化區(qū)域造成的虛警，實現(xiàn)從復雜地貌背景中有效檢測滑坡的目的。

1 研究方法

本文方法的處理流程如圖1所示，包含3個階段。第一階段是預處理過程，對災后PolSAR數(shù)據(jù)進行輻射定標、多視和地形校正。第二階段為滑坡檢測過程，首先通過Yamaguchi分解、H/A/α分解和相關系數(shù)計算得到表面散射ps、二次散射pd、體散射pv、螺旋體散射ph[7]、極化熵H、平均散射角α、反熵A[8]和同極化分量相關系數(shù)的實部Re(ρhh-vv)8個極化特征參數(shù)，然后利用帶權(quán)重的TOPSIS法充分融合上述極化特征參數(shù)信息，經(jīng)閾值分割和形態(tài)學處理得到與滑坡極化散射特性相同的區(qū)域(下文均稱為“疑似滑坡區(qū)域”)。第三階段為虛警抑制過程。利用災前單極化SAR數(shù)據(jù)與災后PolSAR數(shù)據(jù)的相同極化組合圖像構(gòu)建雙時相數(shù)據(jù)，經(jīng)干涉處理生成相干圖以提取出變化區(qū)域，然后將通過邏輯“與”運算將兩個結(jié)果進行融合，從而去除與滑坡表面散射行為一致的未變化區(qū)域造成的大量檢測虛警。

圖1 本文方法實現(xiàn)滑坡檢測的流程圖

本文方法中的形態(tài)學運算為先“開”后“閉”運算，即首先通過“開”運算去除二值化圖像黑色背景中的白色斑點(如椒鹽噪聲)，然后通過“閉”運算來消除白色區(qū)域內(nèi)的黑色空洞、連接鄰近的白色區(qū)域并平滑邊界[9]，從而使檢測結(jié)果圖的區(qū)域連通性更好，邊界特征更加明顯。下文重點闡述極化特征參數(shù)提取、TOPSIS檢測與相干圖計算等主要處理步驟。

1.1 極化特征參數(shù)提取

1.1.1 Yamaguchi分解

Yamaguchi分解是一種基于散射模型的非相干分解算法，可以很好地匹配地物的基本散射機制。該分解方法從極化相干矩陣T中提取出各個散射分量對應的功率大小，即表面散射功率Ps、二次散射功率Pd、體散射功率Pv和螺旋體散射功率Ph。最后，極化相干矩陣T可以表示為各散射功率的加權(quán)組合形式：

T=fsTs+fdTd+fvTv+fhTh

(1)

式中：Ts，Td，Tv，Th分別為表面散射、二次散射、體散射和螺旋體散射分量對應的極化相干矩陣；fs，fd，fv，fh分別為Ps，Pd，Pv和Ph的分解系數(shù)[7]。

1.1.2H/A/α分解

H/A/α分解是一種利用協(xié)方差矩陣或者極化相干矩陣的特征值和特征向量定義極化特征的非相干分解算法。在單站互易條件下，極化相干矩陣T為3×3的半正定厄密特矩陣[10]，因此可以作如下特征值分解：

(2)

式中，特征值λi(1≤i≤3)滿足λ1≥λ2≥λ3>0，λi對應的特征向量ui相互正交，H表示共軛轉(zhuǎn)置。則極化熵H、平均散射角α和反熵A可由特征值λi和特征向量ui計算得到：

(3)

1.1.3 相關系數(shù)計算

同極化分量相關系數(shù)ρhh-vv為目標的兩個同極化散射矢量之間的相關系數(shù)，取值越大表明目標的散射類型越單一，例如Re(ρhh-vv)的正值部分常用于反映地物的表面散射特性[11]。通過極化相干矩陣T計算Re(ρhh-vv)的公式如下：

Re(ρhh-vv)=

(4)

式中，Re(·)表示取實部，T11，T12和T22表示極化相干矩陣T中的相應元素。

1.2 TOPSIS提取疑似滑坡區(qū)域

1.2.1 TOPSIS用于PolSAR圖像滑坡檢測的基本原理

TOPSIS屬于多屬性決策分析中的常用方法，也可以稱為逼近理想解法，其對數(shù)據(jù)分布及樣本含量沒有嚴格限制，能夠充分利用原始數(shù)據(jù)信息，并且評價結(jié)果比較精確可靠。該方法根據(jù)一系列屬性條件(即評價指標)對有限個評價對象與理想解的接近程度進行排序[12]。假設采用M個評價指標對N個評價對象分別評價，其中第i個(i=1,2,3,…,N)評價對象在第j個(j=1,2,3,…,M)評價指標上的數(shù)值為zNM，可以得到評估矩陣如下：

(5)

進一步定義評估矩陣的正、負理想解分別為

(6)

(7)

第i個評價對象與正、負理想解的距離分別為

(8)

(9)

式中wi表示每個評價指標的權(quán)重，可以根據(jù)不同的需求通過不同的方法來確定。

最后計算評價對象的最終得分，即其與正理想解的相對貼近度為

(10)

相對貼近度Si的大小反映了評價對象與正理想解的差距，Si越大說明評價對象越接近正理想解，Si越小說明評價對象越遠離正理想解。

利用TOPSIS將8個極化特征參數(shù)(ps，pd，pv，ph，H，A，α，Re(ρhh-vv))有機融合，可以在充分利用極化信息的基礎上精確計算PolSAR圖像每個像素與滑坡的相似度，即相對貼近度。如圖2所示，將每一個像素點作為一個評價對象，極化特征參數(shù)作為評價對象得分的一個評價指標，然后利用TOPSIS計算出每個像素與理想解的相對貼近度Si，然后選擇合適的閾值對Si進行二值化處理即可得到疑似滑坡區(qū)域。其中，為了確保TOPSIS方法對不同PolSAR圖像數(shù)據(jù)的適用性，確定正理想解為1，負理想解為0。

圖2 TOPSIS提取疑似滑坡區(qū)域的關鍵步驟

1.2.2 極化特征參數(shù)的歸一化和同趨勢化處理

因為8個極化特征參數(shù)(ps，pd，pv，ph，H，A，α，Re(ρhh-vv))作為評價指標的指向性、計量單位、數(shù)值評價方式和評價值范圍存在的差異會對評價結(jié)果產(chǎn)生不利影響，所以需要在構(gòu)建評估矩陣之前對這些極化特征參數(shù)進行歸一化和同趨勢化處理，使各項參數(shù)數(shù)值都在區(qū)間[0,1]內(nèi)，并且與滑坡散射特性成正相關關系，即數(shù)值越大，越能體現(xiàn)滑坡的散射特性。

1) 對Yamaguchi四分量進行歸一化處理[6]。

px=Px/(Ps+Pd+Pv+Ph)

(11)

式中，下標x∈{s,d,v,h}，s，d，v和h分別表示表面散射、二次散射、體散射和螺旋體散射。由于滑坡發(fā)生后地表覆蓋的裸露土壤[11]通常以表面散射的形式所呈現(xiàn)，與其他地物有很大的區(qū)別[6]，所以ps值越大、pd,pv,ph的值越小，越能體現(xiàn)滑坡地貌的表面散射類型，因此本文對歸一化處理后的pd,pv,ph進行正向化處理。

px=1-px

(12)

2) 對于極化熵H、平均散射角α，其中當H屬于0.52～0.63且α屬于29～37的取值范圍時更能體現(xiàn)滑坡的散射特性[13]，所以本文利用區(qū)間型歸一化方式來對H和α的值進行處理。

(13)

rij=1

(14)

3) Re(ρhh-vv)的正值部分常用于識別表面散射地物類型，其值越大表明地物表面散射類型更為顯著，所以將Re(ρhh-vv)的負值部分賦值為0。

1.2.3 AHP方法確定極化特征參數(shù)的權(quán)重

不同的極化特征參數(shù)對于最終評分的影響不同，因此需要對它們賦予不同的權(quán)重。AHP是美國運籌學家Thomas提出的一種分析復雜決策的半定量方法，Niu等[6]將AHP運用到PolSAR圖像滑坡檢測中。在此，本文使用AHP來確定極化特征參數(shù)的權(quán)重。根據(jù)滑坡災后的地貌特征，ps最能體現(xiàn)滑坡的表面散射特性，H和α可以較好地區(qū)分滑坡與其他地貌[13]，Re(ρhh-vv)常用于反應地貌的表面散射特性，pd，pv，ph和A將對識別滑坡作用較小。因此將ps作為最重要的評價指標，H和α作為比較重要的評價指標，Re(ρhh-vv)作為一般重要的評價指標，pd，pv，ph和A作為相對不重要的評價指標，然后根據(jù)文獻[14]中的成對比較矩陣重要性標度表確定不同評價指標之間的重要性比較值，構(gòu)建判斷矩陣如表1所示。為了盡可能消除判斷矩陣中存在的主觀性，需要引入一致性比率(CR)作為參考，當CR<0.1時判斷矩陣被認為是合理的。經(jīng)過計算，表1所示的判斷矩陣CR=0，說明各個評價指標的權(quán)重分配是合理的。

表1 AHP確定權(quán)重的判斷矩陣

1.3 基于相干圖抑制滑坡檢測的虛警

在利用TOPSIS方法從復雜地貌背景中檢測滑坡的結(jié)果中，田地、裸土等其他以表面散射類型為主的地物地貌，可能對滑坡檢測結(jié)果帶來較多的虛警。滑坡作為一種突發(fā)性自然災害，災害區(qū)域具有十分顯著的突變性，而田地、裸土等地物地貌一般不具有這種特征，因此可通過對地物地貌的突變性分析抑制滑坡檢測的虛警。地形地貌的顯著變化通常會在SAR干涉處理中引起嚴重的失相干，即相干性較低，而未發(fā)生變化的區(qū)域?qū)⒈３州^高的相干性，因此可以通過相干圖來判斷某一區(qū)域是否發(fā)生變化，進而抑制未變化地物造成的虛警，使滑坡檢測結(jié)果更加準確。

利用雙時相SAR圖像可以得到相干圖，相干圖中灰度值的大小表示相干性的高低，相干性|γ|計算如下：

(15)

其中，E[ ]表示數(shù)學期望，S1，S2分別為滑坡前后的SAR影像(極化組合相同)，相干性|γ|的取值范圍為[0,1]，|γ|越大說明一個區(qū)域的變化程度越低，|γ|越小說明該區(qū)域的變化越明顯[15]。

在不同極化組合的SAR圖像中，垂直極化后向散射更能體現(xiàn)以水平結(jié)構(gòu)為主的地貌，因此本文采用能突出裸土、裸巖區(qū)域特征[16]的VV極化組合圖像生成研究區(qū)域的相干圖，然后對相干圖進行閾值分割和形態(tài)學處理得出變化區(qū)域。

2 研究區(qū)域與數(shù)據(jù)源

研究區(qū)域包括兩個。區(qū)域A：2015年12月20日由于暴雨引起的山體滑坡發(fā)生于深圳恒泰裕工業(yè)園區(qū)，本次滑坡發(fā)生在復雜的地貌背景中(包含滑坡、水系、森林、居民地、田地五類基本地物)，并且存在其他自然或人類活動引起的變化，因此使用該區(qū)域作為研究對象具有一定的代表性。區(qū)域B：2018年9月6日北海道地震引起的山體滑坡，滑坡主要發(fā)生在植被覆蓋的山區(qū)，其周邊也包含滑坡、水系、森林、居民地、田地五類基本地物類型，但是由于滑坡區(qū)域和其他地物距離較遠，在同一幅圖上展示會因比例尺過小導致目標區(qū)域難以判讀，為了分析不同地貌對滑坡區(qū)域檢測的影響，本文從同一景數(shù)據(jù)中截取不同地貌區(qū)域合并為同一幅圖像。兩個研究區(qū)域均使用日本ALOS-2衛(wèi)星的全極化數(shù)據(jù)，距離向分辨率為5.13 m，方位向分辨率為3.21 m，具體參數(shù)如表2所示。

表2 衛(wèi)星數(shù)據(jù)基本參數(shù)

兩個研究區(qū)域的Yamaguchi分解假彩色圖像和光學圖像如圖3所示。其中假彩色圖像的紅色代表二次散射，綠色代表體散射，藍色代表表面散射。受災前后的光學圖像與PolSAR圖像的采集時間存在一定的差別，但是其仍然能夠作為分析討論滑坡引起地形地貌變化的參考。通過目視判讀后在假彩色圖上利用白色框選出滑坡、水系、森林、居民地、田地五類地物樣本，以便于在下文作出檢測正確率分析。

(a) 研究區(qū)域A的假彩色圖像和光學圖像

(b) 研究區(qū)域B的假彩色圖像和光學圖像圖3 研究區(qū)域基本情況圖

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

3.1 滑坡檢測性能的定量分析

為了定量分析本文方法的有效性，通過圖3標記區(qū)域的不同地物像素樣本計算滑坡檢測正確率和總虛警率(即滑坡以外地物樣本區(qū)域內(nèi)被檢測為滑坡的像素總數(shù)除以樣本區(qū)域像素總數(shù))，與文獻[6]的CD+AHP方法作出對比，并在圖4中描繪出不同方法在兩個研究區(qū)域的正確率和總虛警率關系曲線。其中4條曲線的橫坐標虛警率均小于1，這是由于CD+AHP方法檢測的虛警被限制在變化檢測區(qū)域范圍內(nèi)，本文方法檢測的虛警被限制在相干圖分割區(qū)域范圍內(nèi)。因為地形地貌上的差異導致相同方法在不同研究區(qū)域有著不同的性能表現(xiàn)，對于每個研究區(qū)域，本文方法和CD+AHP方法相比，在相等虛警率的情況下滑坡檢測正確率均有大幅提升，表明了其檢測滑坡的有效性。

為了進一步分析不同方法的滑坡檢測性能，選擇合適的閾值計算本文方法和CD+AHP方法在相同正確率下不同地物造成的虛警率，以及相同總虛警率下滑坡檢測正確率和不同地物造成的虛警率，并與文獻[11]的MT方法作出對比，如表3所示。其中，CD+AHP方法的滑坡檢測正確率和虛警率均由文獻[6]建議的閾值得到，MT方法中的多閾值為文獻[11]設定的經(jīng)驗值，該閾值無法調(diào)整得到與另外兩種方法相同的滑坡檢測正確率和總虛警率。

圖4 不同方法的滑坡檢測正確率和總虛警率關系曲線圖

表3 不同方法滑坡檢測正確率和虛警率對比

對于研究區(qū)域A，MT方法的檢測結(jié)果給出了最高的滑坡檢測正確率，但是田地和水系等以表面散射為主的地物造成了極大的虛警率，使得該方法幾乎失效；CD+AHP方法則以檢測正確率小幅度降低為代價，較好地抑制了各類地物的虛警率，尤其是田地和水系的虛警率大幅度降低；與CD+AHP方法相比，在相同的正確率(73.46%左右)下，本文方法將總虛警率和各類地物的虛警率控制在極低的水平(均低于1%)；在相同的總虛警率(3.45%左右)下，本文方法的滑坡檢測正確率提高了約11%，并且較好地抑制了大部分地物的虛警率，只有田地造成的虛警率較高，這是由于田地區(qū)域的農(nóng)作物非常容易受到自然生長或者生產(chǎn)活動的影響而發(fā)生變化，因此通過相干圖并不能完全剔除田地區(qū)域帶來的虛警。

對于研究區(qū)域B，由于滑坡區(qū)域內(nèi)大量存在以體散射和二次散射為主的殘破樹木、樹干的干擾，導致3個方法的滑坡檢測正確率均有明顯降低。其中，MT方法檢測正確率僅為55.45%，表明該方法的檢測效果較差；CD+AHP方法提供了差強人意的68.61%檢測正確率，同時對虛警的抑制也有明顯提升；在相同的正確率(68.61%左右)下，本文方法與CD+AHP方法相比，總虛警率和各類地物虛警率都降低至理想水平，在與CD+AHP方法相同的總虛警率(8.18%左右)下，本文方法的檢測正確率提高了約16%，達到了較為合理的數(shù)值(84.48%)。

3.2 滑坡檢測結(jié)果圖的分析討論

為了更加直觀地分析不同方法滑坡檢測的效果，將表3對應的結(jié)果圖在圖5中給出，其中本文方法的結(jié)果圖對應表3中上標②所在行的數(shù)值。

(a) 研究區(qū)域A不同方法的滑坡檢測結(jié)果圖

(b) 研究區(qū)域B不同方法的滑坡檢測結(jié)果圖圖5 不同方法滑坡檢測結(jié)果圖對比

通過對比Yamaguchi假彩色圖像和光學圖像目視解譯得到滑坡區(qū)域范圍(在圖5中用黃色標記)。對于研究區(qū)域A，三種方法均能將滑坡與二次散射為主的居民地和體散射為主的森林有效地區(qū)分開，同時在滑坡區(qū)域存在散落分布的小面積黑色空洞，這是由于救災活動中該區(qū)域存在大量的大型機械(比如挖掘機、鏟車、運土車等)，機械臂、車體與地面之間形成二次散射結(jié)構(gòu)，從而對滑坡檢測造成了一定的干擾，但是并不影響大面積滑坡范圍的確定。MT方法檢測結(jié)果中，田地和水系等以表面散射為主的地物造成了大量的虛警，對滑坡區(qū)域的提取造成了很大影響；CD+AHP方法檢測結(jié)果中，森林、居民地等地物類型造成的虛警均較少，但是明顯將少部分田地和水系錯誤地檢測為滑坡，同時滑坡頭部(區(qū)域1)存在明顯空洞，這是由于滑坡前后該區(qū)域均為裸土，主要散射類型沒有發(fā)生變化所致。本文方法得出的滑坡區(qū)域更加完整(包括區(qū)域1)且連通性更強，而且背景中大部分地物的虛警得到了更好的抑制，但是本文方法在田地(區(qū)域2)的虛警稍多于CD+AHP方法，這是由于農(nóng)作物的生長導致了該區(qū)域相干性較低，通過相干圖不能完全抑制該區(qū)域的虛警；區(qū)域3的森林覆蓋率較低，中間夾雜了少量裸土區(qū)域，因為這些區(qū)域在滑坡前后未發(fā)生明顯變化，因此MT方法僅通過滑坡后單時相PolSAR數(shù)據(jù)錯誤地將此類區(qū)域錯判為滑坡，而CD+AHP方法和本文方法的虛警抑制效果較好，這是因為兩種方法采用不同的方式(CD和相干圖)來檢測研究區(qū)域中的變化，進而剔除未變化區(qū)域的虛警。

對于研究區(qū)域B，MT方法和CD+AHP方法均難以在提取出完整的滑坡區(qū)域信息，本文方法則可以在復雜背景中有效地抑制虛警的同時提供較好的滑坡檢測結(jié)果，兩個研究區(qū)域檢測結(jié)果的一致性證明了本文方法可以實現(xiàn)復雜背景下檢測滑坡的目的。

4 結(jié)束語

針對現(xiàn)有方法在復雜背景下滑坡檢測正確率低、虛警率高的問題，本文提出了一種利用TOPSIS融合極化信息的滑坡檢測方法，通過綜合多個極化特征參數(shù)，以不同地貌的散射特性為基礎來檢測滑坡，進一步使用相干圖剔除虛警，可以得到清晰的滑坡區(qū)域范圍，大幅度提高滑坡檢測正確率，而且可以有效抑制大部分地貌造成的虛警。隨著PolSAR技術(shù)的不斷發(fā)展，利用PolSAR數(shù)據(jù)進行災害檢測將會發(fā)揮越來越重要的作用，本文方法在不依賴先驗信息的條件下可以準確快速地提取出滑坡發(fā)生后災區(qū)的信息，為復雜背景下的滑坡檢測提供了具體可行的方案，為災后的救援、二次災害的預警和災后重建等活動提供可靠的依據(jù)，為今后相關方面的研究提供了一定的參考。