張向營，張春山，孟華君，王雪冰，趙偉康,4，鄭滿城,4(1.中國地質(zhì)科學(xué)院地質(zhì)力學(xué)研究所， 北京 100081；2.自然資源部新構(gòu)造運(yùn)動與地質(zhì)災(zāi)害重點(diǎn)實(shí)驗(yàn)室，北京 100081；3.北京城建勘測設(shè)計(jì)院有限責(zé)任公司，北京 100101；4.中國地質(zhì)大學(xué)(北京)地球科學(xué)與資源學(xué)院，北京 100083 )

滑坡是我國山區(qū)環(huán)境中破壞最嚴(yán)重的地質(zhì)災(zāi)害之一[1]，開展滑坡危險(xiǎn)性區(qū)劃工作對實(shí)現(xiàn)防災(zāi)減災(zāi)、災(zāi)害管理等具有重要的現(xiàn)實(shí)意義。

自20世紀(jì)60年代GIS技術(shù)出現(xiàn)后，國內(nèi)外學(xué)者對滑坡危險(xiǎn)性區(qū)劃方法進(jìn)行了積極的探索，至今已發(fā)展了近幾十種評價(jià)模型[2]，主要包括信息量模型[3]、邏輯回歸模型[4]、證據(jù)權(quán)重模型[5]、貢獻(xiàn)率模型[6]、聚類分析[7]、神經(jīng)網(wǎng)絡(luò)模型[8]等客觀評價(jià)模型及層次分析法等主觀評價(jià)模型[9～10]。這些成果豐富了滑坡危險(xiǎn)性的分析與評價(jià)，但也都存在不同程度的缺陷，如主觀評價(jià)模型存在人為主觀性、模型普適性等問題；客觀評價(jià)模型則存在樣本數(shù)據(jù)點(diǎn)要求較高、指標(biāo)選取缺乏科學(xué)依據(jù)等問題。隨機(jī)森林算法[11](Random Forest，RF)具有很高的學(xué)習(xí)能力和預(yù)測精度，已廣泛應(yīng)用于醫(yī)學(xué)、自然生態(tài)學(xué)、市場經(jīng)濟(jì)學(xué)等統(tǒng)計(jì)領(lǐng)域[12]，在滑坡危險(xiǎn)性區(qū)劃領(lǐng)域中也有應(yīng)用[13～14]，但使用者均局限于模型本身，并未克服客觀賦權(quán)法存在的共性問題。

為探索更為實(shí)用的組合評價(jià)形式，選取應(yīng)用較為成熟的層次分析法(Analytic Hierarchy Process，AHP)和隨機(jī)森林算法分別得到致災(zāi)因子的主、客觀權(quán)重，然后基于距離函數(shù)計(jì)算得到組合權(quán)重，構(gòu)建組合評價(jià)模型，并以青海省貴德縣北部地區(qū)為研究區(qū)，對3種模型區(qū)劃結(jié)果進(jìn)行對比分析，探討組合模型的可靠度。

1 研究方法

1.1 隨機(jī)森林算法

隨機(jī)森林算法集成了Bagging與Random subspace兩種主流學(xué)習(xí)方法，與其他統(tǒng)計(jì)學(xué)學(xué)習(xí)方法相比，兩個(gè)隨機(jī)性的引入，使得隨機(jī)森林算法不容易陷入過擬合，在運(yùn)算量沒有顯著提高的前提下提高了預(yù)測精度，而且在異常值和噪聲方面具有較高的容忍度，結(jié)果對缺失數(shù)據(jù)和非平衡的數(shù)據(jù)比較穩(wěn)健[11]。其主要步驟見圖1。

圖1 隨機(jī)森林方法步驟[19]Fig.1 Schematic illustration of the Random Forest method

(1)因子選擇及多元共線性診斷

在對復(fù)雜決策系統(tǒng)充分了解的基礎(chǔ)上，劃分出層次結(jié)構(gòu)及其各因子的從屬關(guān)系，在滑坡危險(xiǎn)性區(qū)劃中，篩選出n個(gè)滑坡危險(xiǎn)性評價(jià)指標(biāo)(x1,x2,…,xn)。為避免多個(gè)變量高相關(guān)，造成模型方程的不穩(wěn)定，構(gòu)建隨機(jī)森林模型之前還應(yīng)對各因子進(jìn)行多元共線性診斷。常用的多元共線性考核指標(biāo)包括因子相關(guān)矩陣的條件數(shù)κ和方差膨脹因子(Variance Inflation Factor)兩類[15]。同時(shí)，為了減少單一取樣方式對模型結(jié)果產(chǎn)生的影響，應(yīng)將總體樣本數(shù)據(jù)(由正樣本和負(fù)樣本等量組成)按一定比例(α︰1-α)隨機(jī)劃分為訓(xùn)練樣本集和檢驗(yàn)樣本集兩類。

(2)模型建立及檢驗(yàn)

隨機(jī)森林模型的創(chuàng)建及結(jié)果驗(yàn)證部分是在R軟件環(huán)境下進(jìn)行的，它可以提供一些集成的統(tǒng)計(jì)工具，如本次建模需要的 “RandomForest”、“caret” 工具包。對于建模效果的評價(jià)，模型運(yùn)算的不同階段中所用的檢驗(yàn)指標(biāo)也不盡相同，通常包括混淆矩陣(Confusion matrix)驗(yàn)算、接受靈敏度曲線(Receiver-Operating Characteristic,ROC)、曲線下面積(Area Under Curve)、泛化能力和預(yù)報(bào)效率曲線[16]等指標(biāo)。

(3)結(jié)果分析及客觀權(quán)重確定

各建模因子對于滑坡危險(xiǎn)性的貢獻(xiàn)和影響規(guī)律各不相同，對各變量貢獻(xiàn)及影響規(guī)律進(jìn)行量化分析，可以為后期滑坡災(zāi)害管理、防災(zāi)減災(zāi)工作提供重要指導(dǎo)。隨機(jī)森林方法的一個(gè)重要特性就是能夠計(jì)算每個(gè)變量的重要性值，RF提供兩種基本的變量重要性值：平均Gini指數(shù)降低度(Mean Decrease Gini)和平均準(zhǔn)確率降低度(Mean Decrease Accuracy)[11]。其中，平均準(zhǔn)確率降低度為袋外數(shù)據(jù)自變量值發(fā)生輕微擾動后的分類正確率與擾動前分類正確率的平均減少量；平均Gini指數(shù)降低度是用來表示所有樹的變量分割節(jié)點(diǎn)平均減小的不純度。本次評價(jià)把平均Gini指數(shù)降低度作為相應(yīng)因子的客觀權(quán)重值ωi(i=1,2, …,n)。

1.2 層次分析法

層次分析法(AHP)是一種定性與定量分析方法相結(jié)合的綜合性評價(jià)方法，目前在地質(zhì)災(zāi)害危險(xiǎn)性評價(jià)領(lǐng)域中是應(yīng)用最為成熟的模型之一[17]，其主要步驟如下：

(1)構(gòu)建判斷矩陣U

在RF模型篩選的滑坡危險(xiǎn)性評價(jià)指標(biāo)基礎(chǔ)上，為了量化各評價(jià)因子的主觀賦值權(quán)重，評價(jià)因子之間需形成一個(gè)兩兩對比的矩陣，利用1～9標(biāo)度表示評價(jià)因子之間的影響強(qiáng)弱(表1)，通過專家賦值法確定各判別因子之間的相對重要性并賦與相應(yīng)的分值，形成行、列相等的對稱矩陣U。

表1 1～9標(biāo)度的意義表Table.1 Meaning of scales 1 to 9 in the judgment matrix

(2)確定主觀權(quán)重

1.3 組合賦權(quán)

層次分析法計(jì)算出的權(quán)重作為主觀權(quán)重，隨機(jī)森林法計(jì)算出的權(quán)重作為客觀權(quán)重，利用特定的數(shù)學(xué)理論公式，使主、客觀權(quán)重有機(jī)地結(jié)合起來稱之為組合賦權(quán)。目前，關(guān)于組合賦權(quán)的方法主要包括基于乘法歸一化的組合賦權(quán)法[18]、線性加權(quán)組合法[19]、距離函數(shù)法[20]等。線性加權(quán)法主要是根據(jù)主觀實(shí)際經(jīng)驗(yàn)判斷分配系數(shù)，人為干擾相對強(qiáng)烈；乘法歸一化則是對各權(quán)重與其對應(yīng)的指標(biāo)值之間的差異程度進(jìn)行優(yōu)化，但當(dāng)數(shù)據(jù)特征差異性較大時(shí)，可靠度會顯著下降；而距離函數(shù)法引入了距離函數(shù)的概念，將主、客觀權(quán)重之前的差異程度與其對應(yīng)的分配系數(shù)間的差異程度一致化，有機(jī)地結(jié)合起來，既兼并了評判者對實(shí)際情況的主觀經(jīng)驗(yàn)，又使組合權(quán)重具有統(tǒng)計(jì)學(xué)意義。因此，筆者采用了距離函數(shù)法進(jìn)行組合賦權(quán)，具體步驟如下：

設(shè)由層次分析法確定的主觀權(quán)重為ai，由隨機(jī)森林法確定的客觀權(quán)重為ωi，兩者的距離函數(shù)為d(ai，ωi)[20]，其表達(dá)式為：

(1)

設(shè)組合權(quán)重為φi，組合權(quán)重值為二者的線性加權(quán)，表達(dá)式為：

φi=μai+νωi

(2)

式中：μ、v——主、客觀權(quán)重的分配系數(shù)。

為了使不同權(quán)重之間的差異程度和分配系數(shù)間的差異程度一致，使分配系數(shù)與式(1)中的距離函數(shù)取等式，其表達(dá)式為：

d(ai,ωi)2=(μ-ν)2

(3)

μ+ν=1

(4)

聯(lián)立式(3)、(4)，即可得組合賦權(quán)的分配系數(shù)μ、v。

(5)

ν=1-μ

(6)

再將μ、v帶入式(2)，即可得到綜合權(quán)重值φi。

1.4 綜合評價(jià)模型

本文滑坡危險(xiǎn)性評價(jià)方法主要采用綜合指數(shù)評價(jià)法，由綜合評價(jià)指數(shù)計(jì)算獲取，其計(jì)算公式為：

(7)

式中：W——滑坡災(zāi)害危險(xiǎn)性綜合指數(shù)；

ai——評價(jià)因子；

φi——綜合權(quán)重。

根據(jù)上述滑坡災(zāi)害相關(guān)性影響因子統(tǒng)計(jì)分析結(jié)果，分別使用主觀權(quán)重、客觀權(quán)重和綜合權(quán)重，在ArcGIS軟件中把各因子?xùn)鸥駡D層進(jìn)行空間疊加運(yùn)算，分別得到了研究區(qū)滑坡危險(xiǎn)性分區(qū)圖，在此基礎(chǔ)上利用自然斷點(diǎn)法將危險(xiǎn)性區(qū)劃圖各自分為四類：低危險(xiǎn)性、中危險(xiǎn)性、高危險(xiǎn)性、極高危險(xiǎn)性。

2 實(shí)例分析

2.1 研究區(qū)概況

貴德縣北部山區(qū)地處青海省海南藏族自治州東部，處于黃河上游龍羊峽與李家峽之間，黃河由西向東橫貫其中，區(qū)內(nèi)溝壑縱橫，中部呈多級河流階地和盆地丘陵地貌，南北兩側(cè)為高山峽谷地貌。氣候?qū)俅箨懶愿咴瓪夂?，特點(diǎn)是“冬寒夏短日照長，風(fēng)多雨少、太陽輻射強(qiáng)”，年降水量在251～559 mm之間。大地構(gòu)造位置處于南祁連塊體和西秦嶺塊體的交接部位，受構(gòu)造運(yùn)動影響，區(qū)內(nèi)主要為NWW向走滑斷裂和逆沖斷裂。巖性包括松散第四系沉積物，新近系鈣質(zhì)泥巖、粉砂巖及砂巖及零星出露的早古生界和中古生界侵入巖，其中，新近系地層為區(qū)內(nèi)主要的易滑地層。

2.2 數(shù)據(jù)來源及因子選擇

選擇地形地貌、年降雨量、地質(zhì)條件、水文、人類工程活動、植被等6大類10個(gè)指標(biāo)作為研究區(qū)滑坡危險(xiǎn)性影響因子(表2)。其中，植被歸一化指數(shù)(NDVI)為利用Landsat8遙感影像在ENVI軟件中提取生成，其余因子都在ArcGIS軟件中處理生成。表2中專家認(rèn)可度這一指標(biāo)主要是在層次分析法中構(gòu)造判斷矩陣時(shí)使用。

表2 滑坡因子說明表Table 2 Explaination of influencing factors

2.3 客觀權(quán)重

(1)多元線性診斷

本文建模所需的樣本集包括正樣本和負(fù)樣本兩部分，正樣本集為野外實(shí)地調(diào)查的158個(gè)滑坡點(diǎn)，負(fù)樣本集是為維持?jǐn)?shù)據(jù)樣本的平衡性，利用ArcGIS隨機(jī)取樣工具，在已知滑坡點(diǎn)100 m外隨機(jī)選取的158個(gè)非滑坡點(diǎn)組成的樣本集。并把樣本總集隨機(jī)劃分為70%的樣本數(shù)據(jù)作為訓(xùn)練樣本集和30%的樣本數(shù)據(jù)作為檢驗(yàn)樣本集。

在獲取總樣本后，基于R軟件平臺進(jìn)行了因子多元線性診斷，得到各因子特征之間的相關(guān)性系數(shù)。發(fā)現(xiàn)坡度與高差之間的特征相關(guān)性系數(shù)高達(dá)0.68，呈較高的相關(guān)性，而其他因子相關(guān)性系數(shù)均小于0.25，因此，為維持模型的結(jié)果客觀性，本次危險(xiǎn)性區(qū)劃剔除了高差這個(gè)因子。

圖2 研究區(qū)滑坡位置點(diǎn)Fig.2 Location of landslide in the study area

(2) 模型建立及評價(jià)

參數(shù)mtry是隨機(jī)森林建模中構(gòu)建決策樹分支時(shí)隨機(jī)抽樣的變量個(gè)數(shù)[21]，通過遍歷設(shè)定mtry為1至9，并進(jìn)行9 次建模，打印出每次建模的錯(cuò)誤率，選擇錯(cuò)誤率最低的mtry取值為3。通過模型錯(cuò)誤率與決策樹數(shù)量的關(guān)系可視化，本次模型選擇以mtry=3，ntree=150進(jìn)行隨機(jī)森林建模。建立模型后，選擇Kappa系數(shù)[13]對模型結(jié)果進(jìn)行了評價(jià)檢驗(yàn)。本次模型的混淆矩陣見表3。

表3 混淆矩陣Table 3 Confusion matrix

從表2可以看出模型對于測試樣本集的預(yù)測情況，實(shí)際未滑動而預(yù)測滑動的樣本個(gè)數(shù)為8，實(shí)際發(fā)生滑動而預(yù)測未滑動的樣本個(gè)數(shù)為6。由Kappa系數(shù)的計(jì)算公式[13]可得，模型的Kappa系數(shù)為0.82，總體來說，模型預(yù)測精度較高。

隨機(jī)森林模型有兩種衡量滑坡因子相對重要性的方式：①將某滑坡因子的取值隨機(jī)打亂。分析其打亂前后隨機(jī)森林的準(zhǔn)確性和錯(cuò)誤率，模型錯(cuò)誤率提升越多的因子越重要，即平均準(zhǔn)確率降低度(Mean Decrease Accuracy)；②計(jì)算某因子對森林中所有決策樹節(jié)點(diǎn)不純度的影響，即平均Gini指數(shù)降低度(Mean Decrease Gini)，該值越大表示該變量的重要性越大。通過R語言的importance命令，得到了9個(gè)滑坡因子的重要性評價(jià)，結(jié)果如圖3所示。

圖3 因子重要性分析結(jié)果Fig.3 Factor importance evaluation

從兩種因子重要性分析方法得出的結(jié)論來看，人類工程活動、河流侵蝕、坡度、降雨量都是對滑坡危險(xiǎn)性貢獻(xiàn)最多的四種因子；排名靠后的幾個(gè)因子在兩種方法中順序雖有所變化，但考慮這幾種因子比例相差不大，加上其相對重要性較低，因此，認(rèn)為這種排序整體上是合理的。此次評價(jià)將平均Gini指數(shù)降低度作為客觀權(quán)重ωi(表4)。

2.4 主觀權(quán)重

(1)構(gòu)建判斷矩陣U

在RF模型篩選的滑坡危險(xiǎn)性評價(jià)指標(biāo)基礎(chǔ)上，為了量化各評價(jià)因子的主觀賦值權(quán)重，評價(jià)因子之間需形成一個(gè)兩兩對比的矩陣，通過專家賦值法，利用1～9標(biāo)度表示評價(jià)因子之間的影響強(qiáng)弱(表3)，確定判斷矩陣U。

(2)AHP法確定主觀權(quán)重

ai即為各因子對應(yīng)的主觀權(quán)重值，然后計(jì)算出最大特征值：

最后利用CR=CI/RI=(λmax-n)/RI(n-1) 進(jìn)行一致性檢驗(yàn)，由文[17]可知n=9時(shí)，RI=1.46。于是可知CR0.0392<0.1，表明主觀權(quán)重計(jì)算合理。

2.5 組合權(quán)重

將上述所求的主觀權(quán)重ai和客觀權(quán)重ωi分別代入距離函數(shù)法計(jì)算公式(5)和(6)，可得：

ν=1-μ=0.42

將μ、v帶入式(2)，即可得到各因子組合權(quán)重值φi(表4)。

表4 權(quán)重賦值表Table 4 Weight assignment table

2.6 危險(xiǎn)性評價(jià)

將主觀權(quán)重、客觀權(quán)重及組合權(quán)重分別代入綜合評價(jià)模型(式7)，得到三種研究區(qū)滑坡危險(xiǎn)性分區(qū)圖(圖4)。

圖4 滑坡危險(xiǎn)性區(qū)劃圖((a)為層次分析法；(b)為隨機(jī)森林法；(c)為組合賦權(quán)法)Fig.4 Landslide susceptibility map(A：AHP；B：RF；C:RF-AHP)

3 危險(xiǎn)性結(jié)果與討論

3.1 危險(xiǎn)性結(jié)果

得到全區(qū)滑坡危險(xiǎn)性區(qū)劃圖后，統(tǒng)計(jì)得到了三種模型的危險(xiǎn)性等級與實(shí)際災(zāi)害分布的對比關(guān)系(表5)。

3.2 討論

對比分析三種方法得出的危險(xiǎn)性評價(jià)結(jié)果可知：

AHP模型中的極高和中等危險(xiǎn)性區(qū)域與已有滑坡災(zāi)害點(diǎn)的分布吻合性最好，但其危險(xiǎn)性和極高危險(xiǎn)性區(qū)域面積也是最高的，達(dá)到了1 365.6 km2。

RF模型結(jié)果中的低危險(xiǎn)區(qū)、中危險(xiǎn)區(qū)及高危險(xiǎn)區(qū)面積均超過26%，且各區(qū)域內(nèi)包括的災(zāi)害點(diǎn)數(shù)量沒有較大變化，評價(jià)結(jié)果過于均值化。導(dǎo)致均值化的原因主要是由于RF模型對模型因子的相對不敏感性，模型往往會對不重要的因子高估，對重要的因子低估，導(dǎo)致危險(xiǎn)性分區(qū)相對均值化[22]。

從RF-AHP區(qū)劃結(jié)果來看：在高危險(xiǎn)區(qū)和極高危險(xiǎn)區(qū)面積占比38.38%的情況下，就包括了60.13%的滑坡災(zāi)害，結(jié)果準(zhǔn)確性相比AHP和RF兩種模型有較大提升；而且隨著危險(xiǎn)性等級的逐步提高，災(zāi)害實(shí)際發(fā)生的比率(b/a，亦即等級中的災(zāi)害密度與研究區(qū)總的災(zāi)害密度的比值)隨之增大，說明RF-AHP不但把RF的高學(xué)習(xí)能力、預(yù)測精度和AHP的專家經(jīng)驗(yàn)修正等優(yōu)點(diǎn)進(jìn)行了有效整合，還將兩種模型各自的缺點(diǎn)進(jìn)行優(yōu)化，使得RF-AHP在評價(jià)準(zhǔn)確性及實(shí)際應(yīng)用方面有了較高的提升。

4 結(jié)論

(1)本文基于層次分析模型和隨機(jī)森林模型，通過引入距離函數(shù)，探索性提出了一種新的組合賦權(quán)方法，并對組合模型進(jìn)行了系統(tǒng)性的推導(dǎo)。以貴德縣北部地區(qū)為研究區(qū)，選取了地形地貌、植被覆蓋率、地質(zhì)構(gòu)造、工程巖組、河流侵蝕、降雨量和人類工程活動等9個(gè)評價(jià)指標(biāo)對組合評價(jià)模型進(jìn)行了驗(yàn)證，以Kappa系數(shù)和災(zāi)害實(shí)際發(fā)生比率兩個(gè)指標(biāo)對模型可靠度進(jìn)行了檢驗(yàn)分析。

表5 劃分的危險(xiǎn)性等級與實(shí)際災(zāi)害分布的對比Table 5 Comparison between the evaluation results and the actual disaster distribution

注：(1)a為本類危險(xiǎn)性等級的面積占研究區(qū)總面積的百分比；(2)b為落在該危險(xiǎn)性分區(qū)內(nèi)的災(zāi)害占災(zāi)害點(diǎn)總數(shù)的百分比；(3)c為落在該類危險(xiǎn)性分區(qū)內(nèi)的災(zāi)害數(shù)量。

(2)把評價(jià)結(jié)果同RF和AHP模型結(jié)果進(jìn)行了比較，得到的危險(xiǎn)性區(qū)劃成果在空間分布上與既有歷史滑坡數(shù)據(jù)一致性較好，說明本方法能夠通過主、客觀方法的結(jié)合和組合賦權(quán)規(guī)則的改進(jìn)，得到相對真實(shí)、客觀的滑坡危險(xiǎn)性區(qū)劃結(jié)果，為此類評價(jià)提供了一種新思路。

(3)基于RF-AHP的滑坡災(zāi)害評價(jià)模型同樣存著適用性未知和研究缺陷，由于模型是基于大量的數(shù)據(jù)統(tǒng)計(jì)來進(jìn)行規(guī)律擬合的，在災(zāi)害點(diǎn)樣本量較少、樣本質(zhì)量較低的地區(qū)，該模型是否具有同樣的可行性和精度尚不得而知。