基于頻率比-投影尋蹤模型的渝東北三峽庫區滑坡敏感性制圖

毛華銳， 孫小飛， 周穎智

(1. 重慶市規劃和自然資源信息中心， 重慶 401120； 2. 西南交通大學地球科學與環境工程學院， 成都 610031; 3. 四川省南充市自然資源與規劃局， 南充 637000)

中國是山地大國，山地特有的能量梯度在降雨和人為活動的影響下，誘發了大量的地質災害，使中國成為世界上地質災害頻發且最嚴重的國家之一[1]。地質災害的頻發造成了巨大的人員傷亡、財產損失和生態破壞，嚴重威脅山區人民生命財產與工程建設安全，制約山區資源開發與經濟發展[2-4]。渝東北是三峽庫區的核心區域，也是長江流域地質災害預警高度關注的地區之一[5]。自蓄水以來，三峽庫區內地質災害隱患明顯增多，具有“點多、面廣”的特點，其中滑坡災害尤為發育[6]。因此，對渝東北三峽庫區的滑坡預測成為災害管理中亟待解決的問題。

滑坡敏感性制圖是一種減少、減輕未來滑坡造成損失的有效方法。滑坡敏感性可理解為滑坡發生的空間概率，應用適當的方法繪制滑坡敏感性分區圖，找出滑坡的高發區，有助于揭示其成因并預測未來發生的可能性，從而加強關注和監測以減少滑坡給人類生命財產造成的損失[7]。近年來，遙感和地理信息系統(geographic information system，GIS)技術為處理大型和復雜的空間數據提供了系統、快速和出色的配置，為大范圍的滑坡敏感性制圖提供了機會[8]。目前，基于遙感和GIS技術的機器學習模型已廣泛應用于滑坡敏感性制圖。例如，齊信等[9]利用GIS技術和頻率比繪制了三峽地區秭歸向斜盆地的滑坡敏感性圖，并通過滑坡敏感性面積累計百分比曲線圖驗證了結果準確度和可靠性。Sun等[10]基于貝葉斯優化算法優化超參數，建立了一個效率和、精度更高的隨機森林滑坡敏感性評估模型，并對模型進行了可靠性評價和應用驗證。高克昌等[11]利用信息量模型評估了萬州地區的滑坡危險性，認為信息量模型能夠很好地為滑坡的敏感性研究服務，可用來解決以往滑坡危險性評價中效率低、精度差的問題。黃發明等[12]將灰色關聯度模型用于浙江南田地區的滑坡敏感性評價，并認為灰色關聯度模型的預測精度略高于支持向量機模型且建模過程較為簡單。Sahana等[13]評估了頻率比、模糊邏輯和邏輯回歸模型在滑坡敏感性制圖中的有效性，認為模糊邏輯模型在評估滑坡敏感性方面更為有效。上述研究科學地預測了滑坡發生的空間概率和發育規律，為滑坡的防治提供了科學依據。總的來看，不同方法各有特點，在進行滑坡敏感性評估與制圖時應根據研究區的實際情況和所建立的評價指標選取合適的評價模型。

進行滑坡敏感性制圖需要結合歷史滑坡頻率，并盡可能避免在分配指標權重時的主觀性/不確定性[14-15]。頻率比可用于分析歷史滑坡和指標之間的相關系[16]；投影尋蹤模型是分析和處理非線性非正態高維數據的有效方法，可有效排除人為因素干擾[17]。現將討論基于頻率比和投影尋蹤模型的滑坡敏感性制圖，并將其應用到渝東北地區。目前，投影尋蹤模型已被廣泛用于區域環境評價研究中[18]，但在滑坡敏感性制圖中的應用還比較少見。因此，現擬利用頻率比和投影尋蹤模型繪制渝東北三峽庫區的滑坡敏感性圖，這對于當地結合國家長江經濟帶發展戰略制定發展規劃、開展地質災害綜合防治，打造長江上游生態保護屏障具有重要的指導意義。

1 研究區概況

渝東北地區地處渝鄂川陜四省市交界地帶，介于東經107°30′40″～110°12′12″，北緯30°02′20″～31°44′00″，面積30 698 km2，是三峽庫區的核心區域(圖1)。區內主要包括萬州、開州、豐都、忠縣、云陽、奉節、巫山、巫溪、石柱9個區縣，是國家重點生態功能區和農產品主產區、長江流域重要生態屏障和特色經濟走廊[19]。研究區屬亞熱濕潤季風氣候，汛期雨量充沛，多災害性天氣和暴雨。區內地形為中低山區，包括中切割低山、中山以及深切割中山，最低海拔73 m，最高海拔2 741 m。區內地質構造復雜且地質環境條件脆弱，是長江流域地質災害高發區之一。

圖1 研究區位置圖Fig.1 Location map of the study area

2 評價指標與數據源

2.1 評價指標

滑坡的發生是地質、地形、氣象和人類活動等許多因素綜合作用的結果。根據研究區的地質環境特點，參照已有研究[20-22]，建立滑坡敏感性評價指標體系。評價指標的選取原則如表1所示。

2.2 數據來源與數據處理

根據評價指標，收集的數據包括歷史滑坡清單、中分辨率成像光譜儀(moderate-resolution imaging spectroradiometer，MODIS)數據、土地利用、數字高程模型(digital elevation model，DEM)、降水、人口分布、地質圖和其他數據。數據源如表2所示。

利用ArcGIS 10.2軟件對DEM數據處理可獲得坡度、高程、地表起伏度等數據。巖性和斷層數據可利用地質圖矢量化獲得。采用最大值合成法[32]將MOD13Q1數據(2020年度)合成最大歸一化植被指數(normalized difference vegetation index，NDVI)分布數據，并利用像元二分模型[33]對研究區的植被覆蓋度進行定量估算。水系和人類工程活動是從土地利用數據中獲取的，并利用ArcGIS 10.2軟件的多重緩沖工具進行處理。年均降水是由原始數據通過空間克里金內插法產生的柵格數據。在評估分析之前需將指標進行重分類處理，結果如表3和圖2所示。

3 方法

3.1 多重共線性分析

指標的多重共線性分析是進行滑坡敏感評估前最重要的步驟之一[34]。指標之間的相關性直接影響到評估結果的準確性。因此，采用皮爾森相關系數[35](Pearson correlation coefficient，PCC)來分析指標之間的相關性。PCC的計算方法如下。

表1 評價指標的選取Table 1 Selection of evaluation indicators

表2 本研究所需數據及來源Table 2 Data and sources of this study

表3 評價指標的重分類Table 3 Reclassification of evaluation indicators

圖2 研究區評價指標分級Fig.2 Classifications of different evaluation indicators in study area

假設存在樣本數據集(Xi，Yi)=(X1，Y1)，(X2，Y2)，…，(Xn，Yn)，則可得到評價指標之間的相關系數PCC為

(1)

3.2 頻率比

在滑坡敏感性制圖中，一般認為未來可能發生的滑坡將在已發生滑坡相同的地質環境條件下發生[36]。因此，將頻率比用于推導歷史滑坡與指標之間的空間關系。頻率比FR的計算公式為

(2)

式(2)中：FR為滑坡在指標類型(或范圍)i中的發生頻率；Ni為類型(或范圍)i中的滑坡數量；N為滑坡總數；Mi為類型(或范圍)i的面積；M為研究區總面積。

3.3 投影尋蹤模型

采用GIS技術和投影尋蹤模型相結合的方法，開展研究區的滑坡敏感性制圖研究。模型構建步驟如下。

步驟1數據標準化。以頻率比分析結果為基礎，根據各指標子類的FR對評價指標進行歸一化處理，歸一化方程為

(3)

式(3)中：Gx為指標x的歸一化數據；gx為指標x中不同子類的FR值；gxmax為指標x中子類的最大FR值；gxmin為指標x中子類的最小FR值。

步驟2線性投影。設第i個樣本的第j個指標為Hij(i=1，2，…，m;j=1，2，…，n)，m為選取評價樣本個數，j為評價指標個數。設k為n維單位投影的方向向量，則Hij在一維線性空間的投影特征值Ii為

(4)

步驟3構造投影指標。滑坡敏感性制圖是以評價指標為基礎，將評價樣本進行合理分類或排序，是進行滑坡敏感性制圖的必要過程。因此，投影指標可依據分類指標構成。指標的分類或排序就是尋求多維數據在一維空間的類間距離W(k)和類內密度V(k)的散布結構同時滿足最大值。故定義投影指標P為

P=W(k)V(k)

(5)

類間距離用指標序列的投影特征值方差來表示[式(6)]，類內密度V(k)是利用投影特征值兩兩之間的距離Fil=|Hi-Hl|，i，l=1，2，…，m來表示[式(7)]。

(6)

(7)

步驟4優化投影向量。投影向量反映的是數據的結構特征，最佳投影向量就是最大限度地反映高維數據的某種特征結構。以投影指標P為依據，當P為最大值時，其對應的投影方向向量k為最優投影向量。因此，最優投影向量的計算可轉化為非線性優化問題。

(8)

步驟5綜合分析。滑坡敏感性指數(landslide susceptibility index，LSI)被認為是表征滑坡敏感性水平的一個重要組成部分，它被定義為若干特征指標加權求和。以特征指標所對應最優向量的最大值作為權重，并根據不同層次的綜合信息進行滑坡敏感性制圖，即

(9)

式(9)中：Ai為指標根據FR值標準化后的數據;Bi為對應的權重，也就是所對應最優向量的最大值。最后，采用自然斷點法將滑坡敏感性指數劃分為低敏感，較低敏感，中敏感，較高敏感和高敏感5個等級。

3.4 滑坡敏感性分布趨勢

采用滑坡敏感性綜合指數(composite index，CI)表征滑坡敏感性的分布趨勢，即

(10)

式(10)中：i為敏感性等級，n為總級別數；Ei為等級i在評級單元j中所占的面積；SFj為評價單元j的面積；Fi為等級i的分級值。根據敏感性區劃將低敏感區賦值1，較低敏感賦值2，中敏感區賦值3，較高敏感區賦值4，高敏感區賦值5。CI的值越大，表示整體滑坡敏感性越高。

3.5 驗證方法

結果驗證是滑坡敏感性制圖中必要的過程。接收者操作特征(receiver operating characteristic，ROC)曲線是一種對于靈敏度進行描述的功能圖像，被廣泛應用于包括滑坡在內的多個領域[37]。ROC曲線繪制了一組閾值或臨界值的真陽性率(true positive rate，TPR)，表明評估結果成功地定位了敏感區；而假陽性率(false positive rate，FPR)則顯示了敏感區不顯著的位置。結果的精度用ROC曲線下面積(area under the curve，AUC)來衡量。一般AUC值在大于0.5的情況下，其值越接近1，則說明模型應用越成功。將歷史滑坡數據劃分為樣本數據(歷史滑坡總數的70%)和驗證數據(歷史滑坡總數的30%)來評價結果的準確性。樣本數據用于計算結果的成功率，驗證數據用于計算結果的預測率。

4 結果

4.1 評價指標之間的相關性

指標之間相關性將影響滑坡敏感性制圖的準確性。通過皮爾森相關系數法計算了評價指標之間的相關系數，得到評價指標之間的相關系數在-0.43～0.58，以弱相關或無相關為主(圖3)。因此，在開展滑坡敏感性制圖時，無需對評價指標進行刪減。

藍色表示負相關；紅色表示正相關；圓圈的大小表示相關性大小圖3 評價指標之間的相關性Fig.3 Correlations between different evaluation indicators

4.2 評價指標與歷史滑坡的關系

選擇研究區內4 322個歷史滑坡(歷史滑坡總數的70%)作為樣本數據來探討指標與歷史滑坡之間的分布關系。基于頻率比模型的分析結果(圖4)，根據各指標子類的FR值對指標進行標準化，

為后續的利用投影尋蹤評估研究區滑坡敏感性做準備。

4.3 滑坡敏感性評估與制圖

采用遺傳算法尋找最優投影向量，優化步驟如下。

(1)隨機產生X組原始單位投影向量，根據式(4)分別計算每個組的特征值。

(2)根據式(6)和式(7)分別計算X組投影方向上的W(k)和V(k)，并帶入式(5)計算X組投影方向上的目標函數P。

(3)利用遺傳算法中的選擇、變異和交叉操作，得到多組新的投影向量。

(4)由于P是具有空間分布且連續的數據，故只需滿足P最大值越大越優的原則，從多組投影向量中選取新的投影向量。

(5)重復上述操作，直到P的最大值不再增大時，得到所對應的投影向量就是最優投影向量。

根據表4中的特征指標及其最優投影向量，帶入式(9)可得到研究區滑坡敏感性指數，并將滑坡敏感性指數劃分為低敏感、較低敏感、中敏感、較高敏感和高敏感5個等級(圖5)。

圖4 滑坡與評價指標之間的關系Fig.4 Relationship between landslide and evaluation indicators

通過對滑坡敏感性區劃統計(表5)，研究區內高敏感區面積為2 823.84 km2，滑坡1 618處，占滑坡總數的26.21%；較高敏感區5 923.30 km2，滑坡1 795處，占滑坡總數的29.07%；中敏感區9 820.13 km2，滑坡1 783處，占滑坡總數的28.88%；較低敏感區7 825.97 km2，滑坡805處，占滑坡總數的13.04%；低敏感區4 304.35 km2，滑坡173處，占滑坡總數的2.80%。通常情況下，在滑坡敏感性分級結果中敏感性越高的區域滑坡密度約大，則表明分級結果越合理。在劃分結果中，存在84.16%的滑坡分布于中、較高以及高敏感區，且滑坡密度隨滑坡敏感程度的增加而增加，這表明研究區域的滑坡敏感性區劃結果是有效合理的。

表4 特征指標及其最優投影向量Table 4 Characteristic indicators and its optimal vector

圖5 研究區滑坡敏感性區劃圖Fig.5 Landslide susceptibility zoning map of study area

表5 研究區滑坡敏感性區劃結果Table 5 Results of landslide susceptibility zoning

4.4 驗證結果

根據ROC曲線驗證結果(圖6)，結果的成功率AUC值為0.844，標準誤差0.003，95%的置信區間內AUC值的上限和下限分別為0.838和0.850，顯著水平P<0.001。預測率AUC值為0.763，標準誤差0.005，95%的置信區間內AUC值的上限和下限分別為0.753和0.773，顯著水平P<0.001。參考Yesilnacar等[38]對AUC的分類，基于頻率比-投影尋蹤模型的渝東北三峽庫區滑坡敏感性制圖具有相當高的精度。

圖6 評估結果的成功率和預測率ROC曲線Fig.6 ROC curve of evaluated success rate and prediction rate

5 討論與結論

5.1 討論

在遙感和GIS技術支持下，基于頻率比和投影尋蹤模型的滑坡敏感性制圖是一項有利于防災減災的工作。其從空間上評估了滑坡發生的可能性，評價結果對減少人民生命財產損失和促進社會發展都具有重要意義。根據滑坡敏感性分布趨勢分析，同一指標在不同范圍或類型中的滑坡敏感性綜合指數存在空間差異性(圖7)。

從地形地貌來看，較高和高敏感區主要分布在海拔<1 000 m的區域；由于隨著坡度的增大，斜坡巖土體抗能力逐漸減弱而下滑力增大，但當坡度達到一定角度，坡體上的松散固體物質切向分力就可以克服摩擦力而向下運動，因此高敏感區的坡度主要在10°～35°[39]。從地質因素來看，地層巖性是控制地質災害的主要因素，高敏感區主要分布在頁巖、泥巖、砂頁巖等軟巖區域，這可能是軟巖具有膨脹性、崩解性、流變性和易擾動性等特性，使得軟巖區斜坡的穩定性較低[40]。高敏感區的植被覆蓋度大多低于60%，可能的原因是高植被區的保持水土、穩固邊坡作用更強。人類工程活動日益增強，改變坡體自然結構，松動巖石，增大坡高坡度，破壞斜坡自身平衡，導致高敏感區主要分布在距人類工程活動1 500 m之內的區域。此外，距離水系越近，滑坡敏感性越高，主要是水流的下切作用使得坡體原來平衡的應力狀態遭到破壞[41]。降水是滑坡產生的重要誘發因子，但由于研究區年均降水豐富但差異不大(1 000～1 200 mm)，以至研究區的滑坡敏感性分布與年均降水沒有明顯的相關性。

此外，將投影尋蹤模型應用于滑坡敏感性評估建模，并完成了渝東北三峽庫區滑坡敏感性制圖與評估。利用ROC曲線證明了應用投影尋蹤模型進行滑坡敏感性制圖與評估是可行的，而且可以定量地描述每個指標對滑坡敏感性的影響和貢獻。投影尋蹤模型是尋找出反映原始高維數據結構特征的投影，以達到研究和分析高維數據的目的。和其他方法相比，該模型不需要人工假設，除了自動設置評價指標權重之外，還避免了嚴重的信息丟失，具有較強的魯棒性和客觀性[42]。

5.2 結論

在遙感和GIS技術支持下，選取了高程、坡度、地形起伏度、水系距離、斷層距離、地層巖性、植被覆蓋度、年均降水、距人類工程活動距離9個指標，結合頻率比和投影尋蹤模型繪制了渝東北三峽庫區的滑坡敏感性圖，并將研究區劃分為低敏感區、較低敏感區、中敏感區、較高敏感區和高敏感區。研究區內滑坡敏感性以中低敏感性為主，制圖結果能夠較好地反映渝東北三峽庫區的滑坡災害狀況。

根據敏感性制圖結果，研究區的較高、高風險區約占22%，主要分布在忠縣、萬州以及巫山等區縣的人類工程活動集中區。同一指標的不同范圍或類型對研究區滑坡敏感性的影響存在空間差異。地形地貌是影響研究區滑坡敏感性的主要環境因素，人類工程活動是影響滑坡敏感性的主要誘發因素。該研究結果可為渝東北三峽庫區的滑坡防治提供參考依據，對區內災害管理具有重要意義。