基于CF-LR組合模型的地質災害易發性評價
——以云南省寧洱哈尼族彝族自治縣為例

秦 紅 富，談 樹 成，2，施 旖 奇，李 紅 梅，汪 旭

(1.云南大學 地球科學學院，云南 昆明 650500； 2.云南省地理研究所，云南 昆明 650223)

0 引 言

地質災害作為一種破壞性的地質事件，制約著人類社會的可持續性發展，對人類的生命財產和生存環境構成嚴重威脅，引發了社會各界的關注，已成為災害學中的研究熱點[1]。地質災害多發生于坡度較大、斷裂帶密集、河網密布、土質疏松、巖土體脆弱、植被稀疏的斜坡上，此外強降雨和過度的人類活動也是滑坡、泥石流等地質災害形成的因素[2]。中國地質災害主要分布在西南、西北地區，該地區區域地質地理環境復雜多變使得地質災害具有類型多、數量多、頻率高、分布廣、危害大等特點，其中尤以四川省、云南省為代表。采用合適的評價模型或方法計算地質災害發生的可能性，并對地質災害易發性結果分區，對區域防災減災具有重要意義。

國內外有關區域地質災害易發性評價的研究工作已取得大量的成果。常用的易發性分析方法有層次分析法[3]、隨機森林法[4]、信息量法[5]、支持向量機[6]、邏輯回歸[7]、確定性系數法[8]、最大熵模型[9]、人工神經網絡和決策樹[10]等。結合各方法的優缺點，進行模型組合，可以有效提高模型精度，采用不同方法結合進行地質災害易發性研究已經成為研究的熱點和趨勢。與此同時，GIS技術的飛速發展為進行地質災害易發性區劃深入研究提供了一個卓有成效的技術平臺與研究途徑[11]。

因此，以云南省寧洱哈尼族彝族自治縣為研究區(簡稱“寧洱縣”)，GIS為技術平臺，采用確定性系數和邏輯回歸模型(CF&LR)相結合的方法開展以滑坡、泥石流、崩塌等為主體的地質災害易發性評價，以期為寧洱縣的防災減災、區域防治、地質災害風險預測等工作提供參考。

1 研究區概況

寧洱縣地處普洱市中部，東經100°42′～101°37′，北緯22°40′～23°36′。研究區東西最大橫距91 km，南北最大縱距101 km，轄6鎮3鄉，國土總面積3 669.77 km2，地處無量山脈西南部，屬紅河水系支流李仙江與瀾滄江的分水嶺地區，雨量充沛，年平均氣溫18.2 ℃，年平均降雨量為1 414.9 mm，夏季降水集中而冬季降水較少。海拔高差最大達2 300 m，相對高差一般500～1 000 m，總體屬中淺切割的中山-低中山地形。研究區山地面積占總面積的96.77%，公路建設和礦山開發等人類社會工程活動對生態環境有一定影響和破壞，誘發了一些地質災害。寧洱縣地質災害類型有滑坡、不穩定斜坡、崩塌、泥石流、地面塌陷等，是云南省滑坡、泥石流等地質災害較為嚴重的地區之一。

2 材料和方法

2.1 數據來源與預處理

災害點數據來源于項目組成員根據2018年高分辨率遙感影像和Google影像聯合解譯得到，并于2019年3月前完成野外調查驗證，其中滑坡36處、崩塌38處、泥石流5處(見圖1)；DEM數據來源于地理空間數據云平臺，用于提取高程、坡度、坡向；降雨量數據來源于中國科學院資源環境科學數據中心；斷裂帶和巖性數據來源于1∶50 000中國地質圖。河流、道路及基礎地理數據來自國家基礎地理信息中心。利用ArcGIS工具將所有的柵格數據的坐標統一為WGS_1984_UTM_Zone_47N，并對各個評價指標進行重分類。

圖1 寧洱縣地理位置及主要災害點分布Fig.1 The geographical location and distribution of major geological disaster points in Ninger County

2.2 評價方法與評價單元的確定

2.2.1評價方法

(1) 確定性系數模型。CF模型是由Shortliffe[12]提出的函數，用來表示概率。后來，Heckerman[13]完善了該模型。CF 模型計算過程如下：

(1)

式中：ppa是地質災害事件在地質環境因子a中發生的概率[0，1]，即地質環境因子a在特定單元中地災點個數與地質環境因子a總面積的比值。pps是地質災害在整個研究區中發生的先驗概率，為地災點總數與研究區總面積之比[14]。CF取值區間為[-1，1]，當CF大于0時，代表發生地質災害可能性高，值越接近于1，發生地質災害得可能性越高；當CF小于0時，代表發生地質災害的可能性較低，值越接近于-1，發生地質災害得可能性越低。

(2) 邏輯回歸模型。在地質災害易發性評價中將評價指標看作獨立變量，地質災害是否發生看作二元因變量(0指代地質災害發生，1指代地質災害不發生)[15]。邏輯回歸模型描述的是二元因變量和獨立變量xi(i=1，2，3….n)的關系。其中，獨立變量不必滿足正態分布，LR模型對數據類型要求較低，可以減少難以量化的數據對研究結果的影響，同時可以將研究步驟化繁為簡。函數如下式：

(2)

式中：P表示地質災害發生的概率，范圍為[0，1]，值越大發生地質災害的概率越大；α為常數項；xi(i=1，2，3，…，n)為因子i中各分類級別的CF值；β為各xi對應的回歸系數[16]。

2.2.2評價單元的確定

合理的評價單元，決定數據的預處理格式和評價結果精度。總的來說，地質災害易發性評價中較常用的評價單元有行政單元、規則柵格單元、斜坡單元等[17]。一方面，研究區地質地貌類型復雜多樣、河流眾多，降雨量充沛；另一方面考慮到比例尺和數據精度，以及研究區地質災害主要以滑坡、崩塌、泥石流為主，本文采用規則柵格單元作為地質災害易發性評價單元，目前較為常用的基于專家經驗[18]的計算公式為

Gt=7.49+0.0006t-2.0×10-9t2+

2.9×10-15t3

(3)

式中：Gt指合適的柵格單元大小；t指等高線精度的分母值。研究區所使用的DEM數據為1∶50 000，則由式(3)計算得到合適的柵格單元大致為30 m，利用ArcGIS工具將研究區劃分為4 075 195個30 m×30 m的評價柵格單元。

2.2.3評價因子的選取和分級

(1) 降雨量。在誘發地質災害的眾多因素當中，降雨發揮著十分重要的作用，為地質災害發育提供了外動力條件。本文共收集近10 a寧洱縣降水量空間插值數據，利用ArcGIS中的柵格計算器工具得到寧洱縣年平均降雨量數據，可以發現近10 a寧洱縣平均降雨量最小值為1 232 mm，最大值為1 798 mm。然后利用自然間斷法將計算結果重分類分為1 230～1 380，1 380～1 458，1 458～1 537，1 537～1 628，1 628～1 800 mm共5個級別(見圖2(a))。將降雨量柵格數據轉換為面文件與地災點數據建立空間連接，進行不同分級范圍內的地質災害數量統計，繪制直方圖(見圖3(a))分析不同分級降雨量與災害點相對密度的關系。結果顯示，降雨量在1 458～1 537 mm范圍內，單位面積災害點比例最大，說明該區域發生地質災害的可能性較高。

(2) 距斷裂帶距離。斷裂構造一直被認為是地質災害發生的重要影響因素。在斷裂帶附近，發育地質災害的斜坡常較為集中的發生，一旦受到誘發因素(降雨、人類工程活動等)的作用就更容易發生地質災害。本次研究在ArcGIS中對寧洱縣斷裂帶數據建立緩沖區，分別構建0～500，500～1 000，1 000～1 500，1 500～2 000，2 000～2 500，>2 500 m共6個等級(見圖2(b))，并對不同等級范圍的地災數量進行統計。然后繪制直方圖(見圖3(b))分析距斷裂帶不同距離與災害點相對密度的關系。結果顯示，在距斷裂帶0～500 m和500～1 000 m范圍內災害點相對密度較其它分級范圍大得多，研究區地質災害主要集中在距斷裂帶0～1 000 m范圍內。

(3) 距道路距離。在誘發地質災害發生的眾多因素當中，人類活動也起著十分重要的作用，包括切坡建路、采石采礦、鄉鎮建設等。一般來說僅考慮單因素對地質災害的影響時，距道路距離越近相對較容易發生地質災害。基于寧洱縣鄉道、省道、國道、高速矢量數據，利用ArcGIS中多環緩沖區工具對道路數據建立緩沖區，得到0～500，500～1 000，1 000～1 500，1 500～2 000，2 000～2 500，>2 500 m共6個等級范圍(見圖2(c))，并對不同等級地災數量進行統計，并繪制直方圖(見圖3(c))分析距河流不同距離與災害點相對密度之間的關系。結果顯示：在距道路距離0～500 m范圍內災害點相對密度為0.034，僅次于距道路距離2 000～2 500 m范圍的0.04，進一步證實了地質災害的發生不僅受到距道路距離的影響還受到其他因素的影響。

(4) 距離河流距離。河流對地質災害的影響在于對河流兩岸斜坡前緣的侵蝕作用，增加了斜坡的臨空面，導致河流邊緣日益陡峭，破壞其穩定性。利用ArcGIS中多環緩沖區工具對河流矢量數據建立緩沖區，得到0～200，200～400，400～600，600～800，800～1 000，>1 000 m共6個等級范圍(見圖2(d))。并對不同等級地災數量進行統計，并繪制直方圖(見圖3(d))分析距道路不同距離與災害點相對密度之間的關系。結果顯示：研究區地質災害主要發生在距離河流0～200 m范圍內，隨著距河流距離的增加，災害點相對密度逐級遞減，說明距離河流越近越容易發生地質災害。

(5) 巖石硬度。巖土體是地質災害的物質組成，也是地質災害發生的重要因素之一。巖石的類型和軟硬程度不同，其抗風化和抗腐蝕能力不同，導致其抗滑動能力也不同。基于寧洱縣地層矢量和巖土體堅硬程度，按照一般巖石堅硬程度分類標準，將巖石劃分為堅硬巖類、較堅硬巖類、較軟巖類、軟巖類共4類(見圖2(e))，統計不同巖石硬度范圍災害點的數量并繪制直方圖(見圖3(e))，分析不同巖石硬度與災害點相對密度之間的關系。結果顯示，較堅硬巖類的災害點相對密度最大，但是僅有包含9個災害點，主要是因為已有地質災害點多分布于堅硬若和較軟巖類交界處導致的，較軟巖類的相對密度次之，但包含了69個地災點，說明該區域地災多發育于相對密度第的較軟巖類，而不是相對密度最大的較堅硬巖類。

圖2 各評價因子分級Fig.2 Grading diagram of each evaluation factor

圖3 各評價因子與災害點相對密度關系Fig.3 The relative density diagram of each evaluation factor and geological disaster points

(6) 歸一化植被覆蓋指數(NDVI)。植被影響著土壤的抗侵蝕能力和地表徑流，進而影響巖土體。因此，NDVI也是地質災害發生的影響因子之一。基于2018年的Landsat8數據，利用ArcGIS中的柵格計算器構建NDVI計算模型計算得到寧洱縣NDVI數據。利用ArcGIS中的重分類工具將計算結果分為-0.2～0，0～0.2，0.2～0.4，0.4～0.6共5類級別(見圖2(f))，統計不同NDVI分級范圍災害點的數量并繪制直方圖分析災害點相對密度與不同NDVI分級范圍之間的關系。結果顯示：NDVI值處于0～0.2之間災害點相對密度值最大，較容易發生地質災害。

(7) 高程。高程對地質災害的作用，一方面是控制巖土體應力值的大小，巖土體的應力值和勢能會隨著所處高程的增高而增大，易在誘發因素的影響下發生地質災害；另一方面是不同海拔高度的地區、氣候、植被類型、降雨量、人類活動均有所不同，可能提供不同的孕災環境。寧洱縣整體高程相差較大，最低高程為534 m，最高為2 830 m。基于ArcGIS對DEM數據進行處理后，將高程按照534～1 000，1 000～1 450，1 450～1 900，1 900～2 350，2 350～2 830 m進行分級(見圖2(g))，統計不同高程分級范圍內災害點數量并繪制直方圖分析災害點相對密度與不同高程之間的關系。結果顯示：高程范圍在1 450～1 900 m之間災害點相對密度值最大，高程范圍在1 000～1 450 m之間災害點相對密度值次之，地質災害主要發生在高程范圍為1 000～1 900 m范圍內，該范圍內災害點個數占總災害點個數的93.67%，主要是因為在該高程范圍內交通線路密集、人類活動頻繁。

(8) 坡度。坡度一般通過決定斜坡體應力大小和方向，影響地表徑流和物質移動等進而控制斜坡巖土體的穩定性。利用ArcGIS中的坡度工具，基于寧洱縣DEM數據計算研究區坡度，得到坡度數據。結果顯示，研究區坡度最大值為69.2°。并利用重分類工具對坡度數據進行分級，劃分為0°～15°，15°～30°，30°～45°，45°～60°，60°～75°共5個坡度等級范圍(見圖2(h))，統計不同坡度分級范圍內災害點數量并繪制直方圖(見圖3(h))分析災害點相對密度與不同坡度范圍之間的關系。結果顯示：在坡度范圍為45°～75°范圍內沒有地災的發生，坡度范圍在30°～45°災害點相對密度值最大地災災害較易發生，在坡度范圍15°～30°內災害點數量最多。

(9) 坡向。坡向對地質災害的影響主要體現在不同坡向的斜坡受到的太陽輻射強度和光照時長有所區別，因此，同一山體的局部氣候會有所不同，也就是溫度、降雨量有所不同，進而影響研究區的植被覆蓋度，巖石風化速率、土壤濕度等。研究表明，在相同巖土體條件下，陽坡面發生地災的概率大于陰坡面，因為陽坡面具有巖石風化速率快、降雨量多等特點。基于寧洱縣DEM數據，利用ArcGIS計算坡向，得到坡向數據，并對其進行劃分，平地(-1°)、東北向(22.5°～67.5°)、正東向(67.5°～112.5°)、東南向(112.5°～157.5°)、正南向(157.5°～202.5°)、西南向(202.5°～247.5°)、正西向(247.5°～292.5°)、西北向(292.5°～337.5°)、正北向(337.5°～22.5°)共9個坡向等級范圍(見圖2(i))。統計不同坡向災害點數量并繪制直方圖(見圖3(i))分析災害點相對密度與不同坡向的關系。結果顯示：平地(-1)內沒有地質災害發生，災害點相對密度值最大的坡向為西南向(202.5°～247.5°)，西南向較其它坡向更容易發生地質災害。

3 易發性評價流程

基于各評價指標因子的柵格數據利用CF模型計算出坡度、坡向、降雨量、巖石硬度、距河流距離、距道路距離、距斷裂帶距離、NDVI、高程9個評價指標因子不同分類級別的確定性系數即CF值，然后將9個評價指標因子各分類級別的CF值看作獨立變量，將地質災害是否發生看作因變量(0表示不發生，1表示發生)。利用python中的邏輯回歸分析方法進行邏輯回歸，得到各評價指標因子的回歸系數，看作各指標的權重，同時進行顯著性檢驗，求得邏輯回歸方程。最后利用多值提取至點工具對研究區4 075 195個格網坡度、坡向、降雨量、巖石硬度、距河流距離、距道路距離、距斷裂帶距離、歸一化植被指數(NDVI)、高程數據進行賦值并代入邏輯回歸方程，得到研究區地質災害易發性P值，最后利用等間距法進行分區即可得到研究區地質災害易發性分區結果。

3.1 計算各評價指標的CF值

在對各評價指標不同分級范圍的地質災害進行統計分析后，根據2.2節中確定性系數模型，計算得到各評價指標因子不同分級范圍的CF值，CF值越大表明地質災害發生的確定性越高，反之越低。各評價指標確定性系數計算結果見表1。

表1 各評價指標確定性系數

由表1可以看出，各評價指標的CF值除了與各分級范圍內的地質災害數量有關還與各分級范圍的面積有關。

(1) 從降雨量5個分級范圍的CF值來看，確定性系數最大的范圍是1 458～1 537 mm，CF值為0.344，確定性系數最小的范圍是1 628～1 800 mm，CF值為-0.586。表明降雨量為1 458～1 537 mm范圍時發生地質災害的確定性最高，為研究區地質災害發生的最佳孕災環境之一。

(2) 從距斷裂帶距離的6個分級范圍的CF值來看，0～500 m范圍CF值為0.331，500～1 000 m范圍CF值次之，表明距離斷層距離越近地質災害更容易發生。

(3) 從距道路距離的6個分級范圍的CF值來看，2 000～2 500 m范圍內CF值最大為0.467，0～500 m范圍內CF值次之為0.373，在0～500 m范圍內發生的地質災害的數量要多于2 000～2 500 m范圍，但0～500 m分級范圍的分級面積大于2 000～2 500 m分級范圍的分級面積，導致0～500 m范圍內CF值小于2 000～2 500 m范圍內CF值。表明2 000～2 500 m范圍地質災害的發生不只受到距道路距離的影響，還受到其他因素的共同作用。

(4) 從距河流距離的6個分級范圍的CF值來看，0～200 m范圍內CF值最大為0.487，表明距離河流越近越容易發生地質災害。

(5) 從巖石硬度的4個分級范圍的CF值來看，較堅硬巖類的CF值最大為0.304，較軟巖類CF值次之為-0.023，這是因為較軟巖類的分級面積為3 278.03 km2遠大于較堅硬巖類的293.49 km2，故即使較軟巖類發生地質災害的數量最多達到69處，但確定性系數卻較低。

(6) 從NDVI的4個分級范圍的CF值來看，0～0.2范圍內CF值最大為0.645，其余分級范圍CF值為負值，表明植被覆蓋越少的地方越容易發生地質災害，反之越不容易發生地質災害。

(7) 從高程的5個分級范圍的CF值來看，在1 450～1 900 m范圍內CF值最大為0.098，1 000～1 450 m范圍內CF值次之，1 450～1 900 m范圍內較其它分級范圍更容易發生地質災害。

(8) 從坡度的5個分級范圍的CF值來看，坡度范圍在30°～45°的CF值最大為0.047，其它坡度范圍CF值均為負值，30°～45°范圍內較其他分級范圍更容易發生地質災害。

(9) 從坡向的分級來看，坡向西南向(202.5°～247.5°)的CF值最大為0.525，西南向較其他坡向更容易發生地質災害。

3.2 評價指標獨立性檢驗

采用相關性分析進行評價指標之間的獨立性分析，剔除相關性較大的指標因子，使用python求得各因子之間相關系數矩陣見表2。

表2 評價指標間的相關系數矩陣

表2中X1，X2，X3，X4，X5，X6，X7，X8，X9分別代表巖石硬度、距斷裂帶距離、距河流距離、降雨量、距道路距離、NDVI、坡度、坡向、高程。結果顯示：各因子之間的相關系數絕對值都小于0.3，表明所選因子之間的相關性較小，9個評價指標因子都可以代入評價模型。

3.3 計算各評價指標權重

研究區一共有79個災害點，在計算得到各評價指標的CF值后，隨機選取研究區內79個沒有地質災害發生的點，利用ArcGIS中的多值提取至點工具和python工具將各評價指標的CF值賦值給這158個樣本點，然后對這158個樣本點基于邏輯回歸模型來進行邏輯回歸分析，結果見表3。

表3 邏輯回歸分析結果

結合式(2)和表3的回歸系數可以得到以下邏輯回歸方程組：

(4)

式中：P為發生地質災害的可能性，值區間為[0，1]；x1為巖石硬度中各分級范圍的CF值；x2為距斷裂帶距離中各分級范圍的CF值；x3為距河流距離中各分級范圍的CF值；x4為降雨量中各分級范圍的CF值；x5為距道路距離中各分級范圍的CF值；x6為NDVI中各分級范圍的CF值；x7為坡度中各分級范圍的CF值；x8為坡向硬度中各分級范圍的CF值；x9為高程中各分級范圍的CF值。

由各評價指標因子的回歸系數可知，代入評價模型的9個評價指標因子對研究區地質災害敏感性由高到低依次為距河流距離、距斷裂帶距離、坡向、高程、NDVI、降雨量、坡度、巖石硬度、距道路距離。

4 結果分析與驗證

4.1 評價結果

寧洱縣地質災害易發性評價結果在0.045～0.861之間。基于ArcGIS使用等間距法將易發性分區評價結果分為：極高易發區(0.657～0.861)，高易發區(0.453～0.657)，中易發區(0.249～0.453)，低易發區(0.045～0.249)4個等級(見圖4)。再結合圖3可知極高易發區主要分布在斷裂帶密度較大、路網密集、人類活動強度較大、年平均降雨量在145～1 537 mm之間的西南和中部區域。高易發區主要分布于研究區西部區域和南部區域。中易發區主要分布于研究區的東部區域以及北部部分區域。低易發區主要分布于研究區的北部區域和東南區域。結合災害點的空間分布來看，災害點越密集的區域，發生地質災害的風險就越大，整體上，研究區地質災害呈現北少南多、東少西多的特征。

4.2 評價結果驗證

得到寧洱縣評價結果后，通過兩個方面對評價分區結果進行驗證，一是通過現狀災害點空間分布對分區結果進行合理性檢驗；二是利用ROC曲線對結果進行準確性檢驗。

基于ArcGIS工具統計現狀災害點的空間分布結果判斷分區結果的合理性及可靠性，對采用CF&LR組合模型得到的研究區易發性分區結果進行統計可得表4。結果顯示：極高易發區分布災害點35個，占災害點總數的44.30%，分區面積占研究區面積的8.23%，單位面積災害點比例最大；高易發區和中易發區均分布災害點22個，占災害點總數的27.85%，但單位面積災害點密度高易發區較中易發區大；低易發區內無現狀災害點分布。現狀地質災害在極高易發區強發育，高易發區中等發育，中等易發區弱發育，低易發區內較弱發育。從驗證結果來看，本次評價分區合理。

ROC曲線是一種不受臨界約束的結果評價方法，能有效的對評價結果的準確性進行檢驗。ROC曲線以假陽性率(1-specificity)為橫坐標，真陽性率(Sensitivity)為縱坐標繪制的曲線。ROC曲線下的面積為AUC值，是衡量模型準確性的指標，AUC的取值區間為[0.5，1]，當AUC=0.5時，模型預測無效；當AUC范圍在[0.5，0.7]時，模型準確性較低；當AUC范圍在[0.7，0.9]時，模型準確性較高；當AUC>0.9時模型準確性特別高。利用python工具計算得到評價結果的ROC曲線 (見圖5)。結果顯示：AUC值為0.74，表明基于CF&LR組合模型的寧洱縣易發性評價結果準確性較高，可較為客觀地對寧洱縣地質災害易發性進行分區評價。

圖5 評價結果ROC的曲線Fig.5 Curve of evaluation results for ROC

5 結 論

(1) 研究區9個評價指標對地質災害發生的敏感性依次為距河流距離、距斷裂帶距離、坡向、高程、NDVI、降雨量、坡度、巖石硬度、距道路距離。

(2) 研究區最易發生地質災害的條件為距河流距離0～200 m、距斷裂帶距離0～500 m、坡向西南向(202.5°～247.5°)、高程1 450～1 900 m、NDVI范圍在0～0.2、降雨量范圍1 458～1 537 mm、坡度30°～35°、巖石硬度為較軟巖類、距道路距離0～500 m。

(3) 極高易發區、高易發區、中易發區、低易發區分區面積分別占研究區面積的8.23%，37.17%，43.39%，11.21%，研究區地質災害呈現北少南多、東少西多的特征。

(4) 對易發性分區結果進行合理性和準確性檢驗，從合理性檢驗來看，基于CF&LR組合模型的分區結果基本合理。從準確性檢驗來看，準確性檢驗結果顯示AUC值為0.74，表明基于CF&LR組合模型的寧洱縣易發性評價結果準確性較高。