基于證據(jù)權法的昆明五華區(qū)地質災害易發(fā)性評價

白光順，楊雪梅，朱杰勇，張世濤，祝傳兵，康曉波，孫 濱，周琰嵩

（1.昆明理工大學國土資源工程學院，云南 昆明 650093；2.自然資源部高原山地地質災害預報預警與生態(tài)保護修復重點實驗室，云南 昆明 650093；3.云南高正地信科技有限公司，云南 昆明 650041；4.云南省地質環(huán)境監(jiān)測院，云南 昆明 650216；5.云南南方地勘工程總公司昆明分公司，云南 昆明 650051）

0 引言

查明與地質災害有關的危險區(qū)域是地質災害管理的重要工作，也是促進研究區(qū)人民生活和基礎設施發(fā)展安全的重要依據(jù)[1]，基于建模評價地質災害易發(fā)性是重要而且有效的途徑。

應用經(jīng)驗式、數(shù)值模擬和統(tǒng)計方法對地質災害易發(fā)性建模和評價，已經(jīng)進行了許多研究[1-10]。其中，經(jīng)驗式方法基于現(xiàn)場觀察和專家經(jīng)驗判斷；數(shù)值模擬計算邊坡的穩(wěn)定性；統(tǒng)計方法部分基于實地觀察和專家的先驗知識，部分基于對地質災害發(fā)生的權重或概率的統(tǒng)計計算，這類方法使用統(tǒng)計技術來評估誘發(fā)地質災害的各種因素的相關作用，每個因素的重要性都是根據(jù)觀察到的與地質災害的關系來確定的。

文中使用基于貝葉斯理論的證據(jù)權法，綜合GIS 技術評價研究區(qū)地質災害易發(fā)性。證據(jù)權法是一種統(tǒng)計方法，最初應用于非空間、定量的醫(yī)學診斷，以結合臨床診斷的證據(jù)來預測疾病[11-12]。在地球科學中，該方法被廣泛應用，如：礦產(chǎn)資源潛力評估和礦床預測[13-16]，公路路基巖溶塌陷危險性評價[17]和滑坡易發(fā)性和危險性[1,3,18-23]。

1 研究區(qū)和數(shù)據(jù)

1.1 研究區(qū)概況

文中選擇云南高原滇中昆明盆地低山丘陵地帶這一云南省地質災害防治重點地區(qū)的典型代表，云南省省會昆明市的主要行政區(qū)之一，昆明市五華區(qū)作為研究對象，該區(qū)地質災害易發(fā)性評價研究具有典型代表性，可向整個云南高原昆明盆地低山丘陵區(qū)和其他低山丘陵區(qū)推廣，具有技術方法和社會經(jīng)濟意義。研究區(qū)面積381.6 km2，地勢西北高東南低，昆明盆地內地形開闊低緩，北部山區(qū)地形崎嶇，溝壑較發(fā)育。區(qū)域年降水量的80%以上集中在6—9月，年平均降水量608.4～887.0 mm。碳酸鹽巖分布最廣，約占全區(qū)面積的38.93%，其次為砂巖、泥巖、頁巖，約占23.11%，巖漿巖主要為玄武巖，約占16.95%，主要分布在昆明盆地和其他小盆地的松散碎石土體約占11.36%，石英砂巖類約占7.56%，還發(fā)育一些巖脈；斷裂構造較發(fā)育，以南北向構造為主[24-25]。

1.2 數(shù)據(jù)收集

通過地質災害風險普查獲得了研究區(qū)地質災害分布數(shù)據(jù)。根據(jù)調查分析，選擇工程地質巖組、斷裂構造、高程、坡度、坡向、坡面曲率、距公路距離和土地利用類型等8 類因素納入評價分析。地質數(shù)據(jù)收集自云南省地質局1∶20 萬昆明幅、武定幅區(qū)域地質調查報告和圖件[24-25]，12.5 m 分辨率DEM（數(shù)字高程模型）收集自ASF，道路數(shù)據(jù)收集自OSM，土地利用類型數(shù)據(jù)收集自ESA（圖1、表1）。

表1 數(shù)據(jù)簡介Table 1 Data introduction

圖1 因素基礎數(shù)據(jù)圖Fig.1 Basic data charts of factors

1.3 地質災害發(fā)育特征

現(xiàn)狀發(fā)育地質災害89 處，滑坡73 處，崩塌11 處，泥石流4 條，地面沉降1 處，為小—中型，無大型，中型14 處，小型75 處，主要分布在研究區(qū)低山丘陵地貌區(qū)，盆地內僅發(fā)育1 處（圖2）。

圖2 地質災害分布圖（底圖為高程和山體陰影渲染）Fig.2 Map of geological hazard distribution (The bottom was rendered by elevation and hillshade)

選擇指標“因子面積百分比A”“地災數(shù)百分比B”和“比率(β=B/A)”表征地質災害的空間分布特征、主控因素和成災特征。β 定義了地質災害點在因素分級中相對于均勻分布的豐度，β＞1 表示相對豐度更高，β＜1則相反。β＞1 的因素分級有（圖3、表2）：高程1 800～1 850 m、1 920～1 950 m 和1 950～2 000 m，坡度15°～25°、25°～35°和＞35°，坡向北東、東、南東和北，坡面曲率-0.75～-0.28（凹形）、-0.28～-0.15（凹形）、-0.15～-0.05(凹形)和0.05～0.15（凸形），石英砂巖巖組和砂巖、泥巖、頁巖巖組，距斷層距離0～50 m、300～500 m 和1 000～2 000 m，距主要公路距離0～50 m 和50～100 m，草地和裸地/稀疏植被區(qū)域。這些因素分級內，發(fā)育了相對于均勻分布豐度更高的地質災害，表征這些因素分級可能是研究區(qū)地質災害的主控因素。

圖3 各因素分級分區(qū)和地災點數(shù)量相關性統(tǒng)計圖Fig.3 Statistical charts of correlation between the factors and the number of geological hazard points

2 方法

2.1 權重計算

把研究區(qū)柵格單元化，利用條件概率計算證據(jù)因素圖層所有單元對地質災害發(fā)生的貢獻權重[13-15,26-27]。定義D為已發(fā)生地質災害的單元，為未發(fā)生地質災害的單元，B為證據(jù)因素區(qū)內的單元，為證據(jù)因素區(qū)外的單元。

證據(jù)因素B條件下D的條件（后驗）概率為：

式中：O(D)——證據(jù)因素B 的先驗概率，O(D)=

P(B|D)、P(B|)——在地質災害發(fā)生（D）和未發(fā)生（）時，證據(jù)因素B的條件概率，取自然對數(shù)即是證據(jù)權法中的正權重（證據(jù)因素存在區(qū)的權重值）W+。

用D和B的單元數(shù)N可表示為：

同式（3）—（6）：

根據(jù)式（11）和（12），使用ArcGIS 空間分析工具執(zhí)行權重W+和W-計算。

W+的大小表明證據(jù)因素的存在與地質災害發(fā)生之間存在正相關關系。W-表示負相關，即證據(jù)因素存在抑制誘發(fā)地質災害的作用。證據(jù)因素原始數(shù)據(jù)缺失區(qū)域的權重值取0。兩個權重之間的差異Wf=W+-W-，即綜合權重，量化證據(jù)因素和地質災害相關性大小。如果Wf為正，則證據(jù)因素對地質災害有利，如果為負，則對滑坡不利。如果Wf接近于零，則表明證據(jù)因素與地質災害的相關性不大。

2.2 因素分級（分類）的顯著性測試

在上述權重值計算及分析的基礎上，實施證據(jù)因素分類的優(yōu)選，選擇類間差異顯著的證據(jù)因素類，歸并不顯著的證據(jù)因素類。選擇近似學生化檢驗（Student-T）統(tǒng)計值進行顯著性測試[15,28]：

式中：σW+、σW-——分別是W+和W-的標準差；

Wf——綜合權重；

σWf——綜合權重標準差。

當測試值的絕對值 |Student-T|為1.96 和2.326 時，置信度達97.5%、99%，文中以 |Student-T|=2作為閾值。先將證據(jù)因素劃分為若干分級（分類），計算權重和標準差、Student-T，將 |Student-T|＜2的各分類視為顯著性低并歸為一類，保留 |Student-T|≥2的因素分類，然后重新計算歸并后各分類的權重值。

2.3 證據(jù)綜合

根據(jù)貝葉斯法則，任一單元K為地質災害的可能性，即對數(shù)后驗概率可表示為[13-15,26,27]：

式中：Bi——第i個證據(jù)因素層；

K(i)——Wi是第i個證據(jù)因素存在或不存在的權重，在第i個證據(jù)因素層存在時是+，不存在時是-。

最后計算后驗概率：

后驗概率的大小作為易發(fā)性高低的指標，值越大表示易發(fā)性越高，值越小表示易發(fā)性越低。

3 結果與討論

3.1 證據(jù)權重計算結果和相關性分析

證據(jù)權重計算結果（表2、圖4）與1.3 節(jié)可相互印證。在地形高程方面，1 800～1 850 m、1 920～1 950 m 和1 950～2 000 m 段利于地質災害發(fā)生，正權重0.555 0、1.175 8 和0.643 9。＞35°和15°～25°的山體斜坡較易于地質災害發(fā)生，正權重0.543 6 和0.378 5。坡向因素各分級權重值均不高，表明坡向對地質災害發(fā)生的驅動作用可能不太顯著。坡面曲率結果顯示，-0.75～-0.28（凹形）和-0.28～-0.15（凹形）兩個凹形坡分級段較易于地質災害發(fā)生，正權重0.569 0 和0.757 7。工程地質巖組各巖組分類的正權重值總體不高，但砂巖、泥巖、頁巖巖組的統(tǒng)計結果仍然表現(xiàn)出對地質災害發(fā)生的較有利性，其正權重0.447 4，高于排在第二位的石英砂巖巖組（正權重值為0.294 7）。距斷層距離和距主要公路距離因素統(tǒng)計結果均顯示出了較明顯的距離效應，即距斷裂或主要公路遠的地區(qū)與地質災害發(fā)生負相關，距斷裂0～50 m 和距主要公路0～50 m、50～100 m 易于地質災害發(fā)生，其正權重0.797 3、0.982 0 和0.511 1。裸地或稀疏植被地區(qū)是易于地質災害發(fā)生的區(qū)域，其正權重0.871 9。

圖4 因素證據(jù)權重計算結果圖Fig.4 Calculation results charts of factor evidence weights

續(xù)表2

3.2 模型ROC 曲線和AUC

采用接受者操作特性曲線（Receiver Operating Characteristic Curve，ROC）和ROC 曲線下與坐標軸圍成的面積（Area Under Curve，AUC）[29-32]評估模型擬合精度。模型擬合精度越好則AUC越接近1，0.7～0.9 時表示較好。文中建立的證據(jù)權法模型的AUC為80.4%，擬合精度優(yōu)異（圖5）。

圖5 模型預測性能ROC 曲線圖Fig.5 ROC curve of model prediction performance

3.3 易發(fā)性評價結果

綜合自然間斷點分級和地質災害分布，圈定了高易發(fā)區(qū)、中易發(fā)區(qū)和低易發(fā)區(qū)（表3、圖6），其中高易發(fā)區(qū)188.55 km2（占研究區(qū)總面積的49.41%），中易發(fā)區(qū)152.21 km2（占研究區(qū)總面積的39.88%），89.9%和9.1%的地災點落入高易發(fā)區(qū)和中易發(fā)區(qū)，顯示易發(fā)性分區(qū)符合已發(fā)地質災害分布，模型預測性能較好。

表3 地質災害易發(fā)性分區(qū)表Table 3 Form of geological hazard susceptibility zoning

結合地質環(huán)境因素特征分析西部高易發(fā)區(qū)（圖6藍色框范圍內、圖7）主要位于砂巖、泥巖和頁巖巖組，斷裂構造較密集，以山谷斜坡地貌為主，坡度15°～25°和＞35°較陡峭斜坡范圍成片發(fā)育且面積較廣，主要公路建于本區(qū)山谷，裸地/稀疏植被和草地連片覆蓋范圍較大。預測圈定的高易發(fā)區(qū)的這些分布特征，與上文分析得到的地質災害控制因素特征吻合，預測結果符合地質災害空間分布特征。

圖6 地質災害易發(fā)性柵格圖Fig.6 Grid map of geological hazard susceptibility

圖7 典型區(qū)因素和地質災害分布圖Fig.7 Factors and geological hazards in typical zone

4 結論

（1）“因子面積百分比A”“地災數(shù)百分比B”和“比率β”，以及各因素各分類地質災害證據(jù)權重可以定量地分析各因素與地質災害發(fā)生的相關性。

（2）圈定高易發(fā)區(qū)188.55 km2（占總面積的49.41%），中易發(fā)區(qū)152.21 km2（占總面積的39.88%），易發(fā)性分區(qū)圖具有較好的等級區(qū)分度。

（3）通過證據(jù)權法繪制的地質災害易發(fā)性圖可以有效地預測該區(qū)地質災害，模型擬合精度AUC=80.4%。89.9%和9.1%的地災點落入高和中易發(fā)區(qū)，建模結果與實際地質災害發(fā)育情況吻合度高，較好地揭示了研究區(qū)地質災害易發(fā)性特征。

（4）證據(jù)權法在研究區(qū)這類云南高原低山丘陵區(qū)有效性高，方法理論清晰，較為成熟，由數(shù)據(jù)驅動，參數(shù)定義明確，易于一線工程師推廣使用。同時，該方法權重的估計和模型預測性能受預測因子選擇、因子數(shù)據(jù)空間分辨率、因子分級影響較大，具體工作中宜對這些問題進行深入研究和統(tǒng)計分析。建議通過對因子分級進行顯著性測試實施優(yōu)選，減小對權重的高估或低估，提高模型效能。