基于Logistic回歸和MCMC方法評價地震滑坡敏感性

賀 倩, 汪 明, 劉 凱

(1.北京師范大學 地表過程與資源生態(tài)國家重點實驗室, 北京 100875;2.北京師范大學 減災與應急管理研究院 地理科學學部, 北京 100875;3.北京師范大學 環(huán)境演變與自然災害教育部國家重點實驗室, 北京 100875)

地震滑坡是最常見的次生災害之一，許多學者對地震滑坡的敏感性，危害和風險進行了研究[1]。在山區(qū)，地震引發(fā)的滑坡會掩埋村莊，阻塞道路和河流，阻礙救災，其造成的損失可能比地震本身造成的損失更加嚴重[2]。因此，準確的評估地震滑坡的敏感性具有重要意義。我國西南地區(qū)地震多發(fā)，地震誘發(fā)的滑坡造成了交通阻塞，生命線中斷，嚴重威脅著人民的生命財產(chǎn)安全[3]。

在滑坡敏感性評估中，Logistic回歸是廣泛應用的模型之一[1,4-5]。邏輯回歸模型的參數(shù)估計以及參數(shù)不確定性的定量分析對于模型的構(gòu)建來說十分重要。經(jīng)典Logistic回歸模型一般利用最大似然法來估計方程中的參數(shù)，最大似然是點估計，其對參數(shù)的不確定性估計較為困難[6-7]。馬爾可夫鏈蒙特卡洛方法是在貝葉斯理論框架下，將蒙特卡洛模擬和馬爾科夫鏈模型結(jié)合，實現(xiàn)抽樣樣本隨模擬的進行而改變的動態(tài)模擬[8]，其可以獲得參數(shù)估計的后驗分布，從而利用參數(shù)后驗分布來評估參數(shù)的不確定性。將MCMC方法應用于Logistic回歸的參數(shù)估計能夠結(jié)合馬爾科夫模型長期預測的優(yōu)勢和Logistic模型模擬復雜系統(tǒng)變量和響應變量變化的能力，具有科學性和實用性[9]。MCMC方法在水文模型參數(shù)不確定性分析中應用廣泛[10-12]，而在滑坡敏感性建模上的研究較少。此外，隨著計算機技術的快速發(fā)展，基于MCMC方法估計Logistic模型中的回歸系數(shù)易于實現(xiàn)[12-13]，特別是R語言的應用使得MCMC方法實際應用更加簡便可行。

本文以2017-8-8九寨溝地震，2014-8-3魯?shù)榈卣鸷?013-4-20蘆山地震為案例進行研究。本文從地震動參數(shù)，地形以及地質(zhì)三方面構(gòu)建滑坡敏感性評價指標體系?；贛CMC和Logistic回歸方法建立地震滑坡敏感性的評價模型，對模型的精度進行評價，然后分別繪制三次地震事件的滑坡敏感性圖。通過MCMC獲得的參數(shù)后驗概率分布來量化參數(shù)估計值的不確定性。最后對區(qū)域影響因子的重要性進行評價。

1 研究區(qū)與數(shù)據(jù)源

1.1 研究區(qū)概況

中國西南地區(qū)地處地震帶，該地區(qū)地殼活動頻繁，地震多發(fā)，再加之地形起伏大，由地震引發(fā)的滑坡災害嚴重。本文選擇了位于該區(qū)域的三次地震事件作為研究案例。九寨溝地震發(fā)生于2017年8月8日，位于九寨溝國家公園，震級為7級，震源深度20 km，地震引發(fā)了各種類型的地質(zhì)災害，包括滑坡、潰壩和大壩毀壞等，共造成25人遇難[2]。魯?shù)榭h位于云南省西北部，為中國地震最活躍的地區(qū)之一[14]。2014年8月3日，云南省昭通市魯?shù)榭h發(fā)生6.5級地震，造成了617人死亡[15]，引發(fā)了大量滑坡[16]。蘆山位于四川平原與青藏高原的交界處，橫跨龍門山構(gòu)造帶和四川盆地兩個構(gòu)造單元。2013年4月20日，四川雅安蘆山縣發(fā)生7.0級地震，震中位于青藏高原最東緣，龍門山斷裂的西南部，地震造成了217人死亡，誘發(fā)了數(shù)以千計的滑坡，造成了嚴重的破壞并嚴重阻礙了救援行動[17-18]。

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

1.2.1 滑坡編錄數(shù)據(jù) 滑坡編錄數(shù)據(jù)對于評估區(qū)域規(guī)模的滑坡災害或風險至關重要。本文的地震滑坡編錄數(shù)據(jù)包括：九寨溝地震誘發(fā)的226處滑坡[19]，魯?shù)榈卣鹫T發(fā)的滑坡1 024處[15]，蘆山地震滑坡15 546處[20]?；聰?shù)據(jù)是研究學者根據(jù)地震前后的高分辨率的衛(wèi)星影像進行系統(tǒng)準確的解譯得到的。本文所使用的滑坡數(shù)據(jù)類型為點。本文采用Zhu等[21]提出的雙緩沖區(qū)采樣方法，在內(nèi)部緩沖區(qū)(1 km半徑)和外部緩沖區(qū)(15 km半徑)之間的區(qū)域進行隨機采樣，分別獲得每個事件等量的非滑坡樣本點。然后將每個事件中滑坡和非滑坡點隨機分成70%的訓練樣本用于模型的構(gòu)建，30%的測試樣本用于模型精度的驗證。

1.2.2 指標選擇 參照前人的研究和數(shù)據(jù)的可用性[4, 22-25]，本文從地形，地質(zhì)和地震動參數(shù)3個方面共選擇了9個指標，包括高程，坡度，地形濕度指數(shù)(Topographic Wetness Index，TWI)，矢量穩(wěn)健性測度(Vector Ruggedness Measure，VRM)，平面曲率，剖面曲率，距離斷層的距離(簡稱斷層距)，巖性，修正麥卡利烈度(Modified Mercalli Intensity，MMI)。地形數(shù)據(jù)為DEM，來源于SRTM (Shuttle Radar Topography Mission) DEM (http:∥srtm.csi.cgiar.org/)，空間分辨率為90 m×90 m。利用DEM數(shù)據(jù)和SAGA GIS平臺提取得到研究區(qū)的高程、坡度、TWI，VRM、平面曲率和剖面曲率。地質(zhì)數(shù)據(jù)包括斷層和巖性數(shù)據(jù)。地震滑坡的空間分布與斷層之間存在著較強的相關性[26]。本文所使用的斷層數(shù)據(jù)來自GEM(Global Earthquake Model)全球活動斷層數(shù)據(jù)庫(the GEM Global Active Faults Database, GAF-DB)，GAF-DB匯總了許多地區(qū)的數(shù)據(jù)集，包含13 500個斷層，空間完備性高[27]。巖性數(shù)據(jù)來源于高分辨率的全球巖性圖(Global Lithology Map, GLiM)數(shù)據(jù)集[28]，該數(shù)據(jù)集的巖性類別分為3個級別，本文使用第一級的分類(16類)，包括松散沉積巖，碎屑沉積巖，火山碎屑巖等。地震動參數(shù)數(shù)據(jù)來自美國地質(zhì)調(diào)查局(USGS)的ShakeMap(https:∥earthquake.usgs.gov/data/shakemap/)，ShakeMap提供了地震后地面運動和震動強度的近實時地圖。本文使用ShakeMap中的修正麥卡利烈度產(chǎn)品。

2 方法與模型

2.1 Logistic回歸

Logistic回歸是分析滑坡敏感性的有效的方法[5,29]，其因變量是二分類變量(滑坡或者非滑坡)，解釋變量(自變量)可以是離散的類別型，也可以是連續(xù)的數(shù)值型。Logistic回歸的優(yōu)勢是變量不需要滿足正態(tài)分布[30]。Logistic回歸通過構(gòu)建二元變量和多個自變量之間的回歸關系來確定事件發(fā)生的概率，公式如下：

(1)

Z=β0+β1x1+β2x2+…+βkxk

(2)

式中:P表示單個柵格像素內(nèi)滑坡發(fā)生的概率；xk表示自變量；β0表示算法的截距(常數(shù)項)；βk表示不同滑坡因子的回歸系數(shù)；k為滑坡影響因子的個數(shù)。在本文中,Logistic回歸的因變量指定為1,0,分別表示滑坡和非滑坡。模型輸出在0～1之間表示滑坡發(fā)生的概率,即本文中的滑坡敏感性。

2.2 馬爾可夫鏈蒙特卡羅

MCMC使用隨機游動的Metropolis算法從Logistic模型的后驗分布中進行模擬。其通過因變量和自變量構(gòu)建的貝葉斯框架來解釋先驗信息。MCMC參數(shù)估計中包含3個關鍵部分：先驗分布、似然函數(shù)和后驗分布。Das等[7]建立了一個簡單的頻率Logistic模型的Bayes等價物，如下所示：

yi～Bernoulli(logit-1(η))

(3)

(4)

βj～N(0,0.00001)j=0,…,k

(5)

式中:yi表示因變量；βj是具有獨立正態(tài)先驗分布的系數(shù);xj表示第j個變量的值；η是線性預測變量。利用貝葉斯理論,可以獲得模型下參數(shù)的后驗分布,利用貝葉斯公式進行擴展,其后驗分布與先驗分布和似然函數(shù)的乘積成比例關系[7,9]。

MCMC結(jié)合參數(shù)的先驗分布通過迭代仿真得到參數(shù)估計值的后驗分布和概率密度函數(shù)。本文基于R語言中的MCMCpack包進行模型參數(shù)的估計[31]。在每個事件中，使用70%的樣本通過MCMClogit函數(shù)來擬合Logistic方程，其余30%來進行模型預測精度的驗證。

2.3 模型精度評價

混淆矩陣是分類問題中模型精度評價的重要方法[32]。在二分類問題中，根據(jù)真實類別以及模型預測結(jié)果的組合劃分為真正例(True positive，TP)，假正例(False positive，F(xiàn)P)，真反例(True negative，TN)以及假反例(False negative, FN)。在本研究中，正例為滑坡，反例為非滑坡，分別用1，0表示。TP和TN分別表示正確分為滑坡和非滑坡的數(shù)量，F(xiàn)P和FN分別表示錯誤地分為滑坡和非滑坡的數(shù)量?；诨煜仃囉嬎愠隹偩?Accuracy)，精確率(Precision)和召回率(Recall)。精確率(Precision)在滑坡二分類問題中指的是在所有滑坡樣本中預測正確的滑坡數(shù)比例，召回率(Recall)指的是預測為滑坡的樣本點中真正的滑坡數(shù)比例，F(xiàn)1-score是基于召回率和精確率的調(diào)和平均定義的，計算公式見(6)—(9)。

(6)

(7)

(8)

(9)

此外，本文還計算了受試者工作特征曲線(Receiver Operating Characteristic，ROC)和曲線下面積(Area Under the Curve, AUC)。ROC反映了特異性和敏感性之間的關系，AUC反映模型的總精度，AUC范圍為0.5～1，AUC值越大表明模型的預測能力越高。

3 結(jié)果與分析

3.1 指標共線性分析

指標之間的多重共線性會降低模型預測能力[33]。在建立模型前，本文利用容差(Tolerance，TOL)和方差膨脹因子(the variance inflation, VIF)[34-35]進行因子共線性的檢驗以確保因子之間相互獨立。通常認為，當VIF值大于5時，因子之間存在強共線性，應該將該變量從模型中剔除[36-37]。VIF值越大，變量與其他變量之間就具有越強的共線性。TOL小于0.2表示因子之間存在多重共線性，小于0.1，表明存在嚴重的共線性。多重共線性分析表明所選擇的指標因子的容差均在0.4以上，高于閾值0.1，并且其VIF值均小于3，低于閾值5(表1)。綜合分析表明本文所選的因子之間不存在多重共線性問題，驗證了評價指標的合理性。

表1 影響因子的多重共線性分析

3.2 指標與滑坡的關系分析

本文分析了滑坡數(shù)量占比與不同指標之間的關系。將各個連續(xù)的指標變量等分成30個區(qū)間，分析不同區(qū)間內(nèi)滑坡數(shù)量的分布情況。由于篇幅有限，本文在此處呈現(xiàn)蘆山地震滑坡的結(jié)果圖。蘆山的地震滑坡主要集中在高程值較小的區(qū)間上(圖1A)，表明此次地震事件中區(qū)域內(nèi)滑坡主要分布在低海拔地區(qū)。隨著坡度的增大，蘆山的滑坡比率與坡度呈現(xiàn)出正相關關系(圖1B)。對于平面曲率和剖面曲率，隨著曲率的增加，滑坡比率增加(圖1C—D)。在TWI和VRM較大時，滑坡占比相對更高(圖1E—F)，而在TWI和VRM較小時，滑坡數(shù)量較多。斷層距越大，滑坡占比越小(圖1G)，大部分滑坡點集中在斷層距較小的區(qū)間內(nèi)。地震參數(shù)MMI與滑坡比率呈現(xiàn)出明顯的正相關關系(圖1H)，隨著MMI的增大，滑坡占比整體升高。

圖1 蘆山地震，滑坡發(fā)生比率與各指標之間的關系

3.3 MCMC的模型參數(shù)估計

本文對馬爾可夫鏈的中burnin參數(shù)設置為2 000，表示棄用前2 000次的抽樣消除初始值的效益影響。mcmc參數(shù)設置為40 000，表示使用后面的40 000次抽樣進行后驗參數(shù)的估計。此外，采樣過程的收斂性檢驗至關重要，只有迭代收斂之后才能夠用來獲得參數(shù)的后驗分布。本文利用樣本的軌跡圖和樣本間的自相關性進行收斂性診斷。圖2展示了蘆山地震滑坡案例中的MMI和坡度參數(shù)估計的軌跡圖和樣本自相關圖。通過樣本路徑圖可以看出，采樣過程基本穩(wěn)定。從樣本的自相關圖可以看出最終樣本達到平穩(wěn)狀態(tài)，收斂性較好。綜合考慮樣本采樣軌跡圖和自相關圖，可以認為參數(shù)樣本的分布已收斂到參數(shù)的后驗分布。

圖2 蘆山地震中馬爾可夫鏈迭代軌跡圖和自相關圖

MCMC參數(shù)估計的結(jié)果見表2。在三次滑坡事件中，9個因子對于滑坡的發(fā)生都有顯著的影響(p<0.05)。 MCMC方法得到的是后驗分布的均值(Mean)作為參數(shù)的估計值。為了比較不同參數(shù)的不確定性，利用變異系數(shù)(Coefficient of variation， CV)進行分析，可以看出，在九寨溝地震事件中，巖性因子的參數(shù)不確定性最高，其次為VRM，而其他因子的CV值都較低。對于魯?shù)榈卣鹗录拭媲实牟淮_定性較高，其次為巖性和平面曲率因子。在蘆山地震案例中，所有因子的不確定性均較低，只有平面曲率因子的CV值絕對值大于1。各個指標因子參數(shù)在蘆山事件中的不確定性最小，分析認為這可能與樣本數(shù)量有關(蘆山地震事件中的滑坡樣本數(shù)量是其他兩個事件總和的12倍之多)。巖性因子在九寨溝和魯?shù)閮蓚€事件中的不確定性都相對較高，分析認為可能是巖性數(shù)據(jù)的分辨率相對較低造成的，而九寨溝地震和魯?shù)榈卣鸬膮^(qū)域面積小，一定程度上導致了其參數(shù)的不確定性較高。

表2 三次地震滑坡事件中基于MCMC的Logistic回歸模型中的參數(shù)估計

3.4 模型預測精度評價

在每個事件中，利用30%的未參與模型訓練的測試數(shù)據(jù)集來驗證模型的預測精度。結(jié)果表明Logistic回歸模型在3個區(qū)域都具有較高的預測精度。本文所使用的滑坡數(shù)據(jù)中，蘆山地震中的滑坡點最多，其次是魯?shù)?，九寨溝滑坡最少。九寨溝地震，魯?shù)榈卣鸷吞J山地震的測試總精度分別為0.875，0.928，0.927，三次滑坡事件中精確率均在0.9以上，召回率均在0.8以上，F(xiàn)1-score均大于0.85(表3)。從各個預測結(jié)果的混淆矩陣圖也可以看出(圖3)，預測正確的滑坡(TP)和非滑坡(TN)占比最大。此外，三次地震滑坡事件中的AUC值均大于0.9。綜合分析表明基于MCMC的Logistic回歸模型具有較好的預測精度。

表3 模型測試精度評價參數(shù)

注：A，D代表九寨溝地震；B，E代表魯?shù)榈卣?；C，F(xiàn)代表蘆山地震。

3.5 滑坡敏感性評價結(jié)果

利用得到的Logistic回歸模型對研究區(qū)的每個像元進行預測，得到滑坡敏感性評價圖(圖4)。將滑坡點疊加在敏感性圖上進行對比，通過對比發(fā)現(xiàn)3個事件中的地震滑坡點都基本上落在0.9～1的高敏感性區(qū)間內(nèi)，僅有少部分的滑坡點落在敏感性較低的區(qū)域。在蘆山地區(qū)，由于滑坡點數(shù)量多，分布密集，造成敏感性高的區(qū)域在空間上成片分布。

3.6 因子相對重要性評價

本文利用隨機森林算法分別對3個案例中的因子重要性進行分析，對比因子相對重要性圖(圖5)發(fā)現(xiàn)，地震參數(shù)MMI，高程和斷層距在三次地震事件中重要性排序均較高。在九寨溝和蘆山地區(qū)，MMI為最重要的因子，而在魯?shù)榈卣鸹轮?，高程為最重要的影響因子。巖性在蘆山地區(qū)的重要性相對較大，而在九寨溝和魯?shù)榈貐^(qū)重要性相對低，分析認為蘆山地震滑坡分布的空間范圍較大，區(qū)域內(nèi)巖性種類較多，巖性指標對模型預測的滑坡敏感性的空間分布有較大影響。而九寨溝地震和魯?shù)榈卣鸹路植嫉姆秶^小，區(qū)域內(nèi)巖性種類較為單一，對滑坡敏感性的空間分布影響不明顯。

4 結(jié) 論

本文采用基于MCMC方法的Logistic回歸模型研究地震滑坡敏感性，利用MCMC方法確定模型參數(shù)，并獲得參數(shù)的后驗分布，對參數(shù)的不確定性進行了評價，并基于所得到的模型進行了區(qū)域的滑坡敏感性空間制圖。在蘆山地震事件中，參數(shù)的不確定性均比較低，在九寨溝事件中巖性的不確定性較高，而在魯?shù)榈卣鸢咐校瑤r性，剖面曲率和平面曲率的不確定性較高。整體而言，三次事件中的大部分因子的回歸系數(shù)具有較低的不確定性。本文所構(gòu)建模型具有較高的預測精度，在三次地震滑坡事件中AUC值均大于0.9，表明了MCMC方法在Logistic回歸模型中進行參數(shù)估計的有效性?；掠绊懸蜃拥南鄬χ匾栽诓煌膮^(qū)域有差異，其中MMI，高程以及斷層距在3個案例的相對重要性都是最大的。

圖4 滑坡敏感性

圖5 因子重要性

本文基于MCMC方法的研究結(jié)果對于Logistic回歸模型在滑坡敏感性上的研究具有重要參考意義。本文所確定的滑坡敏感性因子包括了地形因素，地質(zhì)條件，和地震動參數(shù)，未來的研究中將會考慮更多的因素在滑坡敏感性建模中的應用以及采用更高分辨率的數(shù)據(jù)以獲得區(qū)域更為細致的滑坡敏感性分布圖。