王 冉 王學良 袁鴻鵠 孫娟娟 劉海洋 張如滿

(①中國科學院地質(zhì)與地球物理研究所， 頁巖氣與地質(zhì)工程院重點實驗室， 北京 100029， 中國)(②中國科學院地球科學研究院， 北京 100029， 中國)(③中國科學院大學， 地球與行星科學學院， 北京 100049， 中國)(④北京市水利規(guī)劃設(shè)計研究院， 北京 100048， 中國)

0 引 言

滑坡的最大垂直運動距離H和最大水平運動距離L是滑坡運動學研究的重要參數(shù)。前者能決定滑坡的總能量，后者是滑坡災害防災減災重點考慮的因素。H/L被稱為視摩擦系數(shù)，其能夠衡量滑坡的運動能力(Johnson et al.，2017)。不少學者研究了坡腳角度、滑面長度、偏轉(zhuǎn)角等對滑坡運動的影響，詹威威等(2017)選取了38處汶川地震中觸發(fā)的地震滑坡通過逐步回歸方法建立了滑坡-碎屑流最大水平運動距離L的最優(yōu)多元回歸模型； 楊海龍等(2019)利用90處地震滑坡，基于主要地形參數(shù)和滑坡規(guī)模建立了最大水平運動距離的預測統(tǒng)計模型； 鄭光等(2019)考慮滑坡碎屑流運動距離與勢能的關(guān)系，得到遠程運動距離L和勢能VH間的方程； 汪發(fā)武等(2002)分析了土粒破碎特性在高速滑坡運動過程的作用，結(jié)合視摩擦系數(shù)模型和滑坡運動量方程及連續(xù)性方程，對滑坡的運動范圍進行了預測。

隨著遙感技術(shù)的不斷發(fā)展，高精度DEM影像、無人機航測等為進行崩滑幾何特征參數(shù)的識別和獲取提供了便利。王學良等(2018)采用航空遙感、無人機航測技術(shù)，結(jié)合巖體結(jié)構(gòu)分析與數(shù)值模擬，進行了山區(qū)輸變電快速識別分析崩塌(滾石)； Pradhan et al. (2016)將QuickBird遙感影像與通過LiDAR獲得的DEM數(shù)據(jù)融合構(gòu)建滑坡識別規(guī)則，滑坡識別率達到90%以上。融合遙感技術(shù)的滑坡數(shù)據(jù)獲取彌補了傳統(tǒng)調(diào)查手段無法快速獲取滑坡詳細信息的不足，使批量化的滑坡信息獲取成為可能(劉春玲等， 2010； Scaioni et al.，2014； 王珊珊等， 2015)。

本文利用機載雷達獲取的精細DEM(10m分辨率)識別出藏東南通麥和魯朗區(qū)的大量滑坡，將其與地震滑坡統(tǒng)計對比分析。分析3種因素(滑坡體積V、斜坡角度θ、滑坡最大垂直運動距離H)以及H/L值的水平，并采用單因素方差分析法分析各因素對H/L值影響的顯著性； 進一步建立滑坡的概化模型，采用數(shù)值模擬綜合分析各因素與視摩擦系數(shù)(H/L)的相關(guān)關(guān)系及對H/L的影響特征，同時對滑坡運動中平均速度變化特征做了簡要分析，以期為研究滑坡運動特征提供數(shù)據(jù)參考。

1 研究方法和材料

1.1 實例統(tǒng)計

結(jié)合前期工作(王學良等， 2018)，作者對藏東南區(qū)的滑坡進行識別，共識別出滑坡189處(圖1)。提取滑坡的剖面，利用AutoCAD獲取滑坡幾何信息，假定滑坡物源為柱體，V通過物源區(qū)的側(cè)面積與其寬度相乘的方法獲得，H通過滑前后緣高程減去滑后前緣高程獲得，L為滑前后緣與滑后前緣連線的水平分量，θ通過在ArcGIS中計算出滑坡范圍內(nèi)各柵格的坡度值并取其平均值獲得。提取文獻數(shù)據(jù)(表1)，獲得229處地震型滑坡的相關(guān)數(shù)據(jù)，其中包括2008年中國汶川地震引起的滑坡133處、日本歷史地震造成的地震滑坡77處以及國外其他地震滑坡19處。

圖1 藏東南魯朗—通麥區(qū)解譯滑坡分布

表1 地震滑坡

1.2 數(shù)值模擬

在將基于連續(xù)介質(zhì)力學原理的連續(xù)模型應用于滑坡的模擬時，考慮滑坡運動中，相比于運動距離，其深度很小，因此深度積分模型成為模擬計算滑坡運動的有效工具(宋宜祥， 2018)。Savage et al. (1991)的研究成果將深度積分模型引入到地球科學和工程計算。Ouyang et al. (2017a， 2017b)對經(jīng)典Navier-Stokes方程逐步重新推導，采用深度積分方法，用MacCormack-TVD有限差分方法求解，開發(fā)出山地災害動力模擬軟件Massflow，并對多例滑坡運動特征進行數(shù)值模擬，結(jié)果與實際觀測結(jié)果吻合。Massflow中單元體物質(zhì)必須滿足質(zhì)量和動量守恒方程，基本方程如下：

(1)

(2)

(3)

式中：h，u，v分別為運動體的高度，x和y方向上的速度；gx，gy，gz分別為傾斜坐標系下重力加速度在3個方向上的分量；kap為側(cè)向土壓力系數(shù);τzx，τzy為基底阻力。

采用控制變量法在Massflow中模擬V、θ、H取不同值時對H/L值的影響特征，并分析滑坡運動過程的速度變化，滑坡概化模型見圖2。數(shù)值模擬方案如下：

圖2 滑坡模型

(1)建立地形文件。在ArcGIS中構(gòu)造地形的等值線，再生成TIN文件，接著將TIN文件轉(zhuǎn)換成TIF文件，地形的柵格文件建立完成。

(2)建立物源文件。將2個具有高差的地形TIF文件相減并在ArcGIS中進行裁剪獲得物源的TIF文件。

(3)轉(zhuǎn)換文件格式。將地形TIF文件和物源TIF文件轉(zhuǎn)換成Massflow程序中可執(zhí)行的TXT文件。

(4)運行程序。設(shè)置模擬參數(shù)，并分組記錄數(shù)據(jù)。

Massflow中提供滑坡模擬的庫倫模型，模型中的參數(shù)較少且都有具體的物理意義。庫倫模型只需要3個參數(shù)，分別是黏聚力，內(nèi)摩擦角以及孔隙水壓力系數(shù)。建立的滑坡模型是為重點分析滑坡運動過程中的共性特征，進而在選擇參數(shù)時參考已有案例的參數(shù)范圍取適當值(Ouyang et al.，2017a，2017b; 王學良等， 2019)，最終黏聚力取值20kPa，基底摩擦系數(shù)取值0.34，將物源作干碎屑狀態(tài)處理，孔隙水壓力系數(shù)取值0，邊界條件為墻體邊界，物質(zhì)密度取值2500kg·m-3。

2 結(jié)果分析

2.1 影響因素及H/L值水平分析

按照小型滑坡、中型滑坡、大型滑坡、特大型滑坡以及體積的連續(xù)量級區(qū)間(Qi et al.，2010)，將V設(shè)置為5個區(qū)間；θ、H、H/L的區(qū)間劃分根據(jù)變量取值及統(tǒng)計數(shù)據(jù)分布特征，將θ劃分為5個區(qū)間段，H和H/L劃分為6個區(qū)間段(圖3)。統(tǒng)計各水平區(qū)間數(shù)據(jù)量的占比，為便于描述分析，定義占比最大的水平區(qū)間為優(yōu)勢水平，占比之和大于80%的水平區(qū)間為主要水平。統(tǒng)計結(jié)果顯示，藏東南區(qū)滑坡和地震型滑坡的V和θ的優(yōu)勢水平和主要水平差別并不顯著，其中V的優(yōu)勢水平均為106m3

圖3 影響因素及H/L值水平

綜合分析發(fā)現(xiàn)，藏東南區(qū)滑坡、地震型滑坡以中-大型滑坡為主(陸盟等， 2020)，θ主要集中在20°<θ<40°范圍內(nèi)，相較于藏東南區(qū)的滑坡，地震型滑坡顯著差異表現(xiàn)為低H值下，H/L值普遍小于0.6，以遠程滑坡為主(張明等，2010)。

2.2 影響因素顯著性分析

基于單因素試驗法的方差分析確定V、θ、H對H/L值影響的顯著性水平。分析結(jié)果見表2，

表2 影響因素顯著性分析

在實際應用中一般取置信區(qū)間α=0.05和α=0.01，并且當F≤F0.05時，認為影響不顯著； 當F0.05F0.01時，則認為影響特別顯著。對于藏東南區(qū)滑坡，F(xiàn)θ=139.17>F0.01=4.61，F(xiàn)V=17.07>F0.01=4.61，因此θ和V對H/L值的影響特別顯著，是決定性因素且θ>V，而FH=1.35F0.01=4.61，F(xiàn)H=5.29>F0.01=4.61，因此θ和H對H/L值的影響特別顯著且θ>H，而FV=2.60

綜合分析發(fā)現(xiàn)，影響因素中θ對H/L值的影響特別顯著，θ是控制H/L值的決定性因素。數(shù)理統(tǒng)計方法在滑坡行為預測甚至地質(zhì)災害預測中具有重要作用(Middleton et al.，1994)，而在對滑坡的H/L值進行預測分析時θ具有重要的參考意義。

2.3 數(shù)值模擬與統(tǒng)計結(jié)果綜合分析

采用穩(wěn)健線性回歸方法對V、θ、H與H/L的關(guān)系進行擬合，該方法相比傳統(tǒng)的最小二乘線性擬合方法，能夠有效消除統(tǒng)計中的異常點對回歸模型的影響(Naranjo et al.，1993)。數(shù)值模擬時，綜合統(tǒng)計結(jié)果中各變量的分布特征及確保滑坡概化模型的合理性，變量為V時，θ=30°，H=1500m； 變量為θ時，V=5×105m3，H=500m； 變量為H時，V=2×106m3，θ=30°，將模擬結(jié)果與統(tǒng)計結(jié)果對比驗證分析。

從圖4a可以看出統(tǒng)計結(jié)果中隨著V的增大，藏東南區(qū)滑坡以及地震型滑坡的H/L值具有減小的趨勢，H/L與V的關(guān)系可用冪律關(guān)系H/L=αVβ表示(Lucas et al.，2007)。結(jié)果顯示：藏東南滑坡的α=1.96，β=-0.075，R2=0.189，地震滑坡的α=0.88，β=0.048，R2=0.075。可見不同環(huán)境、不同成因下的滑坡H/L值與V的冪律關(guān)系及相關(guān)性的差異性明顯。模擬結(jié)果顯示：

H/L=1.544V-0.089

(4)

從圖4b可以看出統(tǒng)計結(jié)果中隨著θ的增大，藏東南區(qū)滑坡和地震型滑坡的H/L值具有明顯線性增大的趨勢，關(guān)系表達式分別為：H/L=0.96tanθ+0.03(R2=0.84)，H/L=0.513tanθ+0.145(R2=0.51)。H/L值與tanθ具有較好的線性正相關(guān)關(guān)系。模擬結(jié)果顯示：

H/L=0.387 tanθ+0.105

(5)

從圖4c可知，藏東南區(qū)滑坡及地震型滑坡的H/L與H無線性關(guān)系。將藏東南滑坡和地震型滑坡數(shù)據(jù)綜合線性擬合顯示R2=0.17，線性關(guān)系增強，說明在較大H范圍內(nèi)，H/L與H仍具有一定線性正相關(guān)關(guān)系。模擬結(jié)果顯示：

圖4 H/L值與影響因素關(guān)系圖

H/L=0.214H+0.116

(6)

對比發(fā)現(xiàn)，θ與H/L值的線性相關(guān)關(guān)系最顯著，方差分析也同時表明H/L值伴隨θ的變化具有特別顯著的變化，分析認為H/L值受控于θ明顯，利用單因素分析預測H/L值時，選擇θ可能具有更小的誤差。在各因素與H/L的相關(guān)關(guān)系上，藏東南區(qū)滑坡普遍比地震型滑坡更好，這表明地震作用增加了各因素對H/L值影響特征的不確定性，離散了V、θ、H與H/L的相關(guān)關(guān)系。各因素在相同水平下，地震型滑坡的H/L值普遍更低，在較小的H下，地震型滑坡卻能運動更遠的距離，這表明地震型滑坡具有遠程的特征。模擬結(jié)果中H/L與各影響因素的相關(guān)關(guān)系及變化趨勢與統(tǒng)計結(jié)果吻合較好，但同時H/L值普遍小于統(tǒng)計中的H/L值且低于遠程滑坡的閾值。

2.4 滑坡平均運動速度分析

滑坡的運動速度變化能直觀反應滑坡的沖擊力的變化情況。模擬結(jié)果見圖5、圖6。圖5顯示了當θ=30°、V=5×105m3、H=500m時不同時刻的滑坡速度云圖，圖6顯示了不同因素下滑坡整體的平均運動速度變化情況。從滑坡的速度云圖可知，在斜坡段滑坡前緣運動速度明顯高于后緣，前緣呈舌狀下滑，同時滑坡會在橫向上對稱擴散，直到運動到坡底，此時滑坡前緣速度達到最大。經(jīng)過坡腳的轉(zhuǎn)折點后，滑坡速度開始迅速降低，并在坡底呈扇形堆積與擴散，此時仍停留在物源區(qū)的物質(zhì)在有了臨空條件后開始陸續(xù)往坡腳運動，滑坡最終呈扇形堆積于坡腳。由圖6可知，V、θ、H明顯增大時，滑坡的平均運動速度的最大值也明顯增大?；略谶\動過程中，其平均運動速度先快速增大，達到峰值后迅速減小，隨后平穩(wěn)直到停止運動?；逻\動經(jīng)過坡腳(圖6a中的B、C、D、E點； 圖6b中的C、D、E點； 圖6c中的A、B、C、D點)，其運動狀態(tài)由斜向下滑動，瞬間轉(zhuǎn)變?yōu)樗较虻倪\動，運動狀態(tài)的突變以及滑坡體與水平面的碰撞需耗散大量能量，因而其平均運動速度開始迅速降低。而在θ和V較小時，滑坡物質(zhì)的擴散也意味著滑坡與斜坡段的接觸面積增大，滑坡體受斜坡的摩擦阻力作用顯著的同時，缺少高值θ和高值V提供的速度增益，因而在斜坡段滑坡的平均速度便開始降低，最終滑坡運動到坡腳時速度已經(jīng)很小(圖6a中的A點； 圖6b中的A、B點)。

圖5 滑坡速度云圖

圖6 滑坡平均速度變化曲線

云南頭寨溝滑坡(胡凱等，2016)，其體積1.8×107m3，滑坡垂直下滑高度1000m，平均斜坡角度35°，該滑坡最大平均運動速度28.1m·s-1，在滑坡模擬中變量V=3.5×107m3，θ=30°，H=1500m時的最大平均運動速度為27.5m·s-1(圖6a中D點)，兩者速度基本一致； 都江堰廟壩滑坡(茍富剛等，2012)，其體積1.13×106m3，山體的平均坡度35°，滑坡垂直下滑高度440m，該滑坡最大平均運動速度為20.1m·s-1，在模擬中變量H=500m，V=2×106m3，θ=40°時的最大平均速度19.4m·s-1(圖6b中D點)，兩者速度接近； 貴州關(guān)嶺大寨崩滑流(劉傳正， 2010)，其體積7×105m3，主滑段斜坡角度35°，最大滑動高差150m，該滑坡最大平均運動速度12.4m·s-1，在模擬中變量H=250m，V=5×105m3，θ=30°時的最大平均速度12.3m·s-1(圖6c中A點)，兩者基本一致。這表明在宏觀參數(shù)基本一致的條件下，概化的滑坡模型其運動速度特征對實際滑坡仍具有一定的參考作用及適用性。

3 討 論

藏東南是地震多發(fā)的地區(qū)，統(tǒng)計數(shù)據(jù)中仍有37%的滑坡表現(xiàn)出遠程特征(H/L<0.6)，這在一定程度上佐證了藏東南區(qū)地震可能引發(fā)了部分滑坡。地震型滑坡在較低H的情況下，卻能運動更遠的距離，從而使H/L值相比藏東南區(qū)滑坡的H/L值更低。在地震力作用下，滑坡體更容易被地震力快速拋出，并進一步與下伏巖體撞擊破碎，破碎的巖體碎屑加劇了滑坡的流態(tài)化特點，從而使地震滑坡表現(xiàn)出遠程的特點(殷躍平， 2008； 吳樹仁等， 2010)。

本文所采用的單因素分析和方差分析，所需要的參數(shù)能較方便測量得到，且可以獨立地分析各個因素對H/L值的影響特點，不用先驗的假設(shè)特定的模型或者因素間的關(guān)系，不足之處在于未考慮到因素間是否存在相關(guān)關(guān)系。文中采用滑坡的概化模型嘗試分析滑坡的時空演化共性特征，因而數(shù)值模擬參數(shù)并未取自具體滑坡，但在應用于特定滑坡分析時，參數(shù)的確定仍需根據(jù)反演、現(xiàn)場測試等手段獲得。后續(xù)工作還應就參數(shù)值的改變對滑坡運移特征的影響做進一步研究。

數(shù)值模擬結(jié)果與統(tǒng)計值比較發(fā)現(xiàn)，在同一因素水平條件下，模擬結(jié)果中H/L值均低于遠程滑坡H/L閾值(H/L=0.6)(張明等， 2010)。這表明采用連續(xù)介質(zhì)力學方法對解釋遠程滑坡的運動行為可能具有更好的適應性。值得注意的是真實滑坡的運移過程受眾多因素影響，如滑坡體性質(zhì)、地形條件、誘發(fā)因素等(Crosta et al.，2018)，而滑坡的概化模型并不考慮這諸多因素，所以并不能說明基于連續(xù)介質(zhì)力學法的滑坡數(shù)值模擬僅適用于遠程滑坡。由于控制滑坡運動的不確定因素很多，如降雨、坡面粗糙度、巖性等(Wang et al.，2014)，H/L與各因素之間更深層次的物理機制存在很大的爭議，它們之間的定量關(guān)系還有待更進一步的探索研究。

4 結(jié) 論

藏東南區(qū)滑坡及地震滑坡主要為中-大型滑坡，滑坡規(guī)模主要在105m3

相較于V和H，藏東南區(qū)滑坡和地震滑坡的H/L值受θ控制作用顯著，H/L值與θ具有良好的線性關(guān)系。地震作用會弱化θ、V、H與H/L值的相關(guān)程度，地震誘發(fā)的滑坡遠程特征顯著。

采用連續(xù)介質(zhì)力學方法模擬滑坡的運移特征，結(jié)果表明滑坡的運動速度變化主要分為3個階段，速度陡增階段、速度峰值后陡降階段和速度平緩至停止階段，坡腳的轉(zhuǎn)折段能顯著降低滑坡的運動速度。