基于統計分析和數值模擬方法的滑坡視摩擦系數影響因素及特征研究*

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

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

0 引 言

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

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

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

1 研究方法和材料

1.1 實例統計

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

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

表1 地震滑坡

1.2 數值模擬

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

(1)

(2)

(3)

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

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

圖2 滑坡模型

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

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

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

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

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

2 結果分析

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

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

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

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

2.2 影響因素顯著性分析

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

表2 影響因素顯著性分析

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

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

2.3 數值模擬與統計結果綜合分析

采用穩健線性回歸方法對V、θ、H與H/L的關系進行擬合，該方法相比傳統的最小二乘線性擬合方法，能夠有效消除統計中的異常點對回歸模型的影響(Naranjo et al.，1993)。數值模擬時，綜合統計結果中各變量的分布特征及確?；赂呕Ｐ偷暮侠硇裕兞繛閂時，θ=30°，H=1500m； 變量為θ時，V=5×105m3，H=500m； 變量為H時，V=2×106m3，θ=30°，將模擬結果與統計結果對比驗證分析。

從圖4a可以看出統計結果中隨著V的增大，藏東南區滑坡以及地震型滑坡的H/L值具有減小的趨勢，H/L與V的關系可用冪律關系H/L=αVβ表示(Lucas et al.，2007)。結果顯示：藏東南滑坡的α=1.96，β=-0.075，R2=0.189，地震滑坡的α=0.88，β=0.048，R2=0.075?？梢姴煌h境、不同成因下的滑坡H/L值與V的冪律關系及相關性的差異性明顯。模擬結果顯示：

H/L=1.544V-0.089

(4)

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

H/L=0.387 tanθ+0.105

(5)

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

圖4 H/L值與影響因素關系圖

H/L=0.214H+0.116

(6)

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

2.4 滑坡平均運動速度分析

滑坡的運動速度變化能直觀反應滑坡的沖擊力的變化情況。模擬結果見圖5、圖6。圖5顯示了當θ=30°、V=5×105m3、H=500m時不同時刻的滑坡速度云圖，圖6顯示了不同因素下滑坡整體的平均運動速度變化情況。從滑坡的速度云圖可知，在斜坡段滑坡前緣運動速度明顯高于后緣，前緣呈舌狀下滑，同時滑坡會在橫向上對稱擴散，直到運動到坡底，此時滑坡前緣速度達到最大。經過坡腳的轉折點后，滑坡速度開始迅速降低，并在坡底呈扇形堆積與擴散，此時仍停留在物源區的物質在有了臨空條件后開始陸續往坡腳運動，滑坡最終呈扇形堆積于坡腳。由圖6可知，V、θ、H明顯增大時，滑坡的平均運動速度的最大值也明顯增大?；略谶\動過程中，其平均運動速度先快速增大，達到峰值后迅速減小，隨后平穩直到停止運動。滑坡運動經過坡腳(圖6a中的B、C、D、E點； 圖6b中的C、D、E點； 圖6c中的A、B、C、D點)，其運動狀態由斜向下滑動，瞬間轉變為水平向的運動，運動狀態的突變以及滑坡體與水平面的碰撞需耗散大量能量，因而其平均運動速度開始迅速降低。而在θ和V較小時，滑坡物質的擴散也意味著滑坡與斜坡段的接觸面積增大，滑坡體受斜坡的摩擦阻力作用顯著的同時，缺少高值θ和高值V提供的速度增益，因而在斜坡段滑坡的平均速度便開始降低，最終滑坡運動到坡腳時速度已經很小(圖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點)，兩者速度接近； 貴州關嶺大寨崩滑流(劉傳正， 2010)，其體積7×105m3，主滑段斜坡角度35°，最大滑動高差150m，該滑坡最大平均運動速度12.4m·s-1，在模擬中變量H=250m，V=5×105m3，θ=30°時的最大平均速度12.3m·s-1(圖6c中A點)，兩者基本一致。這表明在宏觀參數基本一致的條件下，概化的滑坡模型其運動速度特征對實際滑坡仍具有一定的參考作用及適用性。

3 討 論

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

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

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

4 結 論

藏東南區滑坡及地震滑坡主要為中-大型滑坡，滑坡規模主要在105m3

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

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