基于統(tǒng)計(jì)分析和數(shù)值模擬方法的滑坡視摩擦系數(shù)影響因素及特征研究*

王 冉 王學(xué)良 袁鴻鵠 孫娟娟 劉海洋 張如滿

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

0 引 言

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

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

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

1 研究方法和材料

1.1 實(shí)例統(tǒng)計(jì)

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

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

表1 地震滑坡

1.2 數(shù)值模擬

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

(1)

(2)

(3)

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

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

圖2 滑坡模型

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

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

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

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

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

2 結(jié)果分析

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

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

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

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

2.2 影響因素顯著性分析

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

表2 影響因素顯著性分析

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

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

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

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

從圖4a可以看出統(tǒng)計(jì)結(jié)果中隨著V的增大，藏東南區(qū)滑坡以及地震型滑坡的H/L值具有減小的趨勢(shì)，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?？梢姴煌h(huán)境、不同成因下的滑坡H/L值與V的冪律關(guān)系及相關(guān)性的差異性明顯。模擬結(jié)果顯示：

H/L=1.544V-0.089

(4)

從圖4b可以看出統(tǒng)計(jì)結(jié)果中隨著θ的增大，藏東南區(qū)滑坡和地震型滑坡的H/L值具有明顯線性增大的趨勢(shì)，關(guān)系表達(dá)式分別為：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)系增強(qiáng)，說明在較大H范圍內(nèi)，H/L與H仍具有一定線性正相關(guān)關(guān)系。模擬結(jié)果顯示：

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

H/L=0.214H+0.116

(6)

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

2.4 滑坡平均運(yùn)動(dòng)速度分析

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

圖5 滑坡速度云圖

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

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

3 討 論

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

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

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

4 結(jié) 論

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

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

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