基于统计分析和数值模拟方法的滑坡视摩擦系数影响因素及特征研究*

2021-07-19王学良袁鸿鹄孙娟娟刘海洋张如满

工程地质学报 2021年3期

王冉王学良袁鸿鹄孙娟娟刘海洋张如满

(①中国科学院地质与地球物理研究所，页岩气与地质工程院重点实验室，北京 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)。数值模拟时，综合统计结果中各变量的分布特征及确保滑坡概化模型的合理性，变量为V时，θ=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/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与各因素之间更深层次的物理机制存在很大的争议，它们之间的定量关系还有待更进一步的探索研究。