鄭 雪，唐章英，宋 超

1.西南石油大學(xué) 地球科學(xué)與技術(shù)學(xué)院，四川 成都 610500；2.四川大學(xué) 華西公共衛(wèi)生學(xué)院，四川 成都 610041

滑坡是中國山地丘陵地區(qū)發(fā)生頻率最高的自然災(zāi)害之一，常造成多人傷亡和巨大的經(jīng)濟(jì)損失[1]。根據(jù)中國國土資源通報(bào)，2013年中國發(fā)生滑坡災(zāi)害數(shù)量占地質(zhì)災(zāi)害總數(shù)量的63.9%，占比最大，并造成近500人死亡以及大約100億元的直接經(jīng)濟(jì)損失(https：//www.mnr.gov.cn/)；同時(shí)，伴隨著更密集的道路建設(shè)，過度放牧，農(nóng)田耕作灌溉，礦山開采，城市向鄉(xiāng)村擴(kuò)張等人類活動(dòng)對地質(zhì)環(huán)境造成了一定程度的破壞[2]，地質(zhì)環(huán)境變得越來越脆弱，再加上氣候變化改變了全球極端降雨的頻率和強(qiáng)度[3]，使得滑坡、崩塌、泥石流等一系列相關(guān)的地質(zhì)災(zāi)害發(fā)生的概率大大增加?；乱装l(fā)性制圖是防災(zāi)減災(zāi)的一個(gè)重要工具，因?yàn)樗梢灾赋龌碌拇嗳醯貐^(qū)以便于我們有指向性的進(jìn)行災(zāi)害防護(hù)。因此，選擇合適的影響因素和研究方法來保證滑坡災(zāi)害制圖和評估的準(zhǔn)確性是減輕滑坡災(zāi)害損失的關(guān)鍵。

區(qū)域滑坡災(zāi)害易發(fā)性評價(jià)，有些研究也將其稱為滑坡敏感性評價(jià)或是危險(xiǎn)性區(qū)劃，旨在區(qū)域尺度上(通常是一個(gè)較大的區(qū)域)根據(jù)當(dāng)?shù)氐牡乩憝h(huán)境因素確定未來容易發(fā)生滑坡的地區(qū)并進(jìn)行分級，量化不同因素對滑坡發(fā)生的影響，進(jìn)而劃分干預(yù)方式區(qū)別對待，這種方法已是主動(dòng)有效預(yù)警滑坡并減輕災(zāi)害損失的重要手段之一[4]。一般來說，主要有基于機(jī)器學(xué)習(xí)算法的模型和統(tǒng)計(jì)方法的模型?；跈C(jī)器學(xué)習(xí)的模型有支持向量機(jī)[5]、相關(guān)向量機(jī)[6]、人工神經(jīng)網(wǎng)絡(luò)[7]、灰色關(guān)聯(lián)度[8]等。但是這類傳統(tǒng)的機(jī)器學(xué)習(xí)方法往往需要大量的先驗(yàn)知識(shí)和假設(shè)，調(diào)參復(fù)雜，深度不足以完全提取潛在的滑坡特征[9]，不能合理地解釋滑坡空間分布與滑坡風(fēng)險(xiǎn)因子之間的聯(lián)系[10]。統(tǒng)計(jì)方法的模型包括邏輯回歸[11]、貝葉斯空間統(tǒng)計(jì)[12]、頻率比[13]、信息量[14]、確定系數(shù)法[15]等?，F(xiàn)如今，邏輯回歸模型、頻率比模型、信息量模型、確定性系數(shù)模型已經(jīng)有很多的學(xué)者進(jìn)行了大量的研究并應(yīng)用，這些方法比較容易快速實(shí)現(xiàn)，且原理較為簡單。同時(shí)，也有一些文章將這幾種模型進(jìn)行組合并同單一模型進(jìn)行對比，得到了相較于前人好一些的結(jié)果[16]。而貝葉斯空間統(tǒng)計(jì)模型由于考慮到數(shù)據(jù)的空間信息，與當(dāng)前許多研究僅僅只利用數(shù)據(jù)的屬性信息有較大區(qū)別，因?yàn)槟軌蝻@著提升滑坡易發(fā)性預(yù)測的精度而得到推崇。

本文就是基于貝葉斯空間統(tǒng)計(jì)模型的分支——貝葉斯空間邏輯回歸模型，使用當(dāng)前主流的隨機(jī)森林模型篩選風(fēng)險(xiǎn)因子后進(jìn)行建模，對四川雅安蘆山地區(qū)震后滑坡易發(fā)性進(jìn)行分級并與普通未考慮空間結(jié)構(gòu)信息的邏輯回歸建模結(jié)果進(jìn)行比較。

1 研究區(qū)域和數(shù)據(jù)

1.1 研究區(qū)概況

研究區(qū)大部分位于四川省雅安市內(nèi)，跨雅安市市轄區(qū)、寶興縣、蘆山縣、名山縣、滎經(jīng)縣、天全縣，小部分位于成都市的大邑縣和邛崍市，其地理覆蓋范圍接近29°87′-30°62′N 緯度和102°40′-103°30′E 經(jīng)度，覆蓋面積約4 888平方公里(圖1)。該地區(qū)為亞熱帶濕潤季風(fēng)氣候，夏季通常漫長、炎熱且潮濕，冬季涼爽溫和，全年降水量分布大致均勻，植被覆蓋較為茂密。研究區(qū)地形自東南向西北，由平原、低丘陵轉(zhuǎn)向高山，位于青藏高原和四川盆地之間的地形起伏很大的區(qū)域，有巨大的高差和陡峭的山坡。研究區(qū)大部分位于的雅安市是中國城市年降水量最高的城市之一，年均降雨量1 800毫米左右，有雨城之稱(四川省公共氣象服務(wù)網(wǎng)，https：//www.scggqx.com/)。

圖1 研究區(qū)概況Fig.1 General situation

1.2 數(shù)據(jù)來源

研究使用的15 546個(gè)滑坡點(diǎn)清單來自于2013年4月20日中國廬山6.6級地震引發(fā)滑坡數(shù)據(jù)庫。此數(shù)據(jù)庫是基于之前不完整的數(shù)據(jù)庫和補(bǔ)充遙感數(shù)據(jù)，包括地震后的高分辨率Rapid-Eye和ZY-3衛(wèi)星圖像，來建立的一個(gè)更詳細(xì)、更客觀、更完整的蘆山事件的滑坡數(shù)據(jù)庫[17]，并且在選定地區(qū)進(jìn)行了現(xiàn)場調(diào)查來驗(yàn)證所得到的滑坡庫存。

同時(shí)，本研究總共收集處理了16個(gè)影響滑坡發(fā)生的風(fēng)險(xiǎn)因素，分為地質(zhì)地形 (高程、坡向、坡度、到斷層的距離、巖性、地形、土地覆蓋)、生態(tài) (歸一化植被指數(shù)、年平均降水量、地形濕度指數(shù)、到河流的距離)，以及社會(huì)經(jīng)濟(jì) (到道路的距離、到居民點(diǎn)的距離)，地震相關(guān)因素(到震中的距離、地震烈度、峰值地面加速度)。具體來源如表1所示。

表1 風(fēng)險(xiǎn)因子數(shù)據(jù)匯總表Table 1 Summary of risk factor

2 風(fēng)險(xiǎn)因子的選取

2.1 共線性檢查

為了使數(shù)據(jù)更容易解釋和比較，對所有候選的風(fēng)險(xiǎn)因素進(jìn)行了標(biāo)準(zhǔn)化。隨后，在這些統(tǒng)計(jì)變量中，使用方差膨脹因子(VIF)結(jié)合斯皮爾曼相關(guān)系數(shù)進(jìn)行篩選。VIF越大，多重共線性越嚴(yán)重，通常以VIF<10為閾值，再綜合斯皮爾曼相關(guān)系數(shù)查看相關(guān)性>0.7的因子，刪除兩個(gè)相關(guān)性較高的因子的其中一個(gè)。最后，經(jīng)過篩選和比較保留了除降水和地面峰值加速度之外的14個(gè)風(fēng)險(xiǎn)因子。

2.2 隨機(jī)森林進(jìn)一步篩選風(fēng)險(xiǎn)因子

由于將貢獻(xiàn)率較低的風(fēng)險(xiǎn)因子納入回歸模型中將會(huì)影響其他相對重要的變量，所以使用隨機(jī)森林這種流行且有效的方法進(jìn)一步篩選對滑坡發(fā)生具有高度貢獻(xiàn)的風(fēng)險(xiǎn)因子[18]。該方法不依賴于特定的模型假設(shè)，而是基于數(shù)據(jù)驅(qū)動(dòng)的閾值來做出決策[19]。具體來說，就是要訓(xùn)練隨機(jī)森林模型，再使用隨機(jī)森林模型計(jì)算各個(gè)風(fēng)險(xiǎn)因子重要性得分。具有高重要性分?jǐn)?shù)的風(fēng)險(xiǎn)因子是結(jié)果的驅(qū)動(dòng)因素，它們的值對結(jié)果值有重大的影響。根據(jù)重要性得分對風(fēng)險(xiǎn)因子進(jìn)行排序，選擇排名靠前的11個(gè)因子作為滑坡是否發(fā)生的潛在預(yù)測因子(圖2)。這些篩選后的風(fēng)險(xiǎn)因子能夠構(gòu)建更簡單、更快速的預(yù)測模型。

圖2 隨機(jī)森林篩選結(jié)果Fig.2 Random forest screening results

3 滑坡易發(fā)性評價(jià)模型及結(jié)果分析

3.1 邏輯回歸模型

邏輯回歸公式可用于根據(jù)各種風(fēng)險(xiǎn)因素對二元結(jié)果進(jìn)行建模，對本研究來說體現(xiàn)為是否發(fā)生滑坡。它可以處理有多個(gè)自變量的情況，這對于模擬滑坡與各種風(fēng)險(xiǎn)因素之間的復(fù)雜關(guān)系非常重要[20]，已經(jīng)在滑坡易發(fā)性評價(jià)領(lǐng)域得到了廣泛的應(yīng)用[21，22]。邏輯回歸通過將滑坡事件的對數(shù)概率作為多個(gè)自變量的線性組合來對滑坡事件發(fā)生的概率進(jìn)行建模，具體計(jì)算公式如下：

logit(P)=β0+β1X1+β2X2+…+βnXn

(1)

(2)

其中，P是滑坡事件發(fā)生的幾率，logit(P)是滑坡事件發(fā)生的對數(shù)幾率，β0是截距項(xiàng)，β1，β2，…，βn是要確定的回歸系數(shù)，X1，X2，…，Xn是自變量。

3.2 貝葉斯空間邏輯回歸模型

1970年，沃爾多·托布勒(Waldo Tobler)提出了地理學(xué)第一定律(空間自相關(guān)定律)：“一切事物都與其他事物相關(guān)，但進(jìn)處的事物比遠(yuǎn)處的更相關(guān)[23]?！庇捎诘乩韺W(xué)第一定律的存在，使用空間回歸模型，考慮數(shù)據(jù)的自相關(guān)性來對滑坡這種具有地理空間位置信息的數(shù)據(jù)進(jìn)行建模是要比普通的回歸模型更有道理的[24]。

空間邏輯回歸是一種將邏輯回歸與空間分析技術(shù)相結(jié)合的統(tǒng)計(jì)方法，它用于對二元因變量與一個(gè)或多個(gè)自變量之間的關(guān)系建模，同時(shí)考慮數(shù)據(jù)的屬性信息和空間結(jié)構(gòu)。在空間邏輯回歸中，因變量是二元的，這意味著它只有兩種可能的結(jié)果(滑坡或者不滑坡)，自變量可以是連續(xù)的或是分類的。簡單來說，貝葉斯空間邏輯回歸模型就是在普通邏輯回歸模型之上考慮地理數(shù)據(jù)空間分布的結(jié)構(gòu)信息，并采用貝葉斯統(tǒng)計(jì)的思想完成的適用于地理數(shù)據(jù)分析和預(yù)測的模型。該模型在考慮了數(shù)據(jù)的空間自相關(guān)性后，估計(jì)自變量的系數(shù)及其與因變量的關(guān)系，具體的計(jì)算公式如下：

Y=β0+β1X1+β2X2+…+βnXn+ρWY+ε

(3)

其中，Y是滑坡事件發(fā)生的對數(shù)幾率，β0是截距項(xiàng)，β1，β2，…，βn是要確定的回歸系數(shù)，X1，X2，…，Xn是自變量，ρ是空間自相關(guān)參數(shù)，W是n×n維的空間權(quán)重矩陣，其中n是滑坡單元的總數(shù)，該矩陣定義了滑坡單元之間的鄰接關(guān)系。ρWY是由空間自相關(guān)引起的空間結(jié)構(gòu)因素。ε是服從高斯分布的誤差項(xiàng)。

3.3 評價(jià)結(jié)果及分析

在進(jìn)行滑坡易發(fā)性評價(jià)之前，首先要將研究區(qū)劃分為一個(gè)一個(gè)小的評價(jià)單元，以便于將風(fēng)險(xiǎn)指標(biāo)進(jìn)行收集與整理。目前主流的劃分單元有兩個(gè)：基于網(wǎng)格的單元和基于斜坡的單元，與基于網(wǎng)格的單元相比，斜坡單元可以貼切的反映研究區(qū)域的地形，更具有地貌意義[25]。本研究是采用改進(jìn)后的曲率流域法將研究區(qū)劃分為5 352個(gè)斜坡單元再進(jìn)行滑坡易發(fā)性評價(jià)[26]。

為了避免過度擬合，把數(shù)據(jù)集劃分為訓(xùn)練集(80%)和測試集(20%)，在訓(xùn)練集上訓(xùn)練邏輯回歸和空間邏輯回歸模型，并在測試集上評估其性能。再將模型得出的滑坡發(fā)生概率分別使用自然斷點(diǎn)法進(jìn)行分類，圖3是邏輯回歸輸出的滑坡易發(fā)性圖，圖4是空間邏輯回歸輸出的結(jié)果。區(qū)別于傳統(tǒng)的邏輯回歸，空間邏輯回歸還可以輸出一個(gè)能夠反應(yīng)研究區(qū)空間自相關(guān)效應(yīng)強(qiáng)弱的空間效應(yīng)圖(圖5)，這是更具有空間自相關(guān)規(guī)律(平滑后的)的滑坡易發(fā)性風(fēng)險(xiǎn)趨勢?？梢杂^察到，空間效應(yīng)圖與空間邏輯回歸輸出的滑坡易發(fā)性圖具有相似的空間格局，說明空間自相關(guān)效應(yīng)提示了滑坡發(fā)生的風(fēng)險(xiǎn)?；赂咭装l(fā)區(qū)和極高易發(fā)區(qū)均集中在研究區(qū)中部，位于龍門山斷裂帶，沿鹽井—五龍斷裂、雙石—大川斷裂發(fā)育。這些滑坡易發(fā)區(qū)處在青藏高原與四川盆地之間高地勢起伏區(qū)，海拔高的地區(qū)往往具有陡坡的特點(diǎn)，這會(huì)使土壤和巖石不穩(wěn)定，更容易發(fā)生滑坡。起伏的地形會(huì)產(chǎn)生壓力脊和凹陷，從而削弱土壤和巖石并使其更容易滑動(dòng)。

圖3 邏輯回歸易發(fā)性圖Fig.3 Logistic regression susceptibility map

圖4 空間邏輯回歸易發(fā)性圖Fig.4 Spatial logistic regression susceptibility map

圖5 空間邏輯回歸得出的空間效應(yīng)圖Fig.5 Spatial effect map obtained from spatial logistic regression

4 精度驗(yàn)證與討論

本研究采用ROC曲線和混淆矩陣評估兩個(gè)模型的性能。ROC曲線是二元分類模型性能的圖形表示，它顯示了靈敏度(真陽性率，TPR)和特異性(1-假陽性率，F(xiàn)PR)之間的權(quán)衡。通過測試每個(gè)可能的閾值并將每個(gè)結(jié)果繪制為曲線上的一個(gè)點(diǎn)來生成曲線，該曲線在 y 軸上繪制了 TPR，在 x 軸上繪制了 FPR(圖6)。閾值是模型預(yù)測正類的概率，越靠近左上角的ROC曲線表示模型性能越好。ROC曲線下的面積 (AUC)是衡量模型在所有可能的分類閾值下的整體性能的指標(biāo)。AUC的范圍從 0到1，完美模型的AUC值為1，0.7～0.8是可以接受的識(shí)別能力，0.8～0.9是具有優(yōu)秀的分類能力，模型的AUC值大于0.9表示模型具有杰出的分類識(shí)別能力[27]。混淆矩陣是顯示分類模型的真正值、真負(fù)值、假正值和假負(fù)值的表格，它評估每個(gè)類的正確和錯(cuò)誤分類實(shí)例的數(shù)量(表2)。

表2 混淆矩陣Table 2 Confusion matrix

圖6 邏輯回歸與空間邏輯回歸ROC曲線圖Fig.6 ROC curves of logistic regression and spatial logistic regression

如圖6所示，空間邏輯回歸模型的AUC值相較于傳統(tǒng)的邏輯回歸模型提升了近14%，大于了0.9，表示此模型具有杰出的分類識(shí)別能力。對比來看，空間邏輯回歸模型的ROC曲線也更靠近左上角，說明對于存在空間自相關(guān)性的滑坡數(shù)據(jù)，空間邏輯回歸模型在模型的預(yù)測效果方面顯著優(yōu)于傳統(tǒng)的邏輯回歸模型。結(jié)合混淆矩陣，空間邏輯回歸的預(yù)測準(zhǔn)確率為84.9%，在傳統(tǒng)邏輯回歸的基礎(chǔ)上提升了11.5%，單獨(dú)比較陽性(滑坡)預(yù)測正確率和陰性(非滑坡)預(yù)測正確率也分別提升了14.3%和9.8%。為了更直觀的進(jìn)行比較，我們利用源數(shù)據(jù)生成了實(shí)際滑坡點(diǎn)密度圖(圖7)，空間邏輯回歸的預(yù)測結(jié)果(圖4)與其高度一致，而傳統(tǒng)邏輯回歸在研究區(qū)東北邊緣將沒有發(fā)生滑坡的區(qū)域預(yù)測為高和極高易發(fā)性，在研究區(qū)西南有滑坡發(fā)生的區(qū)域預(yù)測為低和極低易發(fā)性。出現(xiàn)這種差異的原因還是在于有沒有考慮空間自相關(guān)這個(gè)不可忽視的風(fēng)險(xiǎn)因子，空間自相關(guān)效應(yīng)(圖5)會(huì)導(dǎo)致傳統(tǒng)的回歸模型高估或低估某些區(qū)域的滑坡易發(fā)性。因?yàn)猷徑膮^(qū)域往往具有相似的特征，如果這些特征與滑坡的發(fā)生相關(guān)聯(lián)，則模型可能預(yù)測附近區(qū)域也容易發(fā)生滑坡，即使它們在傳統(tǒng)回歸模型的預(yù)測中并不容易發(fā)生滑坡。

圖7 實(shí)際滑坡點(diǎn)密度圖Fig.7 Actual landslide point density map

5 結(jié)論

總的來說，空間邏輯回歸模型具有在斜坡單元水平上定義的空間潛在效應(yīng)，允許評估模型中的影響因子仍然無法解釋的空間影響，大大提高了傳統(tǒng)邏輯回歸模型的性能和預(yù)測精度，易于繪制出更準(zhǔn)確的滑坡易發(fā)性地圖，從而幫助政府確定高滑坡風(fēng)險(xiǎn)區(qū)域，并因地制宜的采取適當(dāng)措施降低風(fēng)險(xiǎn)，這有助于保護(hù)人民的生命和財(cái)產(chǎn)，減少滑坡對社會(huì)的整體影響。