基于XGBoost 模型的三峽庫區(qū)燕山鄉(xiāng)滑坡易發(fā)性評價與區(qū)劃

吳宏陽，周超,2，梁鑫，袁鵬程，余藍(lán)冰

（1.中國地質(zhì)大學(xué)（武漢）地理與信息工程學(xué)院，湖北 武漢 430078；2.三峽庫區(qū)地質(zhì)災(zāi)害野外監(jiān)測與預(yù)警示范中心，重慶 404199；3.中國地質(zhì)大學(xué)（武漢）工程學(xué)院，湖北 武漢 430074）

0 引言

滑坡作為我國主要的地質(zhì)災(zāi)害類型，具有分布地區(qū)廣、發(fā)生次數(shù)多、威脅嚴(yán)重等特點[1]，對國民經(jīng)濟(jì)建設(shè)和自然資源可持續(xù)利用造成了不可估量的破壞。由于復(fù)雜脆弱的地質(zhì)環(huán)境，三峽庫區(qū)是全國滑坡發(fā)生最多、經(jīng)濟(jì)損失最為嚴(yán)重的地區(qū)之一[2]。自2003 年蓄水以來，三峽庫區(qū)新滑坡的發(fā)生和老滑坡的復(fù)活明顯增多[3-4]。目前，庫區(qū)內(nèi)共有滑坡4,664 個，其中674 個有明顯變形特征[5]。因此開展定量化的滑坡災(zāi)害空間預(yù)測與精細(xì)化風(fēng)險評價研究十分必要。

滑坡易發(fā)性建模的本質(zhì)是研究滑坡災(zāi)害在地質(zhì)、環(huán)境及人類工程活動等因素影響下發(fā)生的空間概率[6]，準(zhǔn)確、可靠的評價結(jié)果能為風(fēng)險防控措施的制定提供可靠的科學(xué)依據(jù)。近幾十年來，國內(nèi)外學(xué)者開展了大量滑坡易發(fā)性建模研究，主要包括知識驅(qū)動模型和數(shù)據(jù)驅(qū)動模型[7-11]。隨著對地觀測技術(shù)的發(fā)展，高質(zhì)量的區(qū)域滑坡數(shù)據(jù)獲取成為可能。由于相對簡單的操作和可靠的性能，數(shù)據(jù)驅(qū)動模型逐漸在滑坡易發(fā)性評價中受到歡迎，主要可以分為數(shù)理統(tǒng)計模型[12-14]和機(jī)器學(xué)習(xí)模型[15-16]兩類。由于具有更強(qiáng)的非線性預(yù)測能力，機(jī)器學(xué)習(xí)模型在易發(fā)性建模中表現(xiàn)出更高的預(yù)測精度，常用的算法有人工神經(jīng)網(wǎng)絡(luò)[17-18]、支持向量機(jī)[19-20]、決策樹（decision tree,DT）[21-22]等。

DT 是一種經(jīng)典的樹型結(jié)構(gòu)分類算法，由于計算速度快、訓(xùn)練簡單、便于理解和解釋性好等優(yōu)點被廣泛運用，但在訓(xùn)練過程中易產(chǎn)生較復(fù)雜的模型，導(dǎo)致數(shù)據(jù)泛化能力差，出現(xiàn)過擬合情況，甚至微小的數(shù)據(jù)變化也會導(dǎo)致預(yù)測結(jié)果出現(xiàn)較大偏差。相比而言，基于DT 和Boosting 集成的梯度提升樹模型（gradient boosting decision tree,GBDT）引入隨機(jī)性，降低模型過度訓(xùn)練的可能性，能夠較好地擬合多維復(fù)雜數(shù)據(jù)，并在相對短的時間內(nèi)對海量數(shù)據(jù)得出較好的結(jié)果[23]，但模型仍存在損失函數(shù)難收斂和難以處理特征缺失樣本等問題。極致梯度提升模型（extreme gradient boosting,XGBoost）是在GBDT 基礎(chǔ)上進(jìn)行優(yōu)化得到，通過在GBDT 損失函數(shù)中加上正則項和二階導(dǎo)數(shù)來降低模型運行復(fù)雜度以及權(quán)衡模型方差，從而學(xué)習(xí)出更簡單的模型，并進(jìn)一步防止模型過擬合。在處理特征值有缺失的訓(xùn)練樣本時，XGBoost 還可以自動學(xué)習(xí)并擬合出數(shù)據(jù)的分裂方向。由于XGBoost 具備預(yù)測精度高、穩(wěn)定性好等特點，現(xiàn)已被廣泛地應(yīng)用到醫(yī)學(xué)預(yù)測、電力估計等領(lǐng)域[24-27]，而在滑坡易發(fā)性評價領(lǐng)域運用較少。

本文以三峽庫區(qū)萬州區(qū)燕山鄉(xiāng)為研究區(qū)，選取坡度、工程地質(zhì)巖組、堆積層厚度等九個指標(biāo)因子構(gòu)建易發(fā)性指標(biāo)體系，應(yīng)用信息量模型定量分析各指標(biāo)與滑坡發(fā)育關(guān)系；分別采用XGBoost、GBDT 和DT 對研究區(qū)開展易發(fā)性評價研究，并從預(yù)測精度和穩(wěn)健性方面對模型性能進(jìn)行綜合對比分析。

1 方法原理

1.1 決策樹模型

決策樹模型是一種對實例進(jìn)行分類的樹形結(jié)構(gòu)，決策樹由節(jié)點和有向邊組成，其中內(nèi)部節(jié)點代表一個特征或者一種屬性，葉節(jié)點代表類別[28]。在模型運算過程中，首先將實例從根節(jié)點開始排列，然后將屬性和特征在中間節(jié)點根據(jù)特定規(guī)則分割為兩個子集，直到在葉節(jié)點得到兩個分類結(jié)果（圖1）。其中基尼系數(shù)構(gòu)成根節(jié)點到葉子節(jié)點的分類規(guī)則序列，基尼系數(shù)值越小，代表構(gòu)建分類標(biāo)準(zhǔn)越好，最終的分類結(jié)果精度越高。

圖1 決策樹模型流程圖Fig.1 Flowchart of decision tree model

1.2 梯度提升樹模型

梯度提升樹模型是一種常用的樹形模型，采用Boosting 方法將多顆決策樹進(jìn)行關(guān)聯(lián)，新決策樹的生成是在上一棵樹損失函數(shù)的梯度下降方向，通過不斷迭代來優(yōu)化模型[29]（圖2）。

圖2 梯度提升樹模型流程圖Fig.2 Flowchart of gradient boosting decision tree model

設(shè)訓(xùn)練樣本為xi，初始損失函數(shù)F0為：

式中：xi——訓(xùn)練樣本。

利用生成的決策樹hj(xi)去擬合損失函數(shù)的梯度下降方向，使損失函數(shù)得到第j輪的最佳擬合值rj：

并利用損失函數(shù)對模型進(jìn)行更新，得到最終預(yù)測結(jié)果計算函數(shù)為：

FM為單顆決策樹預(yù)測結(jié)果，將多棵決策樹預(yù)測結(jié)果求和可得到梯度提升樹模型最終預(yù)測結(jié)果。

1.3 極致梯度提升模型

XGBoost 是一種基于GBDT 的模型，其模型結(jié)構(gòu)與GBDT 類似，都是以決策樹為基礎(chǔ)，通過不斷迭代，集成弱分類器為強(qiáng)分類器。隨著模型迭代次數(shù)的增多，預(yù)測精度也會不斷提高，計算流程如圖3 所示。與GBDT 不同的是它在損失函數(shù)中加入了二階導(dǎo)數(shù)hi和正則項 Ω(fk)來對每一輪的目標(biāo)函數(shù)進(jìn)行優(yōu)化，目標(biāo)函數(shù)值越小，則樹結(jié)構(gòu)越好[30-31]。

圖3 極致梯度提升模型流程圖Fig.3 Flowchart of extreme gradient boosting model

為了求得最小化目標(biāo)函數(shù)，分別進(jìn)行二階泰勒展開、正則化項展開來合并一次項系數(shù)和二次項系數(shù)，經(jīng)過多輪迭代后得到最終預(yù)測結(jié)果計算公式為：

f(x)——其中一棵回歸樹；

Ω(fk)——第k棵樹的正則項。

2 研究區(qū)概況與數(shù)據(jù)準(zhǔn)備

2.1 基本概況

燕山鄉(xiāng)位于三峽庫區(qū)萬州區(qū)西南部長江右岸，面積約56.93 km2。燕山鄉(xiāng)屬構(gòu)造剝蝕中淺切割丘陵河流地貌，地勢東高西低，海拔范圍在120～1 430 m。區(qū)內(nèi)地形總體向西傾向長江，多形成單傾斜坡地形，長年受雨水沖刷切割形成縱向凹槽、沖溝、溪流，匯集于長江。燕山鄉(xiāng)年平均降雨量為1 193.3 mm，日最大降雨為243 mm，最大連續(xù)降雨量為488.7 mm，夏季大雨、暴雨頻繁，極易誘發(fā)滑坡災(zāi)害。區(qū)內(nèi)共發(fā)育滑坡災(zāi)害33 處，其中小型滑坡6 處，中型滑坡22 處，大型滑坡5 處。研究區(qū)地理位置及滑坡分布如圖4 所示。

圖4 研究區(qū)位置及滑坡分布Fig.4 Location of the study area and distribution of landslides

研究區(qū)位于方斗山背斜西側(cè)，區(qū)內(nèi)地層巖性復(fù)雜，主要出露地層為侏羅系新田溝組、自流井組、珍珠沖組和三疊系巴東組、嘉陵江組、大冶組，巖性為砂巖夾泥巖、泥巖、頁巖和灰?guī)r等。區(qū)內(nèi)堆積層整體較薄，多分布于河谷斜坡、山頂侵蝕平臺之上。河谷地帶多為沖洪積亞黏土、砂土、含土卵礫石等，局部具二元結(jié)構(gòu)，而山間斜坡地帶主要為殘積、坡積、崩積等重力堆積的含碎石土，總體上，區(qū)內(nèi)堆積體結(jié)構(gòu)松散，孔隙發(fā)育，為堆積層滑坡發(fā)育提供充分條件。隨著社會經(jīng)濟(jì)的快速發(fā)展，頻繁的人類工程活動（如G69 高速公路等）對周邊自然環(huán)境造成了影響，特別是三峽庫區(qū)周期性蓄降水，造成了許多老滑坡的復(fù)活及穩(wěn)定斜坡的變形。

2.2 滑坡編錄數(shù)據(jù)

基于歷史滑坡編錄數(shù)據(jù)、影像信息和野外實地調(diào)查建立研究區(qū)滑坡災(zāi)害編錄數(shù)據(jù)庫。根據(jù)野外實際調(diào)研結(jié)果，研究區(qū)滑坡主要誘發(fā)因素有大氣降雨、人類工程活動、水庫蓄水等。其中沙榜咀滑坡在暴雨和連續(xù)強(qiáng)降雨下發(fā)生變形，大量雨水匯集，一方面加強(qiáng)了地下水入滲，使滑體本身的重量變大，另一方面當(dāng)雨水入滲到滑動面，也減弱了滑帶土的抗剪強(qiáng)度，增加了滑體的下滑力。桐子林滑坡滑體物質(zhì)為大量后期人工生活及建筑堆積物和崩積含碎塊石粉質(zhì)黏土；五尺壩滑坡為區(qū)內(nèi)典型的受庫水位影響的庫岸滑坡，其中前緣以滑塌變形為主，受到庫水位浸泡和側(cè)向沖刷是誘發(fā)滑坡前緣變形的主要原因。研究區(qū)滑坡從發(fā)育平面形態(tài)特征來看，共發(fā)育舌形滑坡16 個，箕形滑坡7 個，扇形滑坡5 個，橫長形滑坡3 個，不規(guī)則長條形滑坡2 個。從發(fā)育剖面形態(tài)特征來看，直線形所占比例最大，其比例為91%，圓弧形較少分布，其比例為9%。整體來看，研究區(qū)滑坡在東西方向數(shù)量分布存在明顯差異，且東西方向上堆積層厚度、坡度、工程地質(zhì)巖組等指標(biāo)存在明顯差異，根據(jù)主要誘發(fā)因可分為受水庫波動影響較大的緩傾角堆積層滑坡和受降雨影響較大的陡傾角堆積層滑坡（圖5）。

圖5 典型滑坡全貌圖Fig.5 Overview of typical landslide

2.3 易發(fā)性評價指標(biāo)體系

影響滑坡發(fā)育的指標(biāo)因子主要包括有地形地貌、地質(zhì)條件和人類工程活動等[32]。在前人對萬州區(qū)進(jìn)行滑坡易發(fā)性評價基礎(chǔ)上[33-35]，考慮研究區(qū)地質(zhì)背景、滑坡形成條件及其發(fā)育特征，并根據(jù)野外實際調(diào)研中，研究區(qū)東西方向地理條件差異鮮明情況，選取土地利用、斜坡結(jié)構(gòu)、巖組、坡度、距長江距離、堆積層厚度、植被歸一化指數(shù)、斜坡形態(tài)、匯水面積九個指標(biāo)因子構(gòu)建研究區(qū)易發(fā)性評價指標(biāo)體系。為了判別指標(biāo)之間相關(guān)性密切程度，利用ArcGIS 中波段集統(tǒng)計工具檢驗各指標(biāo)因子之間的相關(guān)性，結(jié)果顯示各指標(biāo)因子之間相關(guān)性皆小于0.4，表明指標(biāo)因子之間呈弱相關(guān)性或不相關(guān)（圖6），可用作研究區(qū)滑坡易發(fā)性評價建模，各因子指標(biāo)分級情況如圖7 所示。

圖6 指標(biāo)相關(guān)性Fig.6 The correlation plot of Indicator factors

圖7 研究區(qū)易發(fā)性評價指標(biāo)圖Fig.7 Indicator plot for landslide susceptibility assessment in the study area

3 滑坡易發(fā)性建模

3.1 模型建立

將指標(biāo)因子作為模型輸入數(shù)據(jù)，滑坡預(yù)測相對概率作為模型輸出進(jìn)行易發(fā)性建模。隨機(jī)選擇70%的滑坡和相同數(shù)量的非滑坡數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)，剩余30%滑坡數(shù)據(jù)則作為驗證數(shù)據(jù)[36]。

為了保證模型在對比時達(dá)到一個較優(yōu)的狀態(tài)，需要對各個模型進(jìn)行調(diào)參。最大葉子節(jié)點數(shù)是決策樹模型主要影響參數(shù)，通過試算法確定葉節(jié)點個數(shù)。如圖8（a）所示，DT 模型的預(yù)測精度隨著最大葉子節(jié)點數(shù)的增大呈現(xiàn)出上升趨勢，當(dāng)葉子節(jié)點數(shù)目達(dá)到160 時，預(yù)測精度達(dá)到峰值，呈穩(wěn)定狀態(tài)，因此設(shè)置決策樹模型最大葉子節(jié)點數(shù)目為160。

圖8 參數(shù)與預(yù)測精度關(guān)系曲線Fig.8 Relationship curve between parameters and prediction accuracy

GBDT 和XGBoost 為集成樹模型，主要通過創(chuàng)建新決策樹，由多顆決策樹的預(yù)測值得到最終結(jié)果，而影響集成樹模型精度的主要參數(shù)為決策樹數(shù)量。一般來說，模型決策樹數(shù)量較小時，容易欠擬合，而決策樹數(shù)量過多，計算量快速增加，且決策樹數(shù)量增加到一定時，模型逐漸趨于穩(wěn)定，預(yù)測精度不會隨決策樹數(shù)量增加出現(xiàn)大幅度變動[37]。同樣采用試算法來確定適合模型的決策樹數(shù)量，GBDT 和XGBoost 模型中決策樹數(shù)量和預(yù)測精度的關(guān)系曲線分布如圖8（b）和圖8（c）所示?？梢钥闯鯣BDT 和XGBoost 在決策樹數(shù)量分別大于145 和大于75 預(yù)測精度無明顯提升，且趨于穩(wěn)定，因此設(shè)置GBDT和XGBoost 決策樹數(shù)量分別為145 和75，此時最大決策樹深度為7，能夠最多生成的葉子結(jié)點數(shù)為49。此外，設(shè)置XGBoost 學(xué)習(xí)率、L1 和L2 正則項權(quán)重分別為0.1、0 和1。

將評價指標(biāo)輸入到訓(xùn)練好的DT、GBDT 和XGBoost模型中，計算研究區(qū)滑坡易發(fā)性相對概率值，并分為四個易發(fā)性等級區(qū)間，分別為低（77%）、中（7%）、高（7%）和極高（9%）。結(jié)果如圖9 所示。

3.2 結(jié)果對比與分析

3.2.1 滑坡空間發(fā)育規(guī)律分析

信息量模型是一種二變量統(tǒng)計方法，統(tǒng)計結(jié)果能夠很好地表征影響因素對滑坡空間發(fā)育的影響作用與程度。信息量值為正，說明指標(biāo)對滑坡發(fā)育有促進(jìn)作用；信息量值為負(fù)，則說明有抑制作用，且絕對值越大表明作用越強(qiáng)，各評價指標(biāo)的信息量計算結(jié)果如表1 所示。影響區(qū)內(nèi)滑坡發(fā)育的主要指標(biāo)有距長江距離（0～400 m）、堆積層厚度（＞2.4 m）和工程地質(zhì)巖組（砂巖夾泥巖和砂巖），信息量值分別為2.96、2.54 和1.58。受三峽水庫蓄降水的影響，研究區(qū)長江沿岸滑坡面積約占總滑坡面積的53%，并且長江支流附近水系發(fā)育，受各類水流的長期侵蝕和沖刷，斜坡整體穩(wěn)定性低；研究區(qū)第四系堆積層厚度大于2.4 m 的區(qū)域主要在沖溝和長江沿岸，此類斜坡土體結(jié)構(gòu)松散，吸水能力較強(qiáng)且持水能力差，在降雨誘發(fā)下快速形成暫時的飽水帶，易失穩(wěn)滑動；區(qū)內(nèi)泥巖風(fēng)化程度高，結(jié)構(gòu)破碎，在降雨等不利條件下極易失穩(wěn)滑動。

表1 各因素狀態(tài)信息量表Table 1 The weighted information values of each factor state

3.2.2 模型性能分析

滑坡比率為該易發(fā)性等級滑坡柵格點數(shù)量占總滑坡柵格點數(shù)量之比與分級柵格數(shù)占研究區(qū)柵格數(shù)之比的比值。統(tǒng)計各模型滑坡比率，各模型滑坡比率皆是從低易發(fā)性到極高易發(fā)性依次增大，且在極高易發(fā)區(qū)滑坡比率最高。XGBoost 極高易發(fā)性等級的滑坡比率分別10.19，高于為DT 和GBDT 的8.07、8.59。準(zhǔn)確率為分類正確的滑坡樣本個數(shù)占總滑坡樣本個數(shù)的比例[38]，通過計算得到各模型準(zhǔn)確率分別為92.63%、93.13%和94.32%，表明XGBoost 預(yù)測結(jié)果準(zhǔn)確性最高，優(yōu)于DT 和GBDT（圖10）。繪制全區(qū)受試者工作特征曲線（receiver operating characteristic curve，ROC）[39]（圖11），可以看出XGBoost 全區(qū)預(yù)測精度為0.973，優(yōu)于DT 和GBDT。

圖10 各易發(fā)區(qū)災(zāi)害點比例Fig.10 Proportion of disaster points in different susceptibility zones

圖11 模型 ROC 曲線圖Fig.11 ROC curves of the different models

為進(jìn)一步探究模型預(yù)測的穩(wěn)健性，隨機(jī)生成100 組試驗數(shù)據(jù)進(jìn)行滑坡易發(fā)性預(yù)測，得到預(yù)測精度和抽樣次數(shù)關(guān)系曲線、標(biāo)準(zhǔn)差和變異系數(shù)如圖12 和表2 所示。

表2 標(biāo)準(zhǔn)差和變異系數(shù)Table 2 Standard deviation and coefficient of variation

圖12 抽樣次數(shù)與預(yù)測精度關(guān)系曲線Fig.12 The correlation curve between sampling times and prediction accuracy

在多組抽樣數(shù)據(jù)下，XGBoost 擁有最好的預(yù)測精度，且預(yù)測精度上下起伏不大，穩(wěn)定性較好，預(yù)測精度平均值為97.28%，優(yōu)于GBDT（95.61%）和DT（90.30%），其標(biāo)準(zhǔn)差為0.010，小于GBDT（0.062）和DT（0.734），變異系數(shù)為0.010，小于GBDT（0.065）和DT（0.813）。

4 結(jié)論

以燕山鄉(xiāng)為例，選取坡度、工程地質(zhì)巖組、堆積層厚度等九個影響因子構(gòu)建易發(fā)性評價指標(biāo)體系，應(yīng)用信息量模型分析滑坡發(fā)育與指標(biāo)之間的關(guān)系。研究結(jié)果發(fā)現(xiàn)：距長江距離（0～400 m）、工程地質(zhì)巖組（砂巖夾泥巖、砂巖和泥巖夾砂巖、泥巖）和堆積層厚度（＞2.4 m）是研究區(qū)滑坡空間發(fā)育的主要控制因素；燕山鄉(xiāng)東部出露抗風(fēng)化能力較強(qiáng)的灰?guī)r且基本無第四系堆積物，不具備發(fā)育滑坡的條件，研究區(qū)西部長期受到庫水位壓力變化，斜坡失穩(wěn)概率極大。因此研究區(qū)東部與西部地區(qū)之間滑坡易發(fā)性評價結(jié)果有明顯的差異，與野外實際勘察情況一致。

通過對100 組訓(xùn)練/驗證數(shù)據(jù)集開展易發(fā)性評價探究模型預(yù)測準(zhǔn)確性、穩(wěn)定性。結(jié)果表明XGBoost 模型的標(biāo)準(zhǔn)差和變異系數(shù)均為0.01，優(yōu)于GBDT 和DT，說明該模型在多次重復(fù)中具有較好的穩(wěn)健性。對評價結(jié)果進(jìn)行驗證得到準(zhǔn)確率和預(yù)測精度分別為94.3%和97.3%，優(yōu)于GBDT 和DT。DT 是一個訓(xùn)練簡單、可理解性好的模型，但模型易過擬合且預(yù)測精度不高，通過模型集成的方法可以解決DT 易過擬合的不足；將DT 和Boosting 集成的GBDT 能夠較好地擬合多維復(fù)雜數(shù)據(jù)并降低模型過擬合可能。為解決損失函數(shù)難收斂問題，在GBDT 基礎(chǔ)上添加正則項和二階導(dǎo)數(shù)的XGBoost 擁有模型穩(wěn)定性和預(yù)測準(zhǔn)確性。XGBoost 是一種優(yōu)秀的滑坡易發(fā)性預(yù)測模型，具有較高的預(yù)測精度，并能為后續(xù)滑坡風(fēng)險評價和分析提供技術(shù)支撐。