基于證據權法的北流市南部地質災害易發(fā)性評價

磨英飛，潘宏堅，鄧蔭萬

(廣西壯族自治區(qū)地質環(huán)境監(jiān)測站，廣西 南寧 530201)

廣西北流市是典型的桂東南花崗巖丘陵區(qū)，滑坡、崩塌地質災害多發(fā)，開展地質災害易發(fā)分區(qū)，在地質災害預警和防治規(guī)劃、國土空間規(guī)劃中具有重要應用。目前，地質災害易發(fā)性評價方法有綜合指數法、證據權法、信息量法、多元統計分析法、邏輯回歸法等，綜合指數法、信息量法、多元統計分析法等在選取評價因子和權重中存在較大主觀性和人為因素，不同人員選取評價因子和賦予的權重可能會有差別，同時評價模型可靠度直接取決于地質災害原始數據的精度和發(fā)育數量，存在評價因子組合狀態(tài)比較多、樣本需求大、實際統計數量受限等缺點，具有一定的局限性。證據權法是采用貝葉斯條件概率原理、綜合各種證據層來預測某種事件發(fā)生概率的一種定量方法，在地質災害易發(fā)性評價中，充分考慮了地質災害與評價因子之間的相關關系以及評價因子之間的相互關系，評價結果較其他方法具有客觀性和可靠性。國內學者和工作人員應用證據權法[1-6]開展地質災害易發(fā)性評價，均取得了較好的效果。如胡燕等[7]采用證據權法開展巴東縣城滑坡災害易發(fā)性評價，將巴東縣城劃分為滑坡極高易發(fā)區(qū)、高易發(fā)區(qū)、中易發(fā)區(qū)、低易發(fā)區(qū)與極低易發(fā)區(qū)5類；羅可[8]應用證據權法進行大埔縣崩塌滑坡地質災害易發(fā)性評價分區(qū)，98.37%的災害點位于中、高易發(fā)區(qū)內。本文采用證據權法對廣西北流市南部石窩鎮(zhèn)、六靖鎮(zhèn)、清灣鎮(zhèn)[9]等地區(qū)進行地質災害易發(fā)性評價和分區(qū)。

1 研究區(qū)地質環(huán)境概況

1.1 氣象水文

研究區(qū)位于廣西東南的玉林市北流市南部，面積397.9 km2，多年平均降雨量約2 200 mm，降雨多集中在4—9月，占年均降雨量80%。河流發(fā)育有5條小河流。

1.2 地形地貌

研究區(qū)屬構造剝蝕丘陵地貌，分為高丘、低丘2種類型。

(1)高丘區(qū)山頂高程+300.0～+616.8 m，自然坡度一般30°～40°，局部45°～50°，溝谷切割深200～300 m，多呈“V”形，分布面積約192.6 km2。

(2)低丘區(qū)山頂高程+100.0～+300.0 m，自然坡度一般20°～30°，局部35°～40°，溝谷切割深50～150 m，多呈“U”形，分布面積205.1 km2。

植被以速生桉樹為主，覆蓋率較高，村莊多分布在低丘區(qū)的丘陵坡腳。

1.3 地層巖性

研究區(qū)地層巖性有侵入巖和變質巖(圖1)。

侵入巖侵入時期為寒武紀—志留紀、白堊紀，巖性以黑云(鉀長)二長花崗巖為主，少量黑云花崗片麻巖、片麻狀黑云二長花崗巖。

變質巖包括寒武系淺變質巖及晚元古代平政組花崗變粒巖。其中，寒武系淺變質巖巖性為千枚巖、片巖、板巖，由含礫長石石英砂巖、不等粒砂巖、粉砂巖、頁巖淺變質形成；平政組花崗變粒巖巖性為黑云變粒巖、黑云二長變粒巖、黑云斜長變粒巖及少量黑云石英片巖、含硅線堇青黑云片巖等。

花崗巖、花崗變粒巖區(qū)丘陵山坡分布第四系坡殘積粉質黏土、砂質黏土，厚2～3 m；下部為全風化層砂土，以礫砂、粗砂為主，局部含黏性土、中砂、粉砂，呈散體狀，遇水易軟化崩解，厚度10～30 m。淺變質巖區(qū)丘陵山坡第四系坡殘積黏土、粉質黏土，厚度一般在0.5～3.0 m，局部厚5.0～8.0 m。

1.4 構造

研究區(qū)構造以南南西—北北東向為主，主要發(fā)育有陸川米場—石窩斷裂、大坡外—六靖斷裂、上珍斷裂、中龍斷裂、侯山斷裂、六侖斷裂、黃田斷裂7條次級斷裂。

1.5 人類工程活動

研究區(qū)破壞地質環(huán)境條件的人類工程活動主要有城鎮(zhèn)工程建設、居民臨坡建房、公路建設形成高切坡、高填方等。集鎮(zhèn)及村屯建設挖、填方形成高5～15 m、坡度50°～80°的高陡人工邊坡，道路建設高挖低填形成高5～20 m、坡度45°～70°的高陡邊坡，多數人工邊坡無支護、排水等防護措施，強降雨作用下易引發(fā)崩塌、滑坡地質災害。經調查發(fā)現，研究區(qū)有553處高陡人工邊坡。

1.6 地質災害發(fā)育概況

研究區(qū)發(fā)育地質災害點155處。其中，崩塌30處，滑坡17處，不穩(wěn)定斜坡108處。崩塌滑坡物質以殘坡積土體和全風化花崗巖為主，花崗巖、花崗變粒巖區(qū)共發(fā)育地質災害106處，占災害點總數的68.4%，淺變質巖區(qū)共發(fā)育地質災害49處，占災害點總數的31.6%。地質災害以小型為主?；麦w厚1～5 m，崩塌體厚2～3 m。地質災害主要發(fā)育在居民屋后人工邊坡和道路兩側的人工邊坡。

2 地質災害易發(fā)性評價

2.1 證據權法簡介

證據權法是一種條件獨立假設的前提下，基于貝葉斯條件概率的定量預測方法。證據權法最初是作為醫(yī)療診斷支持的方法，20世紀80年代末，加拿大數學地質學家Bonham-Carter G F等[10-15]將該方法引入到礦產資源定量預測與評價。該方法作為一種人工智能模型，已被國內外學者廣泛用于多元信息綜合評價，近年來又被引入到地質災害的評價中。其基本原理為：假設在一特定的研究區(qū)內，已知有n種二值證據圖層，1代表因子對災害發(fā)生的證據存在，0代表不存在，通過證據權模型給出該二值化的證據因子圖層的權重，最終疊加多元圖層，實現地質災害易發(fā)性評價。證據權法的分析流程如下。

圖1 研究區(qū)地質簡圖Fig.1 Geological sketch map of the study area

2.1.1 權重計算

(1)

式中，W+為證據因子存在區(qū)的權重值；W-為證據因子不存在區(qū)的權重值，其大小表示證據因子與地質災害發(fā)生和不發(fā)生的關系密切程度。證據因子和災點正相關表示為W+>0，W-<0，負相關為W+<0，W->0，不相關時權重為0。P為不同條件發(fā)生的條件概率。

已知地質災害點的先驗概率為Po=D/T。證據因子權重由落入特定證據因子圖層的災點數和全部災點數之比與證據因子圖層面積和調查區(qū)總面積之比的比值決定。為表示證據因子對于地質災害發(fā)生和不發(fā)生的區(qū)分能力，可計算相對系數C=W+-W-，用來度量證據因子和地質災害之間的相關性大小，C值越大越有利于地質災害的出現，C值越小越不利于地質災害的發(fā)生。

2.1.2 證據綜合

在權重值計算及分析的基礎上，通過證據層的優(yōu)選，選擇權重較大、與地質災害關系密切的證據層，剔除權重較小、與地質災害關系不密切的證據層；進一步進行證據因子相對災點的條件獨立性檢驗，剔除地質災害權重相對較小而與其他證據因子相關性大的證據層。對最終篩選出的n個關于地質災害點條件獨立的證據因子，根據貝葉斯法則，研究區(qū)任一單元K發(fā)生地質災害的可能性，即對數后驗概率見式(2)。

(2)

式中，O為D的概率，O(D)=D/(T-D)；D為存在地質災害的單元網格數；Bi為第i個證據層；K(i)為在第i個證據因子層存在時為+，不存在時為-；Wi為第i個證據因子存在或不存在的權重。

最后計算后驗概率:

(3)

后驗概率值的大小表示易發(fā)性的高低，其值在0～1。后驗概率值越大，表示易發(fā)性越高；后驗概率值越小，表示易發(fā)性越低。

2.2 證據因子選取和分級

不同地域地質環(huán)境條件和地質災害發(fā)育特征不同，選取的地質災害易發(fā)評價因子[16-19]存在差異。研究區(qū)地貌類型單一，巖組較單一，植被茂盛且類型較單一，地質災害主要形成于第四系殘坡積土體和全風化花崗巖中。根據地質災害發(fā)育特征和形成機理，選擇地形地貌、地層巖性、斜坡結構類型、土體厚度、構造、人類工程活動等證據因子進行地質災害易發(fā)性評價，各證據因子細分不同的證據層。

2.2.1 地形地貌

地形地貌選取坡度、坡面曲率等證據因子，將坡度劃分為 ≤10°、(10°，15°]、(15°，20°]、(20°，25°]、(25°，35°]、(35°，45°]、(45°，55°]、>55° 共8個證據層；將坡面曲率劃分為≤-7、(-7，-3]、(-3，3]、(3，7]、>7共5個證據層。

2.2.2 地層巖性

研究區(qū)地質災害主要發(fā)育在丘陵斜坡的第四系殘坡積粉質黏土和全風化花崗巖中，溝谷及河岸分布的第四系全新統沖洪積層地質災害不發(fā)育，各地層分布面積和地質災害發(fā)育情況如圖2所示。易崩易滑地層為晚志留世扶新單元(S3F)、六楊單元(S3Ly)侵入黑云(鉀長)二長花崗巖、寒武紀古桑單元(∈G)、平政組(Pt3p)黑云變粒巖、片巖、寒武系黃洞口組第一段(∈h1)淺變質長石石英砂巖與頁巖，共劃分為11個證據層。

圖2 研究區(qū)各地層巖性分布面積及發(fā)育災害數量柱狀圖Fig.2 Histogram of lithologic distribution area and development disaster of each stratum in the study area

2.2.3 斜坡結構類型

研究區(qū)的淺變質巖由沉積巖變質形成，繼承了沉積巖的巖層層理，按斜坡坡向與巖層層理傾向間的夾角，劃分為順向坡(夾角≤30°)、同向斜交坡(30°，60°]、橫向坡(60°，120°]、反向斜交坡(120°，150°]、逆向坡(150°，180°]共5個證據層，侵入巖及平政組花崗變粒巖按順向坡考慮。

2.2.4 土層厚度

土層厚度考慮斜坡殘坡積層和全風化層的厚度，劃分為≤ 3 m、(3，6 m]、(6，9 m]、(9，12 m]、>12 m共5個證據層。

2.2.5 構造

包括與斷裂的距離和巖體節(jié)理裂隙等結構面對地質災害的相關性。

與斷裂的距離劃分為[0，200 m]、(200，500 m]、(500，1 000 m]、>1 000 m共4個證據層。結構面考慮不同傾向的結構面對地質災害發(fā)育的控制或影響，劃分為[0°，30°]、(30°，60°]、(60°，90°]、(90°，120°]、(120°，150°]、(150°，180°]、(180°，210°]、(210°，240°]、(240°，270°]、(270°，300°]、(300°，330°]、(330°，360°] 共12個證據層。研究區(qū)結構面傾向在120°～130°、170°～190°、225°～245°、265°～280°等方位區(qū)段上有比較明顯的優(yōu)勢，與120°～140°、170°～210°、225°～270°等方向的地質災害發(fā)育數量較多基本吻合(圖3)，在0°～100°、300°～360°方位區(qū)段上地質災害分布零散。

圖3 研究區(qū)結構面傾向與地質災害分布統計Fig.3 Structural inclination and distribution statistics of geological hazards in the study area

2.2.6 人類工程活動

人類工程活動考慮丘陵坡腳房屋建設強度、切坡強度，由于研究區(qū)內道路沿線未發(fā)育地質災害點，同時切坡強度也考慮了道路的人工邊坡分布，因此不考慮交通建設強度指標。

丘陵坡腳房屋建設強度可以采用房屋面積密度表示，選用最新的1∶1萬測繪信息和遙感影像資料，統計房屋面積密度，利用自然間斷點法分為4個等級，劃分為[0，23 780]、(23 780，74 312]、(74 312，157 540]、>157 540共4個證據層。

切坡強度采用高陡斜坡和地質災害的分布點密度表示，研究區(qū)調查了553處高陡斜坡和155處地質災害點，大部分地質災害點都是建房切坡引發(fā)的，高陡斜坡在強降雨等因素作用下可能失穩(wěn)發(fā)生崩塌滑坡，因此高陡斜坡和地質災害的分布點密度劃分為[0，11.41]、(11.41，28.29]、(28.29，55.67]、>55.67共4個證據層。

2.3 證據因子權重計算

坡度、坡面曲率、斜坡結構類型、土層厚度、與斷裂的距離、結構面、丘陵坡腳房屋建設強度、切坡強度等均為連續(xù)型數據，地層巖性為分類型數據。對于連續(xù)型證據因子數據，采用計算累計證據權重判斷拐點的方法進行分級。先按照較小的級距將連續(xù)性證據因子分級，然后在ArcGIS平臺下初步計算W+、W-累計權重和相對系數C，初步判斷累計權重值曲線的拐點作為分級的斷點值，將證據因子劃分為若干級，重新計算各個因子等級的權重，再繪制累計權重曲線，判斷曲線拐點和分級的合理性，并進行多次迭代，直至分類趨于合理；分類型數據可直接分級。最終各評價因子的證據層權重及后驗概率計算結果見表1—表9。

表1 坡度證據層權重及后驗概率計算結果Tab.1 Calculation results of slope evidence layer weight and posterior probability

表2 坡面曲率證據層權重及后驗概率計算結果Tab.2 Calculation results of slope curvature evidence layer weight and a posteriori probability

2.4 證據因子獨立性檢驗

證據權法的前提是各證據因子相互之間應滿足條件獨立性假設。為確保評價指標的相互獨立性，需對評價指標進行相關性分析，此次采用卡方檢驗方法進行檢驗。首先，將各因子圖層二值化，即將相對系數C<0和C=0對應的柵格賦值為0，將C>0的柵格單元賦值為1，在GIS軟件內將地質災害點與二值模式柵格圖層疊加分析；然后，用地質災害點圖層匯總二值的數量，計算期望值和卡方值(二聯表計算過程簡略)，各證據因子相互的卡方值見表10。

表3 地層巖性證據層權重及后驗概率計算結果Tab.3 Stratum lithology evidence layer weight and posterior probability calculation results

表4 斜坡結構證據層權重及后驗概率計算結果Tab.4 Calculation results of evidence layer weight and posterior probability of slope structure

表5 土層厚度證據層權重及后驗概率計算結果Tab.5 Calculation results of evidence layer weight and posterior probability of soil layer thickness

表6 與斷裂的距離證據層權重及后驗概率計算結果Tab.6 Calculation results of weight and posterior probability of evidence layer of distance from fault

根據統計學概率，當P(K2≥K)=0.001時，查表得卡方臨界值為10.828，即當卡方值小于10.828時(可信度99.9%)，兩因子之間不存在相關性，此時2組證據因子層可以同時參與易發(fā)性指數計算。

表7 結構面傾向證據層權重及后驗概率計算結果Tab.7 Calculation results of evidence layer weight and posterior probability of structural plane tendency

表8 房屋建設強度證據層權重及后驗概率計算結果Tab.8 Calculation results of evidence layer weight and posterior probability of housing construction intensity

表9 切坡強度證據層權重及后驗概率計算結果Tab.9 Calculation results of evidence layer weight and posterior probability of slope cutting strength

根據表10計算結果，地形坡度、坡面曲率、地層巖性、土層厚度、切坡強度5個證據因子相互之間的卡方值小于臨界值10.828，滿足各證據因子相互之間獨立性假設，滿足基于貝葉斯條件概率的證據權計算方法。研究區(qū)地質災害主要發(fā)育于第四系殘坡積土體和全風化花崗巖、花崗變粒巖土體中，斜坡結構類型、結構面對地質災害的控制和影響僅屬個別現象，斷裂對地質災害的影響微弱。因此，選擇地形坡度、坡面曲率、地層巖性、土層厚度、切坡強度5個證據因子組合進行易發(fā)性評價，也基本符合研究區(qū)的地質災害發(fā)育現狀。

2.5 易發(fā)性評價

研究區(qū)易發(fā)性評價最終選取地形坡度、坡面曲率、地層巖性、土層厚度、切坡強度等5個證據因子，采用證據權法計算各證據因子的權重和后驗概率，利用ArcGIS提取各證據層所在斜坡單元的后驗概率值，按斜坡單元對各證據因子的后驗概率值進行疊加計算，采用自然間斷法，將研究區(qū)劃分為地質災害極高易發(fā)區(qū)、高易發(fā)區(qū)、中易發(fā)區(qū)、低易發(fā)區(qū)4個等級(圖4)，其中極高易發(fā)區(qū)面積62.16 km2、高易發(fā)區(qū)135.22 km2、中易發(fā)區(qū)171.23 km2、低易發(fā)區(qū)29.29 km2。

表10 各證據因子相互卡方值計算結果Tab.10 Calculation results of mutual chi-square values of various evidence factors

圖4 北流市南部易發(fā)性評價分區(qū)Fig.4 Zoning map of vulnerability assessmentin the south of Beiliu City

根據易發(fā)性分區(qū)結果，極高易發(fā)區(qū)分布有69處地質災害點，高易發(fā)區(qū)分布有45處，中易發(fā)區(qū)分布有23處，低易發(fā)區(qū)分布有18處，極高、高易發(fā)區(qū)分布的災害點占73.5%。通過繪制ROC曲線(圖5)，計算AUC值為 0.718，說明證據權法在北流市南部易發(fā)性評價的可信度較高。

3 結語

采用證據權法開展北流市南部地質災害易發(fā)性評價，根據北流市南部地質災害發(fā)育和分布特征，選擇合理的證據因子，通過檢測各證據因子獨立性，選擇地形坡度、坡面曲率、地層巖性、土層厚度、切坡強度5個證據因子組合作為評價指標，劃分出地質災害極高、高、中、低易發(fā)區(qū)。經驗證，73.5 %的地質災害點落在極高、高易發(fā)區(qū)中，劃分結果具有較高的可信度。

圖5 易發(fā)性分區(qū)與地質災害數量ROC曲線Fig.5 Prone zoning and ROC curve of geological disaster quantity