基于確定性系數(shù)與Logistic回歸模型的地質(zhì)災(zāi)害易發(fā)性評價

鐘畝鋒

（廣西地礦建設(shè)集團有限公司，廣西 南寧 530023）

中國是世界上地質(zhì)災(zāi)害發(fā)育最齊全的國家之一，發(fā)生的地質(zhì)災(zāi)害以崩塌、滑坡和泥石流為主。據(jù)不完全統(tǒng)計，全國共發(fā)育較大型崩塌3000多處、滑坡2000多處，中小規(guī)模的崩塌、滑坡、泥石流多達(dá)數(shù)十萬處。根據(jù)自然資源部最新統(tǒng)計數(shù)據(jù)，2018年全國共發(fā)生地質(zhì)災(zāi)害2966起，造成105人死亡、7人失蹤、73人受傷，直接經(jīng)濟損失14.7億元。地質(zhì)災(zāi)害已經(jīng)成為造成人民生命財產(chǎn)損失的最重要因素，開展地質(zhì)災(zāi)害易發(fā)性評價研究愈發(fā)迫切和重要。

地質(zhì)災(zāi)害易發(fā)性評價是根據(jù)區(qū)域的地形地貌、地質(zhì)背景、水文條件、人類活動、土地覆蓋/利用等因素進(jìn)行分析，對區(qū)域地質(zhì)災(zāi)害發(fā)生的可能性進(jìn)行評估和劃分，為地質(zhì)災(zāi)害預(yù)測、防治提供技術(shù)支持。目前，國內(nèi)外學(xué)者常用的地質(zhì)災(zāi)害易發(fā)性評價模型主要有：模糊綜合評判法、層次分析法、信息量法、確定性系數(shù)法、人工神經(jīng)網(wǎng)絡(luò)等[1-5]。Schicker等[6]采用多元統(tǒng)計分析方法對區(qū)域災(zāi)害進(jìn)行易發(fā)性分析并制圖；Lee等[7]利用RS技術(shù)、GIS平臺及其空間分析工具對韓國龍仁地區(qū)進(jìn)行滑坡易發(fā)性評價研究；Floris M等[8]利用IFSAR數(shù)據(jù)進(jìn)行滑坡易發(fā)性研究。國內(nèi)學(xué)者彭珂等[9]利用信息量模型對贛州市地質(zhì)災(zāi)害進(jìn)行易發(fā)性分析；劉艷芳等[10]采用確定性系數(shù)方法研究秭歸縣滑坡易發(fā)情況；許沖等[11]利用GIS與確定性系數(shù)分析方法對汶川地震滑坡進(jìn)行易發(fā)性評價。單一的評價模型不能客觀確定評價因子權(quán)重，難以定量、準(zhǔn)確地進(jìn)行區(qū)域地質(zhì)災(zāi)害易發(fā)性評價。

本文以廣西馬山縣為例，選擇統(tǒng)計學(xué)方法來進(jìn)行地質(zhì)災(zāi)害易發(fā)性評價，主要有三個步驟：①選取評價因子并分級，通過已知災(zāi)害點計算地質(zhì)災(zāi)害評價因子分級指標(biāo)對應(yīng)的確定性系數(shù)值。②利用Logistic回歸模型分析災(zāi)害點和非災(zāi)害點樣本，計算各因子回歸系數(shù)。③將各評價因子的確定性系數(shù)值和相關(guān)系數(shù)帶入到Logistic回歸模型，計算研究區(qū)所有單元格地質(zhì)災(zāi)害發(fā)生的概率，以此來進(jìn)行地質(zhì)災(zāi)害易發(fā)性評價。

1 研究區(qū)概況

馬山縣位于廣西壯族自治區(qū)中部，居紅水河中段南岸、大明山西北麓。地理范圍在東經(jīng)107°41′～108°29′，北緯23°24′～24°2′之間，東接忻城、上林縣，西接大化瑤族自治縣、平果縣，南與武鳴縣為鄰，北隔紅水河與都安瑤族自治縣相望，總面積2345 km2，如圖1所示。馬山縣屬于南亞熱帶濕潤氣候，夏季降雨量充沛，冬季較少。境內(nèi)的河流主要屬于紅河水系與武鳴河水系。馬山縣地貌包括巖溶地貌和非巖溶地貌二大類，巖溶地貌分布面積占全區(qū)的72.1%，非巖溶地貌分布面積占全區(qū)的27.9%。馬山縣境內(nèi)出露巖性以沉積巖為主，地層巖性由下到上依次為寒武系、泥盆系、石炭系、二疊系、三疊系和第四系。地質(zhì)構(gòu)造屬于廣西山字型構(gòu)造前弧西翼的一部分，主要構(gòu)造線呈西北——東南方向。境內(nèi)已查明的地質(zhì)災(zāi)害有崩塌、滑坡、危巖、地面塌陷和地裂縫，主要以崩塌和危巖為主，地質(zhì)災(zāi)害規(guī)模多為小型、中型。

圖1 馬山縣位置圖

2 數(shù)據(jù)源

本研究所用到的樣本分為地質(zhì)災(zāi)害點數(shù)據(jù)和非地質(zhì)災(zāi)害點數(shù)據(jù)。地質(zhì)災(zāi)害點數(shù)據(jù)包括崩塌81處、滑坡3處、危巖28處、地面塌陷2處、地裂縫1處，共計115個地質(zhì)災(zāi)害點。非地質(zhì)災(zāi)害點數(shù)據(jù)是根據(jù)研究需求，利用ArcGIS隨機生成的非地質(zhì)災(zāi)害樣本115個。除了樣本數(shù)據(jù)以外，還包括地質(zhì)災(zāi)害易發(fā)性評價因子的基礎(chǔ)數(shù)據(jù)：數(shù)字高程模型（30m分辨率）、Landsat-8 OLI影像、馬山縣1：5萬地質(zhì)災(zāi)害詳細(xì)調(diào)查工作部署圖、馬山縣1：10萬水文地質(zhì)圖、土地覆蓋/利用和路網(wǎng)數(shù)據(jù)。其中，土地覆蓋/利用是利用ArcGIS支持向量機監(jiān)督分類得到的分類數(shù)據(jù)；路網(wǎng)數(shù)據(jù)是在BIGEMAP地圖下載器下載獲取的。

3 評價基本方法

3.1 Logistic回歸模型

Logistic回歸模型(Logistic Regression Model)是一種因變量滿足二項分布的多元回歸模型，通過建立多個自變量與一個因變量之間的多元回歸關(guān)系，預(yù)測某一事件發(fā)生的概率。在對地質(zhì)災(zāi)害進(jìn)行易發(fā)性評價時，選取評價因子作為自變量，將災(zāi)害是否發(fā)生作為因變量，因變量取值為1和0，分別代表災(zāi)害發(fā)生與不發(fā)生。設(shè)為災(zāi)害發(fā)生的概率，為災(zāi)害不發(fā)生的概率。若用來擬合，的取值可能不在[0，1]范圍內(nèi)，因此采用odds函數(shù)對進(jìn)行約束，odds函數(shù)為事件發(fā)生與事件不發(fā)生的比值，即，并取其自然對數(shù)作為因變量，建立線性回歸方程：

利用公式（1）和（2）便可對災(zāi)害發(fā)生概率進(jìn)行預(yù)測[12]。

3.2 確定性系數(shù)模型

確定性系數(shù)（Certainty Factor，CF）本質(zhì)上是一個概率函數(shù)，由Shortliffe[13]提出，Heckrman[14]進(jìn)行改進(jìn)，為了把評價因子量化到以0為中心，變化區(qū)間在[-1，1]范圍內(nèi)，采用確定性系數(shù)對各評價災(zāi)因子進(jìn)行量化。其計算公式如下：

式中：Pa為地質(zhì)災(zāi)害在評價因子分類a中發(fā)生的條件概率，此處表示為評價因子分類a中發(fā)生的地質(zhì)災(zāi)害個數(shù)（或面積）與評價因子分類a面積的比值；Ps為地質(zhì)災(zāi)害在整個數(shù)據(jù)中發(fā)生的先驗概率，此處表示為整個研究區(qū)內(nèi)已發(fā)生的地質(zhì)災(zāi)害個數(shù)（面積）與研究區(qū)域總面積的比值[15]。

公式（3）中CF的取值在[-1，1]之間，當(dāng)CF取正值時，表示區(qū)域容易發(fā)生地質(zhì)災(zāi)害；當(dāng)CF取值為負(fù)值時，表示區(qū)域不易發(fā)生地質(zhì)災(zāi)害。

3.3 確定性系數(shù)與邏輯回歸模型結(jié)合

圖2 評價方法流程圖

用災(zāi)害點樣本計算出每一個評價因子指標(biāo)的確定性系數(shù)值，提取災(zāi)害點樣本與非災(zāi)害點樣本對應(yīng)的每一個評價因子的確定性系數(shù)值，建立回歸方程，計算邏輯回歸常量α和每個評價因子的回歸系數(shù)，計算出區(qū)域所有單元格的概率值，以此為依據(jù)進(jìn)行地質(zhì)災(zāi)害易發(fā)性評價。具體流程如圖2所示。

4 評價因子選擇與分級

馬山縣的地貌條件主要為巖溶地貌，致使馬山縣容易發(fā)生崩塌和形成崩塌隱患點。因此，馬山縣的地質(zhì)災(zāi)害更多的是側(cè)重研究崩塌，本文在進(jìn)行馬山縣地質(zhì)災(zāi)害易發(fā)性評價時，將不同的地質(zhì)災(zāi)害類型合并統(tǒng)計分析。

根據(jù)收集的區(qū)域統(tǒng)計資料和實地現(xiàn)場踏勘，初步確定的馬山縣地質(zhì)災(zāi)害評價因子有：地形地貌、地質(zhì)背景、水文條件、人類活動、土地覆蓋/利用等五類。具體細(xì)分，地形地貌分為高程、坡度、坡向，地質(zhì)背景分為地層巖性、距斷層距離，水文條件為地下水類型，人類活動用距道路距離來表示，以及土地覆蓋/利用共8個因子，每個因子的分級指標(biāo)如表1所示。

表1 地質(zhì)災(zāi)害易發(fā)性評價因子分級指標(biāo)和CF值

距道路距離（m） 4 0～300 0.17528 300～600 0.31997 600～900 -0.33452＞900 -0.44896土地覆蓋/利用 5水域 -1裸地 -1建設(shè)用地 -0.53957耕地 0.13601林地 0.01893

高程與地質(zhì)災(zāi)害的發(fā)育無直接關(guān)系，但高程與人類活動和土地覆蓋/利用密切相關(guān)，不同高程植被的發(fā)育類型不同，人類活動的強度也不一樣，形成的土地覆蓋/利用類型也完全不一樣。坡度與產(chǎn)生崩塌、滑坡的勢能、下滑力有關(guān)，也會確定災(zāi)害的影響范圍。坡向不同，山坡受到的光照條件也不一樣，南面朝向的為陽坡，北面朝向的為陰坡。光照條件會影響到植被發(fā)育、坡面水分蒸發(fā)量、風(fēng)化程度等。斜坡巖土體的巖性影響著巖石強度、風(fēng)化程度、節(jié)理裂隙發(fā)育、地下水類型等，是地質(zhì)災(zāi)害發(fā)育的內(nèi)在控制因素。地質(zhì)構(gòu)造活動會造成區(qū)域巖體破碎，節(jié)理裂隙發(fā)育，為加速巖體風(fēng)化、降低巖土體強度提供了內(nèi)地質(zhì)營力。地下水沿裂隙面滲流，使巖土體粘合力降低，巖土體強度下降。地表河流對坡腳沖刷、侵蝕易使邊坡失穩(wěn)，水位上升會浸濕坡腳，坡岸巖土體穩(wěn)定性降低。道路的分布反映人類活動頻繁程度，道路密集說明該區(qū)域人類活動強烈。人類活動：修筑道路、人工切坡、工程爆破等都會增加區(qū)域發(fā)生地質(zhì)災(zāi)害的可能性。土地覆蓋/土地利用在一定程度可以反映區(qū)域植被類型、人類活動的情況。

5 地質(zhì)災(zāi)害易發(fā)性評價

5.1 評價因子分析

通過收集準(zhǔn)備的資料，將所有的因子圖層矢量用統(tǒng)一的矢量范圍進(jìn)行裁剪，然后全部轉(zhuǎn)化為柵格數(shù)據(jù)。同時，設(shè)置每個因子圖層的柵格大小為30m×30m。根據(jù)115個災(zāi)害點在各個評價因子分級指標(biāo)中的分布，結(jié)合公式（3）計算出每一個評價因子分級指標(biāo)的確定性系數(shù)值，確定性系數(shù)值如表1所示。利用ArcGIS的多值提取到點，就可以得到每一個災(zāi)害點對應(yīng)的評價因子分級指標(biāo)的確定性系數(shù)值。

為了確保地質(zhì)災(zāi)害易發(fā)性分析的合理性，要使因子之間保持相互獨立。本文利用相關(guān)分析對8個因子進(jìn)行獨立性檢驗，各因子之間的相關(guān)系數(shù)如表2所示。從表2可以看出，高程與坡度、地層巖性、土地覆蓋/利用的相關(guān)系數(shù)大于0.3，其余因子之間的相關(guān)性基本都小于0.3或只與某個因子相關(guān)性略偏高，綜合考慮后去除高程因子，保留其余7個因子。

利用ArcGIS在研究區(qū)范圍內(nèi)隨機生成115個非災(zāi)害點，與已知的115個災(zāi)害點一起作為地質(zhì)災(zāi)害易發(fā)性評價的統(tǒng)計樣本。將各個影響因子的確定性系數(shù)值作為自變量，是否發(fā)生地質(zhì)災(zāi)害作為因變量（取0代表地質(zhì)災(zāi)害不發(fā)生，1代表地質(zhì)災(zāi)害發(fā)生）。把330個樣本點導(dǎo)入到SPSS中進(jìn)行二項邏輯回歸分析，邏輯回歸分析的分類表如表3所示，觀測值與預(yù)測值一致的可能性為64.8%，說明了方程回歸的可靠性為64.8%?；貧w系數(shù)（B）是在回歸方程中表示自變量對因變量影響大小的參數(shù)，自變量回歸系數(shù)的絕對值越大，自變量對因變量的影響也越大，回歸分析結(jié)果如表4所示。因此，可以用評價因子回歸系數(shù)（B）的絕對值來分析評價因子的易發(fā)性，從表4分析可以得出，斷層、地下水、地層巖性、降雨對區(qū)域地質(zhì)災(zāi)害發(fā)生的作用較大，其中斷層的控制作用最顯著。

表2 評價因子相關(guān)系數(shù)表

表3 分類表

5.2 模型分析

5.2.1 合理性分析

式中：x1為坡度分區(qū)的CF值；x2為坡向分區(qū)的CF值；x3為地層巖性分區(qū)的CF值；x4為距斷層距離分區(qū)的CF值；x5為地下水分區(qū)的CF值；x6為距道路距離分區(qū)的CF值；x7為土地覆蓋/利用分區(qū)的CF值。由公式（2）計算出所有單元格發(fā)生地質(zhì)災(zāi)害的概率P，利用自然斷點法將研究區(qū)分為了5級：極高易發(fā)區(qū)、高易發(fā)區(qū)、中易發(fā)區(qū)、低易發(fā)區(qū)和穩(wěn)定區(qū)。結(jié)果如圖3所示，圖中顏色由綠到紅易發(fā)程度逐漸增加。

表4 邏輯回歸分析結(jié)果表

表5 地質(zhì)災(zāi)害易發(fā)性分區(qū)統(tǒng)計表

經(jīng)過計算得出：極高易發(fā)區(qū)占比10.48%，高易發(fā)區(qū)占比24.26%，中易發(fā)區(qū)占比27.84%，低易發(fā)區(qū)占比17.46%，穩(wěn)定區(qū)占比19.96%。從地質(zhì)災(zāi)害易發(fā)區(qū)空間分布上看，馬山縣易發(fā)生災(zāi)害的地區(qū)主要集中在東部和西部，不易發(fā)生災(zāi)害的地區(qū)主要集中在中部。將災(zāi)害點與易發(fā)性結(jié)果圖疊加分析，獲得不同程度的易發(fā)區(qū)災(zāi)害點分布情況，如表5所示。穩(wěn)定區(qū)的災(zāi)害點占比最小，災(zāi)害點密度也遠(yuǎn)小于其他易發(fā)性分區(qū)。極高易發(fā)區(qū)災(zāi)害點雖然占比不大，但災(zāi)害點密度較大，災(zāi)害點密度最大的為高易發(fā)區(qū)。研究區(qū)易發(fā)性評價的結(jié)果基本符合災(zāi)害點分布的實際情況。

圖3 地質(zhì)災(zāi)害易發(fā)性分區(qū)

5.2.2 精度分析

易發(fā)性評價結(jié)果的有效性可以用受試者工作特征曲線（ROC）進(jìn)行定性分析，AUC是曲線下方的面積，可以用來度量分類模型的好壞，AUC越接近1，模型的模擬值越接近于樣本值。本文用ROC曲線和AUC值對耦合模型進(jìn)行檢驗，ROC曲線的橫軸表示易發(fā)性面積百分比累積量，縱軸表示災(zāi)害點面積百分比累積量，結(jié)果如圖4所示。AUC的值為0.688，表明該方法能較為有效地對馬山縣地質(zhì)災(zāi)害進(jìn)行評價。

圖4 ROC曲線和AUC值

6 結(jié)論

本文以廣西壯族自治區(qū)南寧市馬山縣為例，計算高程、坡度、坡向、地層巖性等地質(zhì)災(zāi)害影響因子的確定性系數(shù)，并結(jié)合邏輯回歸模型計算單元格地質(zhì)災(zāi)害發(fā)生的概率，在此基礎(chǔ)上進(jìn)行馬山縣地質(zhì)災(zāi)害易發(fā)性評價，主要結(jié)論有：

（1）從計算得到的回歸系數(shù)可以確定，馬山縣地質(zhì)災(zāi)害受地下水、坡度、斷層和巖性的影響比較顯著。研究區(qū)地質(zhì)災(zāi)害易發(fā)性評價結(jié)果表明，馬山縣地質(zhì)災(zāi)害主要分布在碳酸鹽巖裂隙溶洞水區(qū)域，離斷層距離小于400m區(qū)域，二疊系、泥盆系碳酸鹽巖區(qū)域。馬山縣地質(zhì)災(zāi)害防治工作應(yīng)該重點關(guān)注有斷層經(jīng)過的碳酸鹽巖區(qū)域，這些區(qū)域是地質(zhì)災(zāi)害詳查的重點調(diào)查對象。

（2）基于耦合模型的易發(fā)性評價方法可以準(zhǔn)確客觀地計算出評價因子與發(fā)生地質(zhì)災(zāi)害之間的確定性關(guān)系和各個因子的權(quán)重。利用GIS平臺可以很好的解決不同數(shù)據(jù)轉(zhuǎn)換、疊加分析的問題。精度分析中ROC曲線和AUC值用來檢驗耦合模型易發(fā)性評價的有效性，檢驗結(jié)果表明該方法可以較為準(zhǔn)確的對馬山縣地質(zhì)災(zāi)害進(jìn)行評價。

本文選取的評價因子主要為一些符合區(qū)域災(zāi)害特點且容易獲取的因子，還有一些其它方面重要的影響因子未被納入其中分析，如降雨、地震烈度等。對研究區(qū)進(jìn)行單元格劃分時，還需要結(jié)合區(qū)域地質(zhì)災(zāi)害的特點，確定更加合適的單元格大小。上述兩個問題還需收集更多的數(shù)據(jù)，進(jìn)行更加深入的研究。