周天游,劉 暢,薛 鵬,楊 豹,舒建冬

(1. 中國地質(zhì)科學(xué)院探礦工藝研究所,四川 成都 610031; 2. 萬融建工有限公司,四川 成都 610074;3. 中冶成都勘察研究總院有限公司,四川 成都 610023; 4. 四川省華地建設(shè)工程有限責(zé)任公司,四川 成都 610081)

0 引言

滑坡易發(fā)性評價是指在對區(qū)域滑坡調(diào)查的基礎(chǔ)上,進(jìn)行區(qū)域滑坡易發(fā)條件和空間發(fā)生概率的預(yù)測[1],主要包括滑坡編錄與基礎(chǔ)環(huán)境因子獲取、訓(xùn)練與測試集的確定、易發(fā)性評價模型的選擇、模型參數(shù)的確定、易發(fā)性評價結(jié)果分析等[2]?，F(xiàn)階段,滑坡易發(fā)性評價方法可被劃分為啟發(fā)式、基于力學(xué)模型的確定性方法、數(shù)據(jù)驅(qū)動模型(包含數(shù)理統(tǒng)計和機(jī)器學(xué)習(xí)模型)以及集成模型等[3]。其中,常見的機(jī)器學(xué)習(xí)模型包括邏輯回歸、神經(jīng)網(wǎng)絡(luò)、多層感知器、模糊數(shù)學(xué)、支持向量機(jī)、決策樹、隨機(jī)森林以及深度學(xué)習(xí)算法等。研究表明,機(jī)器學(xué)習(xí)相比常規(guī)數(shù)理統(tǒng)計模型,具有樣本數(shù)據(jù)少、泛化能力強(qiáng)、建模過程簡單高效和精度高等優(yōu)點[2],因此,越來越多的學(xué)者在區(qū)域面積較大的滑坡易發(fā)性研究中視其為更好的選擇。

不同機(jī)器學(xué)習(xí)模型的易發(fā)性建模過程、預(yù)測精度甚至易發(fā)性指數(shù)分布特征結(jié)果往往差異較大,難以確定哪種模型效果最優(yōu),例如黃發(fā)明等[5]采用邏輯回歸、多層感知器、支持向量機(jī)和C5.0決策樹等方法,分別構(gòu)建了基于原始滑坡因子和改進(jìn)因子的機(jī)器學(xué)習(xí)模型以預(yù)測滑坡易發(fā)性,結(jié)果表明決策樹模型預(yù)測性能最好;李文斌等[6]以江西省瑞金市為例,通過5種聯(lián)結(jié)方法耦合4種機(jī)器學(xué)習(xí)模型,開展了20種模型下的滑坡易發(fā)性預(yù)測研究和易發(fā)性指數(shù)的不確定性分析;謝明娟等[7]開展了機(jī)器學(xué)習(xí)模型預(yù)測滑坡易發(fā)性的對比分析,認(rèn)為支持向量機(jī)模型具有預(yù)測結(jié)果較為穩(wěn)定和認(rèn)可度較高等優(yōu)點。值得說明的是,不同區(qū)域內(nèi)的地質(zhì)條件和環(huán)境背景差異較大,模型性能受到由特定研究區(qū)域內(nèi)生成的建模數(shù)據(jù)的影響,因此,所得到的滑坡易發(fā)性分區(qū)結(jié)果的準(zhǔn)確性及合理性也有所差異。

在構(gòu)造活躍的西南山區(qū),地震誘發(fā)滑坡是重要的次生災(zāi)害類型之一,地震波使得斜坡存在“裂”而未“滑”、“松”而未“動”的震裂山體,在強(qiáng)降雨和重力作用下會形成“滑坡/崩塌—碎屑流—泥石流”等鏈?zhǔn)酱紊鸀?zāi)害,因此,在震后及時準(zhǔn)確地基于GIS平臺開展滑坡易發(fā)性預(yù)測建模,有助于滑坡災(zāi)害防治及震后災(zāi)害演化機(jī)制研究工作的開展[8]。吉日伍呷等[9]以魯?shù)榈卣馂槔?選取邏輯回歸、K近鄰、樸素貝葉斯和隨機(jī)森林算法分別構(gòu)建了4種滑坡易發(fā)性評價模型,并對模型的預(yù)測精度進(jìn)行了對比。針對九寨溝Ms7.0地震同震滑坡災(zāi)害,張玘愷等[10]比較了信息量、確定性系數(shù)、邏輯回歸模型及其耦合模型的快速評估方法對于九寨溝縣范圍內(nèi)滑坡災(zāi)害開展易發(fā)性評價的適用性;羅路廣等[11-14]基于層次分析法、確定性系數(shù)、邏輯回歸模型、信息量模型、廣義加權(quán)統(tǒng)計模型等不同評價方法和評價單元開展了評價模型的比較與耦合、歷史地震情景滑坡發(fā)生概率反演、因子組合選取及風(fēng)險評價等研究。上述研究結(jié)果均表明,數(shù)理統(tǒng)計模型能較好地評價九寨溝地區(qū)滑坡災(zāi)害易發(fā)性,然而機(jī)器學(xué)習(xí)模型是否具有適用性仍有待驗證。針對此問題,本文選取4種代表性機(jī)器學(xué)習(xí)模型,選取地震因子、地形地貌因子、地質(zhì)因子、水文因子以及人類工程活動等五個方面共14個滑坡初始環(huán)境因子,在因子共線性診斷的基礎(chǔ)上,采用10折交叉驗證數(shù)據(jù)取樣方法、受試者特征ROC曲線和易發(fā)性指數(shù)分布特征(均值和標(biāo)準(zhǔn)差)等來探討基于機(jī)器學(xué)習(xí)的滑坡易發(fā)性建模及其不確定性。以期對九寨溝地區(qū)震后短期或長期防災(zāi)減災(zāi)及推廣機(jī)器學(xué)習(xí)模型在其他地區(qū)的易發(fā)性預(yù)測建模具有一定的參考價值。

1 機(jī)器學(xué)習(xí)模型

本文利用邏輯回歸(LR)、支持向量機(jī)(SVM)、隨機(jī)森林(RF)、人工神經(jīng)網(wǎng)絡(luò)(ANN)4種機(jī)器學(xué)習(xí)模型預(yù)測輸出震后滑坡易發(fā)性指數(shù)p(范圍介于0和1之間),并分別以p≥0.5或p<0.5對滑坡發(fā)生和滑坡不發(fā)生進(jìn)行劃分?；赗統(tǒng)計軟件強(qiáng)大的函數(shù)功能和編程算法構(gòu)建機(jī)器學(xué)習(xí)分類模型,并通過10折交叉驗證(10-fold cross validation)取樣方法對數(shù)據(jù)樣本進(jìn)行模型訓(xùn)練和測試,其原理詳見文獻(xiàn)[13]。

1.1 邏輯回歸(LR)

邏輯回歸(logistic regression,LR)利用因變量與多個自變量之間的線性回歸關(guān)系,分析預(yù)測滑坡發(fā)生概率,是應(yīng)用較為廣泛的滑坡區(qū)域評價模型之一[2,5,11]。該模型描述二元因變量(滑坡是否發(fā)生)和自變量(x1,x2,…,xk)之間的關(guān)系,自變量可以是連續(xù)的也可以是離散的,不需要滿足正態(tài)分布,表達(dá)式為:

ln[p/(1-p)]=β0+β1x1+β2x2…+βkxk

(1)

(2)

式中:p為滑坡發(fā)生概率;β0為截距常數(shù)項;β1,β2,…,βk為各滑坡因子所對應(yīng)的回歸系數(shù);x1,x2,…,xk為各滑坡因子指標(biāo)值。

1.2 支持向量機(jī)(SVM)

支持向量機(jī) (support vector machine,SVM)通過在高維空間中找到分類間隔最大化的分離超平面對數(shù)據(jù)進(jìn)行分類[5,15]。假設(shè)輸入數(shù)據(jù)為xi(i=1, 2, 3, ……,n),對應(yīng)二分類問題的輸出為y(y=±1)。

(3)

yi=(ω·xi)+b≥1

(4)

(ω·xi)+b≥1+ξi

(5)

式中: ‖ω‖2為超平面法向量的范數(shù);b為常數(shù);L為拉格朗日定義的損失函數(shù);λi為拉格朗日乘子;ξi為松弛因子。下一步可將式(3)表示為式(6),其中υ(0,1)為誤分類問題。另外,本文選擇徑向基函數(shù)作為SVM的核函數(shù):

(6)

1.3 隨機(jī)森林(RF)

隨機(jī)森林(random forest,RF)模型是由多個決策樹集合而成的集成分類模型,利用隨機(jī)方法構(gòu)建決策樹,其原理是首先利用bootstrap抽樣從原始訓(xùn)練集中有放回地抽取K個樣本,且各樣本的特征數(shù)都與原始訓(xùn)練集相同;再分別對K個樣本建立決策樹模型,得到K種分類結(jié)果;各樣本種隨機(jī)選取n(n≤m)個特征作為分裂特征集,從中選擇最優(yōu)特征對節(jié)點進(jìn)行生長,當(dāng)n≤m時,每一棵決策樹之間又存在差異性。節(jié)點分裂以基尼系數(shù)GI作為雜質(zhì)函數(shù),如式(7)所示:

(7)

式中:c為分類類別個數(shù);t為決策樹的節(jié)點;p為c的相對頻率。最終,形成隨機(jī)森林,根據(jù)K種分類結(jié)果進(jìn)行投票表決以決定其最終分類。由于該模型每個樹的訓(xùn)練樣本及節(jié)點分裂屬性均為隨機(jī)選取,使得在一定程度上避免了過擬合問題。因此,利用RF方法進(jìn)行滑坡易發(fā)性評價逐漸受到重視[2]。

1.4 人工神經(jīng)網(wǎng)絡(luò)(ANN)

人工神經(jīng)網(wǎng)絡(luò)(artificial neural network,ANN)以數(shù)學(xué)模型模擬神經(jīng)元活動,是基于模仿大腦神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)和功能而建立的一種信息處理系統(tǒng)[16]。ANN模型不需要知道輸入、輸出之間的確切關(guān)系,不需大量參數(shù),只需要知道引起輸出變化的非恒定因素,即非常量性參數(shù)。神經(jīng)網(wǎng)絡(luò)技術(shù)在處理模糊數(shù)據(jù)、隨機(jī)性數(shù)據(jù)、非線性數(shù)據(jù)方面具有明顯優(yōu)勢,對規(guī)模大、結(jié)構(gòu)復(fù)雜、信息不明確的系統(tǒng)尤為適用。

2 研究區(qū)概況與滑坡環(huán)境因子

2.1 研究區(qū)概況

本文以2017年8月8日21時19分46秒四川省北部阿壩州九寨溝縣漳扎鎮(zhèn)發(fā)生的Ms7.0地震作為研究案例,九寨溝是長江水系嘉陵江西源白龍江右支流白水江上游的一條支溝,如圖1(a)所示,因溝內(nèi)有9個藏族村寨而得名。它地處青藏高原與四川盆地過渡的深切割高山峽谷地帶,位于青海-西藏和揚(yáng)子板塊的邊緣,平均海拔3000 m以上,地勢南高北低,地質(zhì)背景復(fù)雜,褶皺斷裂發(fā)育,新構(gòu)造運(yùn)動強(qiáng)烈,地殼抬升幅度大,形成了形態(tài)各異的冰斗、U形谷、懸谷、喀斯特等多種地貌和以高山湖泊群、鈣華瀑群及鈣華灘流為特色的自然景觀。九寨溝以自然、文化獨(dú)特的科研價值、保護(hù)價值及旅游觀賞價值等突出普遍價值,于1992年作為自然遺產(chǎn)被聯(lián)合國教科文組織世界遺產(chǎn)委員會列入《世界遺產(chǎn)名錄》。研究區(qū)總面積655.49 km2,受九寨溝“8·8”地震影響,區(qū)內(nèi)斜坡巖體松動、破碎,在降雨作用下極易失穩(wěn)產(chǎn)生滑坡災(zāi)害。

圖1 研究區(qū)概況

2.2 數(shù)據(jù)來源與處理

滑坡易發(fā)性預(yù)測研究的基礎(chǔ)工作包括建立滑坡數(shù)據(jù)庫、劃分評價單元以及提取滑坡環(huán)境因子等。本文的主要數(shù)據(jù)來源包括:①空間分辨率分別為0.2 m的Gaofen-1(日期:2017年8月16日)和3.0 m的Landsat-8(日期:2020年9月20日)衛(wèi)星遙感影像,經(jīng)過影像裁剪、幾何校正、正射矯正、色彩和飽和度調(diào)整等,在ArcGIS軟件中對比分析進(jìn)行目視解譯,建立了研究區(qū)九寨溝“8·8”地震震后滑坡數(shù)據(jù)庫,圈定364處滑坡,如圖1(a)所示,主要分布于研究區(qū)西北部震中附近及沿溝谷發(fā)育,滑坡總平面面積為2.11×106m2,最大、最小面積分別為5.39×104m2和3.89×102m2,規(guī)模以中、小型為主,經(jīng)野外調(diào)查驗證發(fā)現(xiàn)該滑坡編錄數(shù)據(jù)具有較高的質(zhì)量和完整性;②從美國地質(zhì)調(diào)查局網(wǎng)站下載的分辨率為12.5 m的數(shù)字高程數(shù)據(jù)(digital elevation model,DEM)。目前常用的滑坡區(qū)域評價單元主要有柵格單元和斜坡單元[17],斜坡單元以山脊和山谷線為水文分界線,是真實反映滑坡發(fā)生的基本單元。根據(jù)DEM水文分析模型和斜坡坡向?qū)⒀芯繀^(qū)劃分為1234個斜坡單元,如圖1(b)所示,其中有163個斜坡單元內(nèi)有滑坡發(fā)生,剩余斜坡單元則沒有滑坡發(fā)生;③滑坡環(huán)境因子和地形地貌數(shù)據(jù)在ArcGIS軟件中利用表面分析工具從DEM中提取;地質(zhì)因子通過矢量化1∶50000地質(zhì)圖得到地層巖性、巖層產(chǎn)狀及斷層信息;降雨量數(shù)據(jù)來自于研究區(qū)氣象水文觀測站多年平均數(shù)據(jù),采用克里金插值法獲得;公路、水系均通過衛(wèi)星遙感影像矢量化獲取。

2.3 滑坡環(huán)境因子

基于數(shù)據(jù)的可獲取性和前人研究結(jié)果[10,13],本文從地震因子、地形地貌因子、地質(zhì)因子、水文因子以及人類工程活動等五個方面共選取了14個滑坡初始環(huán)境因子,其空間分布特征如圖2所示。地震因子包括地震動峰值加速度PGA(圖2(a))、發(fā)震斷層距(圖2(b))、震中距(圖2(c))等,能夠反映九寨溝地震動參數(shù)對斜坡穩(wěn)定性產(chǎn)生的影響,其中以虎牙斷裂的北西延伸段為發(fā)震斷層[18];地形地貌因子主要包括海拔高程(圖2(d))、斜坡坡度(圖2(e))、坡向(圖2(f))、剖面曲率(圖2(g))以及平面曲率(圖2(h));地質(zhì)因子有斷層距(圖2(i))、坡體結(jié)構(gòu)(圖2(m))、地層巖性(圖2(n))等,其中坡體結(jié)構(gòu)可劃分為橫向坡、逆向坡、逆斜向坡、順向坡和順斜向坡等5個類型;水文因子有距水系距離(圖2(j))和降雨量(圖2(l));道路反映了人類工程活動活躍范圍,以距公路距離(圖2(k))來表示。相關(guān)因子距離分析均利用了ArcGIS的歐式距離分析工具。所有因子圖層分辨率均被重采樣為12.5 m×12.5 m,各指標(biāo)值通過計算斜坡單元內(nèi)柵格值得平均值獲取。

圖2 研究區(qū)滑坡環(huán)境因子

3 滑坡易發(fā)性評價結(jié)果

3.1 因子相關(guān)性分析結(jié)果

滑坡易發(fā)性預(yù)測中為防止因子間存在共線性,可能增加模型的運(yùn)行時間并降低模型精度,因此,在訓(xùn)練機(jī)器學(xué)習(xí)模型之前,需要對各因子進(jìn)行相關(guān)性分析以驗證其獨(dú)立性[19]。圖3為基于R統(tǒng)計分析軟件得到的滑坡易發(fā)性評價因子皮爾遜相關(guān)系數(shù)矩陣,相關(guān)系數(shù)r取值介于-1和1之間,其中r>0為正相關(guān),表明一個因子隨另一個因子的增加而增加;r<0為負(fù)相關(guān),情況則剛好相反。從圖3可知,各因子間基本不存在共線性或相關(guān)性較弱可以忽略不計,然而地震動峰值加速度與震中距、降雨量均呈現(xiàn)較高的負(fù)相關(guān)性,相關(guān)性系數(shù)r分別為-0.96和-0.75,震中距與降雨量則呈現(xiàn)一定的正相關(guān)性(r=0.73)。綜合考慮后,剔除震中距和距發(fā)震斷層距離,僅保留地震動峰值加速度這一地震誘發(fā)因素。最終采用地震動峰值加速度、高程、坡度、坡向、剖面曲率、平面曲率、坡體結(jié)構(gòu)、地層巖性、斷層距、距公路距離、距水系距離和降雨量進(jìn)行后續(xù)滑坡易發(fā)性建模。

圖3 滑坡易發(fā)性評價因子皮爾遜相關(guān)系數(shù)矩陣

3.2 震后滑坡易發(fā)性制圖

圖4(a)～(d)分別為基于LR、SVM、RF和ANN模型得到的研究區(qū)震后滑坡易發(fā)性指數(shù)空間分布結(jié)果,其中LR模型反演得到的易發(fā)性指數(shù)為0～0.966,SVM模型的易發(fā)性指數(shù)介于0和0.671之間,RF模型的易發(fā)性指數(shù)為0～0.952,ANN模型得到的易發(fā)性指數(shù)分布于0～0.766。從圖4可知,研究區(qū)大部分地區(qū)易發(fā)性指數(shù)均較低,易發(fā)性指數(shù)較高區(qū)靠近震中附近和沿溝谷發(fā)育分布,與震后滑坡編錄結(jié)果一致,而年降雨量對滑坡災(zāi)害的影響較小,表明地震波對斜坡結(jié)構(gòu)產(chǎn)生震損是滑坡易發(fā)的主要因素。

圖4 4種評價模型易發(fā)性指數(shù)分布結(jié)果

3.3 易發(fā)性指數(shù)分布規(guī)律

圖4(e)為基于LR、SVM、RF和ANN模型得到的研究區(qū)震后滑坡易發(fā)性指數(shù)分布規(guī)律,為分析模型的不確定性,本文采用均值和標(biāo)準(zhǔn)差分別反映易發(fā)性指數(shù)分布的平均水平和離散程度。LR模型和RF模型的統(tǒng)計結(jié)果較為一致,在易發(fā)性指數(shù)值較低的區(qū)間分布集中,隨著值增大而逐漸減少,SVM模型的易發(fā)性指數(shù)分布呈現(xiàn)先增加后減小的特征,而ANN模型則表現(xiàn)為“跳躍式”分布。各模型易發(fā)性指數(shù)p均值相差不大,其值從小至大依次為:SVM(0.093)

3.4 模型精度比較

本文基于ROC曲線(receiver operating characteristic curve, ROC)和曲線下面積(area under curve,AUC)來進(jìn)行機(jī)器學(xué)習(xí)模型性能評價。ROC曲線的橫、縱坐標(biāo)分別為1-特異性和敏感度,當(dāng)曲線越靠近左上角,AUC值越接近1,表明該模型精度越高。根據(jù)已有的評判標(biāo)準(zhǔn)[20],當(dāng)0.50.9則表明模型性能出色。

4種機(jī)器學(xué)習(xí)模型在研究區(qū)滑坡易發(fā)性評價的ROC曲線如圖5所示。從圖5可知, LR、SVM、RF、ANN模型均具有較好或出色的性能,準(zhǔn)確率AUC值均超過0.88,除RF模型外,其余3個模型的AUC值均超過0.9。準(zhǔn)確率AUC值從大至小依次為:AUCANN>AUCSVM>AUCLR>AUCRF?；?0折交叉驗證取樣方法對模型進(jìn)行10次測試,獲得10次預(yù)測精度CVs,如表1所示。圖6為各模型預(yù)測精度ROC曲線和CVs分布箱型圖,各模型的預(yù)測精度均值CVmean從大至小依次為:SVM(0.925)>ANN(0.920)>LR(0.894)>RF(0.877),可以看出,SVM模型的預(yù)測性能最好,RF模型的預(yù)測性能最差。根據(jù)已有評價標(biāo)準(zhǔn),LR和RF模型具有較好的預(yù)測精度,而SVM和ANN模型的預(yù)測性能被判定為出色。此外,CVs標(biāo)準(zhǔn)差CVstd從小至大依次為:SVM(0.023)

表1 基于10折交叉驗證的各模型預(yù)測精度

圖5 各評價模型ROC曲線與準(zhǔn)確率AUC值

圖6 基于10折交叉驗證的各評價模型ROC曲線與預(yù)測精度CV

4 結(jié)論

1)通過皮爾遜相關(guān)系數(shù)分析各滑坡易發(fā)性評價因子間的相關(guān)性,剔除震中距和距發(fā)震斷層距離,僅保留地震動峰值加速度這一地震誘發(fā)因素,最終采用12個滑坡影響因子進(jìn)行滑坡易發(fā)性建模。

2)利用4種機(jī)器學(xué)習(xí)模型和斜坡單元開展九寨溝地區(qū)震后滑坡,結(jié)果表明,易發(fā)性指數(shù)較高區(qū)靠近震中附近和沿溝谷發(fā)育分布,與震后滑坡編錄結(jié)果一致,表明地震波對斜坡結(jié)構(gòu)產(chǎn)生震損是滑坡易發(fā)的主要因素。

3)基于10折交叉驗證取樣方法的ROC曲線和易發(fā)性指數(shù)分布特征結(jié)果表明,SVM模型是在九寨溝地區(qū)震后滑坡預(yù)測中性能表現(xiàn)最好的模型,ANN模型次之,LR模型緊隨其后,RF模型效果最差。該結(jié)果可為震后滑坡減災(zāi)工作和推廣機(jī)器學(xué)習(xí)模型在其他地區(qū)的易發(fā)性預(yù)測建模提供理論依據(jù)與參考。