聯合SBAS-InSAR 和PSO-BP 算法的高山峽谷區地質災害危險性評價

周定義，左小清※，喜文飛,2，肖 波,3，游 洪

（1. 昆明理工大學國土資源工程學院，昆明 650093； 2. 云南師范大學地理學部，昆明 650500；3. 云南交通職業技術學院公路學院，昆明 650500）

0 引 言

中國是地質災害高發的國家，各類地質災害給人民的生命和財產造成了巨大的損失[1]。滑坡、泥石流、崩塌作為主要的地質災害，具有隱蔽性強、危害性大、突發性高等特點，廣泛分布于中國山區和峽谷地帶[2-4]。近年來，高山峽谷區地質災害頻頻發生，給人民生命和財產安全帶來嚴重威脅。因此，對高山峽谷地區地質災害開展危險性評價研究，能夠為防災減災事業及政府部門決策提供參考。

目前，國內外對于地質災害危險性評價主要分為三類：經驗模型[5-6]、統計分析模型[7-9]、機器學習模型[10-12]。以上三類方法證明了其在地質災害評價中的有效性。但存在以下兩個弊端：一是時效性差、識別不準確：基于收集的歷史災害數據，以前是滑坡或泥石流等地質災害區域，可能已經被治理好的情況，如果再將其定義為地質災害點，容易導致評價結果不準確。運用RS（Remote Sensing）技術解譯得到的地質災害點范圍，無法識別潛在地質災害，特別是在高山峽谷區。二是過多權值計算：三類評價模型都需要運用權重和信息量來進行地質災害統計，這無疑加大專家參與的計算量。為此，亟需提出一種可行有效方法解決此弊端。

在以往研究中，筆者曾提出一種基于小基線集合成孔徑雷達干涉（Small Baseline Subsets Interferometric Synthetic Aperture Radar，SBAS-InSAR）和粒子群優化-反向傳播（Particle Swarm Optimization-Back Propagation，PSO-BP）神經網絡算法的礦區地表沉降監測及預測模型[13]，最終證明了SBAS-InSAR 在礦區地表長時間沉降監測中的優勢以及 PSO-BP 模型在礦區地表沉降預測中的有效性和合理性。近年來，眾多學者運用小基線集合成孔徑雷達干涉在地質災害識別方面取得大量成功范例[14-15]，運用SBAS-InSAR 技術在地質災害識別探測中的應用潛力已經得到國內外眾多學者的認可。誤差反向傳播（Back Propagation，BP）學習算法作為人工神經網絡（Artificial Neural Network，ANN）最為廣泛運用的一種神經網絡，已經成功應用于各類地質災害危險性評價及預測等方面并證明了其有效性[16-18]。與傳統的基于梯度下降的優化算法相比，粒子群（Particle Swarm Optimization，PSO）算法具有魯棒性強、擴充性好、不易陷入局部最優等優點；而與進化規劃、遺傳算法等基于自然界進化過程的優化方法相比，粒子群算法的信息共享機制使得種群可以更快地收斂到最優值[19]。綜上所述，利用SBAS-InSAR和PSO-BP 算法為地質災害危險性評價提供了理論和實踐基礎。

為此，針對現有地質災害危險性評價方法存在的弊端，提出聯合SBAS-InSAR 和PSO-BP 相結合的方式對高山峽谷區地質災害進行危險性評價。利用SBAS-InSAR 獲取現有地質災害點和潛在災害的形變量，根據形變量值和引入高分辨率光學影像等作為輔助識別；利用PSO-BP 算法構建模型，通過學習、訓練、測試、驗證得到危險性指標，以期在減少專家參與進行大量統計計算的前提下，解決地質災害評價時效性和準確性差的弊端。

1 研究方法與數據

1.1 研究區概況

東川小江流域為世界典型暴雨泥石流區，被稱為“泥石流的天然博物館”，該地區地質災害類型主要以滑坡、崩塌、泥石流為主。泥石流的發生常伴隨滑坡，屬于暴雨型滑坡泥石流，即先滑坡，再經暴雨沖刷形成泥石流。東川區地質災害活動由來已久，至今呈加劇發展趨勢，截至2019 年4 月，東川區各類地質災害隱患點數為269，潛在經濟損失290 961 萬元。本文選取小江沿線兩側高山峽谷作為研究區，如圖1 所示。以小江為界，河谷凹陷，形成“V”字型，屬于典型的深切割高山峽谷，東側為牯牛寨山，最高峰海拔4 017.3 m；西部為拱王山，最高峰海拔4 344.1 m。該區域地勢陡峭，獨特的地形和地質構造致使局部區域暴雨多，土壤疏松，水土流失嚴重，進而導致滑坡、崩塌、泥石流等地質災害頻發。

1.2 研究數據

監測數據為2018 年至2019 年的40 期Sentinel-1A 雷達數據，每期1 景，用于監測研究區形變信息。Sentinel-1A的軌道類型為近圓形的太陽同步軌道，其軌道高度為693 km，傾角為98.18°。衛星重訪周期為12 d，選用干涉寬模式（Interferometric Wide，IW）下的S1 TOPS-modes單視復數SLC（Single Look Complex）數據，幅寬250 km，升降軌數據的入射角分別為34.17°和39.35°。該數據為C 波段的升降軌數據，SLC 數據處理后的升降軌數據距離向分辨率為5 m，方位向分辨率為20 m。

輔助數據為谷歌衛星影像數據、坡度、坡向和曲率數據用于主被動遙感數據協同監測識別地質災害；精密定軌星厲數據POD（Precise Orbit Ephemerides）用于提高雷達數據軌道精度；DEM（Digital Elevation Model）數據用于去除地形相位影響；降雨、斷層、河流水系、地層巖性、土地利用和歸一化植被指數（Normalized Difference Vegetation Index，NDVI）作為PSO-BP 模型的輸入層評價因子。以雷達數據序列為基準數據集，數據源見表1。

表1 數據源一覽表Table 1 List of data sources

1.3 研究方法

本研究的總體技術方案如圖2 所示，主要技術流程為：1）地表形變信息獲取及地質災害識別；2）格網單元劃分及評價因子選取；3）模型構建及危險指數評定。

1.3.1 地表形變信息獲取及地質災害識別

首先，依據歐空局下載的40 景升降軌道SCL 影像，利用精密定軌星厲數據POD 修正軌道信息，減小定時及軌道誤差對干涉相位的影響。再次，選擇同時位于時間基線中部及時序影像多普勒質心頻率中心時段的數據作為主影像，對主影像和輔影像兩兩進行粗配準，再用增強譜分集算法進行精配準，通過設置臨界基線和時間基線，生成像對，設置多視數Range Looks 和Azimuth Looks為4∶1，這樣設置可以較好地抑制斑點噪聲，采用Minimum Cost Flow 解纏方法和Goldstein 濾波方法做干涉工作流，最終生成干涉圖，調整刪除不理想的數據，對多視處理后的主從影像進行共軛相乘獲得干涉圖，利用DEM 對和精密定軌數據對干涉相位進行差分處理移除平地效應及地形效應，生成時間序列的干涉相位。具體包括感興趣區選擇、干涉處理、地形相位去除、差分模型構建4 個步驟。主輔影像干涉相位可表示為公式（1）[13]

式中φtopo為地形相位；φdef為形變相位；φatm為大氣延遲相位；φflat為平地相位；φnoise為噪聲相位。通過相位解纏提取有效的形變信息，運用最小二乘（Least Squares，LS）或奇異值分解（Singular Value Decomposition，SVD）算法得到變形速度，最后，對序列信息進行地理編碼獲得研究區域內升降軌下的雷達視線方向（Line of Sight，LOS）的形變信息。

利用不同軌道數據可以互補，使得監測結果更為準確全面，能夠避免單一軌道帶來的幾何畸變等問題。研究區內植被覆蓋度高，失相干嚴重，導致監測地面形變信息的能力較差，對形變監測精度的影響較大[22-23]，因此，引入歸一化植被指數（Normalized Difference Vegetation Index，NDVI）對研究區植被覆蓋數據進行分析，以剔除植被覆蓋對形變監測精度的影響。同時，僅憑獲取到的研究區形變監測結果，無法有效區分研究區形變監測結果是否為地質災害，為此，引入高分辨率光學影像，即借助谷歌衛星影像對獲取的形變區域進行人工識別輔助識別，結合形變范圍、高程、坡度、植被覆蓋和坡體是否具有地質災害特征等進行識別，避免過度依賴形變結果導致的誤判等問題。為驗證地質災害識別結果的準確性，選擇野外實地勘察驗證。

1.3.2 格網單元劃分及評價因子選取

常用的地質災害評價單元包括網格單元、斜坡單元、地形單元和地域單元4 種類型。對于高山峽谷地區，網格單元和斜坡單元是最好的選擇。但斜坡單元具有不連續性，精細程度不高，劃分的正確度沒有保證。網格單元在評價過程中雖然對于地表起伏度不能很好顯示，但是相較于斜坡單元具有更大優勢，網格單元使整個研究區的評價結果過渡平穩，評價結果更加真實可信[24-25]。為此采用網格單元對研究區進行劃分。網格單元大小劃分根據薛凱凱等[26]采用DEM 試驗進行回歸分析而得到的經驗公式（2）：

式中sG為網格單元大小，s為地形比例尺倒數。

查閱資料可知，用于地質災害危險性評價的評價因子主要有高程、坡度、坡向、曲率、降雨量、地層巖性、土地利用類型、距斷層距離、距河流水系距離[7,8,24-25]。本文在上述10 種評價因子的基礎上引入NDVI 和升降軌數據獲取的形變速率作為評價因子。根據識別結果選取部分地質災害區域和選取部分為非地質災害的區域，通過ArcGIS 多值到點工具，提取出地質災害范圍和非地質災害的區域對應格網單元內的12 個評價因子，通過運用SPSSAU 軟件中的灰色關聯和因子分析，驗證以上12 個因子作為評價因子的有效性和可靠性。

1.3.3 模型構建及危險指數評定

傳統的BP 神經網絡存在學習效率低，收斂速度慢等問題，利用粒子群算法來優化BP 神經網絡的權值和閾值，可以減小訓練誤差，提高訓練速度和精度。將識別出來的災害和非災害點集區域定量化處理，用“1”代表高危險區，“0”代表低危險區，運用ArcGIS 提取對應的評價因子，將其作為輸入層，高低危險區作為輸出層，構建PSO-BP 模型，其中隱含層和輸入層激活函數分別選擇tansig、purelin，訓練函數選擇trainlm。選取高危險區和低危險區共4 756 個點進行模型訓練和驗證，從4 756個點中隨機選取其中的4 200 個點作為訓練集，對隨機抽取的556 個點進行驗證，通過訓練得到最優網絡并保存網絡模型，最終得到高、低危險區對應的輸出范圍。運用保存好的PSO-BP 模型對剩余的38 749 個點進行危險性指數預測，結合上述4 756 個點對應的危險性指數，將其得到的危險性指數導入到對應點屬性表內，利用自然間斷點分級法結合專家參與進行危險性分級，得到研究區危險性評價圖。不論是經驗模型、統計分析模型還是機器學習模型，都是在統計評價因子權重信息的基礎上進行，同時均是基于收集的地質災害數據進行危險性評價，為此，綜合考慮，定量角度采用信息量和綜合賦權法，信息量法對數據樣本要求少，可解決多源數據差異性。地質災害是在多種因素相互作用下而形成的，通過計算各個因素的信息量，并加權疊加各個因素的信息量，建立地質災害易發性評價模型，信息量越大，地質災害易發性越高，其計算方法如下：

式中I為評價區某單元信息量預測值；Ni為分布在因素Xi內特定類別內的災害點單元數；N為研究區含有災害點分布的單元總數；Si為研究區內含有評價因素Xi的單元個數；S為研究區總單元個數；Yi為致災因子指標值。組合賦權法通過將主觀賦權法與客觀賦權法得到的權重值進行耦合，能夠避免單一評價方法的缺點和不足。利用乘數歸一法（公式（5））將層次分析法和熵權法得到的權重值進行耦合，得到各評價指標因子的綜合權重值。

式中F總表示評價單元的綜合權重值，Wj表示層次分析法求得的各評價因子的權重值；Vj表示熵權法求得的各評價因子的權重值。

最后，利用定性和定量的方法驗證本文方法的可靠性和有效性，根據文獻[21,27]通過ArcGIS 配準獲取到研究區的地質災害點，定性角度：選取姑海村、大白泥溝和大營盤3 個區域進行對比，定量角度：運用 AUC 值和準確率對3 種不同危險性評價方法得出的結果進行驗證。ROC（Receiver Operating Characteristic Curve）曲線是目前應用最廣泛的驗證評價結果的方法之一，它表示擬合數據與實際數據的關系AUC（Area under Curve）可以反映和比較模型的評價預測精度，AUC 即ROC 曲線下的面積值介于0 到1 之間，AUC 越接近1 說明效果越好。通過SPSSAU 軟件和Excel 計算得到3 種方法對應的AUC值和準確率。

2 結果與分析

2.1 地表形變信息及地質災害識別

研究區域2018 年4 月25 日至2019 年4 月20 日雷達視線方向（Line of Sight，LOS）的形變速率圖如圖3、圖4 所示。從圖3 中可以看出，采用升軌獲取的形變主要分布在小江東側，最大形變速率為?162.05 mm/a，降軌獲取的形變主要分布在小江西側西北方向，最大形變速率為?120.43 mm/a。使不同軌道得到的形變結果不同，是由于升軌數據飛行方向大致從南到北，雷達視線方向位于右側，能夠很好地將峽谷兩側由西向東的地表形變監測出來，相反，降軌數據的飛行方向與之相反，能夠將峽谷兩側由西向東的地表形變監測出來，利用不同軌道數據可以互補，使得監測結果更為準確全面，避免單一軌道帶來的幾何畸變等問題。

借助谷歌衛星影像對獲取的形變區域進行人工識別輔助識別，結合形變范圍、高程、坡度、植被覆蓋和坡體是否具有地質災害特征等進行識別，避免過度依賴形變結果導致的誤判等問題。識別結果如圖3、圖4 所示。為更為直觀地反映形變速率的實效性和準確性，將部分識別出來的地質災害區形變速率與三維影像疊置（藍色實線指向的圖），從圖中可以看到，利用升降軌對高山峽谷區地質災害識別具有良好的互補優勢，采用升軌識別的地質災害主要分布于小江沿線東側，利用降軌得到的結果相反。識別并不要求一次性將所有地質災害識別出來，僅需要盡可能均勻地識別出地質災害，為后期利用PSO-BP 模型做準備。本次識別的地質災害包含了滑坡、泥石流、崩塌和潛在災害點等。

通過實地勘察發現，識別出來的結果與實地勘察結果相吻合。如圖3、圖4 所示，灰色實線指向的圖片為實地調研勘測到的，可以看出識別出的地質災害多以滑坡、泥石流、危巖和崩塌為主。研究區地層巖性從元古界昆陽群至新生界第四系均有分布，包括沉積巖、巖漿巖、變質巖三大類。識別出來的地質災害主要巖組類型為松散巖組。主要位于小江斷裂帶兩側，斷裂兩側巖石擠壓極為強烈，常形成數百米寬的碎粉巖、碎裂巖或糜棱巖帶，兩側碎裂巖帶寬約1 000～7 000 m，鈣質膠結，斷面傾向東，傾角為55°上下。小江深大斷裂為境內地質主構造帶，沿斷裂擠壓劇烈，巖層破碎，東支為重力負異常區，西支便于地應力集中，形成地殼脆弱地帶，歷史上具有長期活動的趨勢，致使該區域地質災害多發。

2.2 格網單元劃分及評價因子選取

依據公式（2），經計算網格單元大小應選擇30 m×30 m，考慮到制圖需求以及自身電腦性能等問題，最終選擇網格單元大小為100 m×100 m，共43 505 個網格單元。選取的部分地質災害區和非地質災害區，其中地質災害區格網單元為2 094 個，非地質災害區2 662 個，利用灰色關聯分析，得到灰色關聯排序，按關聯度從大到小依次為降雨量（0.986）、地層巖性（0.982）、土地利用（0.981）、坡向（0.973）、距斷裂帶距離（0.972）、降軌形變速率（0.966）、高程（0.959）、距河流水系距離（0.953）、歸一化植被指數（0.934）、升軌形變速率（0.926）、曲率（0.917）、坡度（0.902），關聯度最高和最低的評價因子分別是降雨量和坡度，其關聯度分別為0.986 和0.902。可以看出，選取的12 個評價因子與是否為地質災害區具有高度關聯。但還不能選取出最終所需要的評價因子，為此，運用SPSSAU 軟件中的因子分析方法，通過因子分析來確定上述12 個因子之間是否存在共性，分別對12 個因子進行不同組合，當選取12 個因子進行組合時，其檢驗統計量（Kaiser-Meyer-Olkin，KMO）最佳。驗證了選用以上12 個因子作為評價因子具有有效性和可靠性。

2.3 模型構建及危險指數評定

2.3.1 模型構建及訓練驗證

將識別出來的災害和非災害點集區域定量化處理，用“1”代表高危險區，“0”代表低危險區。根據多次試驗得到隱含層為10 時結果最佳，構建的PSO-BP 模型訓練參數如下：動量因子0.9、學習速率0.05、訓練目標0.001、訓練次數500 000、進化代數80、學習因子C1=2.8，C2=1.3，最大速度0.2、種群規模25。對隨機抽取的556個點進行驗證，其輸出值與期望值如圖5 所示。

從圖5 中可以看到，驗證結果總體可靠，雖然存在少量“異常”數據（紫色區域），并不影響其整體的預測[28]。利用構建的PSO-BP 模型得到高、低危險區的危險性指數[11]，在（0.85，1.02]為高危險區，在區間（?0.15，0.25]為低危險區。選用均方誤差（Mean Squared Error，MSE）和平均絕對百分比誤差（Mean Absolute Percentage Error，MAPE）作為模型的評價指標，最終計算出MSE和MAPE 分別為0.002 和1.01%，驗證了該模型精度可靠，效果良好，保存訓練好的網絡。

2.3.2 地質災害危險性評價

將保存好的網絡運用于整個研究區，得到每個格網單元對應的危險性指數，按照文獻[9-10]的危險性分級方式分為5 個等級，采用自然間斷點分級法結合專家參與進行危險性分級處理，將其劃分為[?0.92，?0.15]、（?0.15，0.25]、（0.25，0.85]、（0.85，1.02]、（1.02，1.5]共5個等級，分別對應極低危險區、低危險區、中危險區、高危險區、極高危險區。最終得到SBAS-InSAR 和PSO-BP模型下的研究區地質災害危害性評價區劃圖如圖6a 所示。利用信息量法和組合賦權法得到研究區地質災害危險性評價，如圖6b 和6c 所示。

選取姑海村、大白泥溝和大營盤3 個區域，如圖6 所標示位置，采用本文方法、信息量法[7]和組合賦權法[27]進行定性比對可知姑海村分別為高危險區、高危險區、極高危險區；小白泥溝分別為極高危險區、極高危險區、低危險區；大營盤分別為低危險區、高危險區、極高危險區。大營盤屬于村落，其危險性一般為低危險地區，但信息量法和組合賦權法得到的評價結果均為極高危險區，本文方法判別為低危險區。姑海村屬于村落，其危險性一般為低危險地區，但信息量法和組合賦權法得到的評價結果均為極高危險區，本文方法判別為高危險區，致使本文方法判別為高危險區的原因，可能是由于姑海村位于白泥河泥石流溝對岸，所以其危險性較高。小白泥溝是典型的暴雨型泥石流溝，既有滑坡也有泥石流，其周圍應為極高危險區，本文方法和信息量法均判為極高危險區，但組合賦權法為低危險區。

運用 AUC 值和準確率對3 種不同危險性評價方法得出的結果進行驗證。如表2 所示。

表2 3 種方法AUC 值和準確率Table 2 AUC (Area Under Curve) values and accuracy of the three methods

從表2 可以看出，3 種方法均能有效的對地質災害進行評價，相較于信息量法，組合賦權法和本文方法效果更優，本文方法優于其他兩種方法。信息量法和組合賦權法計算量較大，給危險性評價增加了大量的專家參與與統計。本文方法從預測角度構建危險性評價，有效減少了計算量，其精度也較為可靠。

4 討 論

目前，雖然國內外學者對于地質災害危險性評價的研究已經取得大量成功范例，但是還存在地質災害數據源時效性差以及需進行大量的專家參與統計等弊端。針對地質災害數據源時效性差這一弊端，本文利用InSAR、谷歌衛星影像和其他輔助數據對高山峽谷地區地質災害進行識別，利用升降軌獲得的InSAR 形變速率，能夠識別出高山峽谷區不同方向的已有活躍滑坡和潛在滑坡，利用谷歌衛星影像和其他輔助數據可以克服InSAR 識別結果不全和不準確。針對需大量專家參與統計這一弊端，本文提出利用PSO-BP 模型從預測角度，通過不同格網單元評價因子的輸入來學習，得到對應格網單元的數值（危險性指數），有效的避免了大量的專家參與統計。

采用信息量法和已有學者利用組合賦權法獲取的危險性評價圖存在在部分評價錯誤的地方，導致評價錯誤的原因是，兩種方法均是基于已有災害點數據作為數據源進行危險性評價，其時效性差，存在以前是滑坡或泥石流等地質災害區域，可能已經被治理好，如果再將其定義為地質災害點，容易導致評價結果不準確，數據源的精確性直接影響著最終的評價結果。本文所用方法從定性和定量角度驗證了其評價的有效性和準確性，所用方法從災害點數據源獲取角度就有本質區別，利用該方法可以識別出已有危險和潛在危險區域，提高了災害數據源的時效性，同時運用PSO-BP 方法從預測角度構建危險性評價模型，避免了大量的評價因子權重和信息量統計，有效減少了專家參與的計算量。

5 結 論

本文通過分析現有地質災害評價模型存在的弊端，提出利用SBAS-InSAR 和PSO-BP 算法對高山峽谷區地質災害危險性進行評價，經過對研究區試驗，得出以下結論：

1）與現有地質災害數據收集方法相比，利用SBAS-InSAR 技術可以獲取到研究區地表形變量，結合高分辨率影像等輔助信息，可有效識別出活躍的泥石流、滑坡、崩塌和潛在地質災害。解決了現有地質災害數據時效性差、不準確等弊端。

2）在高山峽谷地區，利用升降軌結合的方式對研究區地質災害進行識別，其識別結果更為準確性和全面，避免了單一軌道存在SAR 成像幾何畸變造成部分地質災害不能識別、識別結果不全面等問題。同時，通過試驗驗證了升降軌形變速率可作為地質災害危險性評價的重要評價因子。

3）PSO-BP 模型適用于地質災害危險性評價，通過試驗驗證了該模型精度可靠，效果良好，其MSE 和MAPE分別為0.002 和1.01%。運用該模型避免了傳統地質災害評價需進行權值統計分級等弊端，從預測角度解決了過多的專家決策參與。

4）通過對比可知，本文方法相比信息量法和組合賦權法在地質災害危險性評價中更具有準確性和有效性，利用本文方法能夠為防災減災事業及政府部門決策提供一定參考價值。5）綜合考慮選擇信息量法和組合賦權法進行定量和定性比較，比較結果表明，本文方法有效的提高了地質災害危險性評價的準確度，信息量法、組合賦權法和本文方法的AUC（Area Under Curve）值分別為0.694、0.721、0.785，準確率為73.3%、76.2%、79.8%。