多元逐步回归模型在滑坡变形预测分析中的应用

莫中平，顾功开，伍中华

(中国三峡建设管理有限公司，云南 昆明 651500)

在各类滑坡体的监测数据分析中，为了解滑坡的变形规律，把握其发展趋势，识别出影响边坡变形的主要因素，研究边坡在何种工况下可能出现异常情况，并找出对其进行预测和预警的方法，往往需要利用较长时期的安全监测数据建立监控模型。

滑坡变形分析和预测中，较多采用突变理论、模糊数学、小波分析等方法研究边坡变形，或者用灰色理论构建预测模型等。但以上方法一般更多关注边坡变形观测数据的变化趋势，较少研究环境变量如降雨及温度等因素与滑坡变形之间的关系。在长期实测变形数据的基础上，考虑相应的环境变量，建立数据模型分析其与滑坡变形的相关性，研究其对滑坡变形的作用大小，预测滑坡的变形趋势，成为滑坡变形数据分析的重要方法[1]。

本文对西南地区某滑坡体的变形监测资料进行了分析，纳入时效、温度、降雨因子，应用SPSS统计分析软件构建了滑坡变形的多元逐步回归模型，分析该滑坡的变形规律，探讨了各个影响因子对滑坡变形的作用。

1 回归模型的构成

影响边坡稳定性的因素众多，有内因，有也外因。地表水和地下水的作用、地震、风化作用、温度变化、开挖施工及工程荷载等是外因；边坡的地形地貌、地质构造、岩土体结构、滑体材料重度以及力学性质等是内因，它们通常起主要的控制作用。本文主要讨论降雨和气温对滑坡变形的影响，其他难以界定的影响因素归为时效作用，所有的分析主要基于这3个因素展开和进行。因此，滑坡体的位移统计回归模型[2]表示为

δ=δθ+δT+δU

(1)

式中：δ为位移量；δθ为时效分量；δT为温度分量；δU为降雨分量。

1.1 时效分量δθ

滑坡变形的时效分量构成复杂，它包含了边坡滑动产生的不可逆位移和自身形变，还包含岩土体的塑性变形、徐变等。根据文献[3]，时效分量的表达式为

δθ=β1(θ-θ0)+β2(lnθ-lnθ0)

(2)

式中：β1、β2为时效因子回归系数；θ为观测日至始测日的累计天数t/100；θ0是建模数据集第一个测值日至始测日的天数t0/100。

1.2 温度分量δT

温度变化会影响边坡裂隙的开合度和岩土应力，可能会对边坡稳定产生一定的影响。因为气温一般表现为周期性变化，为简化计算，通常用正余弦函数的周期项来模拟温度[4]，将温度变化产生的变形分量表示为

(3)

式中，t为观测日至始测日的累计天数。

1.3 降雨分量δU

考虑到滑坡体在雨季的变形量大于旱季，在模型中引入降雨分量。由于降雨入渗抬高地下水位，改变滑坡体含水量，从而影响边坡稳定。降雨对边坡变形的影响会有一定时间的滞后，一般在15 d内[5]。考虑到本滑坡规模巨大，滑坡体厚达30～70 m，入渗时间更长，因此增加相应的滞后雨量因子：选择观测日前1～60 d的降雨量，再引入观测日前5～15 d、观测日前16～30 d的降雨量作为因子。降雨引起的变形分量可表示为式(4)

(4)

式中，Um为前m天的平均降雨量，m=1，2，3，5，10，15，30，60，5～15，16～30。

1.4 确定模型表达式

综上所述，将式(2)至(4)代入式(1)，滑坡体变形的统计模型表示为

(5)

式中，β0为常数项，表示未被自变量解释的部分。

2 应用SPSS实现多元逐步回归算法

多元逐步回归模型用来揭示被解释变量与其他多个解释变量之间的线性关系，其经验回归方程是：

y=β0+β1x1+β2x2+……+βpxp

(6)

式中：y为被解释变量；β0为模型常数项；βp为回归系数；xp为解释变量。

逐步回归计算是建立在F检验的基础上，将显著因子逐个引进回归方程。当回归方程中引入一个因子后，由于因子间具有相关性，可能使已在方程中的其他因子变得不再显著，此时应从回归方程中剔除该不显著因子。因此在接纳一个因子后，应对已在回归方程中的全部因子进行F检验，逐一剔除不显著的因子，再继续对下一个因子用F检验来决定是否引入回归方程(一次仅引入一个)。通过反复运用F检验，逐一剔除和引入，直到得到所需的最佳回归方程[4]。

常用统计分析软件有：

1)Excel。Excel是常用的数据处理软件，其自带的“回归分析”工具，可以使用“最小二乘法”对数据进行线性回归分析，但不能进行逐步回归分析。

2)SAS。是大型集成应用软件系统，数据处理功能强大，但纯编程界面，兼容的数据格式少，操作难度大，适用统计专业研究人员。

3)SPSS。是著名的统计分析软件之一。软件统计分析功能较齐全，绘制表格、图形较方便，输出结果比较直观。人机对话界面友好，操作方便，只要对统计分析原理有一定了解就可以使用。

本文使用SPSS进行统计分析，依次选择[分析]-[回归]-[线性]-[步进模式]，应用其内置的逐步回归分析模块，得到最佳的回归方程。SPSS可以根据需要输出相应的回归系数、拟合值、残差等数据及各类图形，揭示影响因变量的主要成分，进而对滑坡变形规律进行分析。

3 工程实例

3.1 工程概况

西南地区某滑坡体位于金沙江岸，呈弧形，滑坡区域平面面积4.5 km2。前缘高程800～820 m，延伸至金沙江面，后缘高程1 900～2 200 m，前后缘高差超过1 000 m。滑坡体地形复杂，根据滑坡地质构造、地形地貌和成因划分为5个区域。其中Ⅱ区变形较大，呈牵引式蠕滑变形。

Ⅱ区滑体厚度30～70 m不等，滑带位于千枚岩碎屑土内，滑带厚度2.0～9.0 m。滑带土为紫红、灰黑色粉质粘土夹砾石，土为紫红色粉质粘土，硬～可塑状，结构紧密，具明显挤压错动特征，可见光面及擦痕；砾石、碎石岩性为紫红色千枚岩，大小一般0.2 cm×0.5 cm～1 cm×2 cm，多呈次棱角～次圆状、圆状，少量棱角状，含量20%～40%。

2004年开始对该滑坡体实施安全监测。监测项目有地表水平位移、垂直沉降、深部侧向位移、地下水位、渗流监测等。在变形较大的Ⅱ区共布置水平位移监测点16个。

该滑坡体地处亚热带、热带高原性季风气候区，冬夏温差不大。所处干热河谷地区，在气候上明显存在干季和雨季之分，降水量年内分配不均：干季(11～4月)雨水稀少,降水量仅占全年降水量的4%～13%；雨季(5～10)降水量可占全年降水量的87%～96%。同时降水年际分布不均，自2012年至2019年，年平均降雨量586.6 mm，最大年份降雨779.6 mm，最小年份337.3 mm。降雨量时间分布见图1。

图1 滑坡体区域降雨量-时间分布图

3.2 模型应用及分析

该滑坡体的观测频次：5月至10月间每月观测2次，11月到次年5月每月观测1次。选取2013年1月至2020年2月的实测变形、降雨等资料，对Ⅱ区共16个变形测点建立回归模型。16个测点回归模型复相关系数R都大于0.96，精度较高。限于篇幅，表1只列出部分水平位移监测点的回归系数、标准化系数，复相关系数(R)等。

3.2.1 影响因子分析

1)时效分量δθ。表1的结果显示，在各测点的回归模型中都包含了时效因子。根据标准化系数分析，时效因子β1、β2的标准化系数之和均大于0.97，说明时效因素引起的变形占总变形量的绝大部分。回归分析结果表明时效是影响滑坡体表面变形的主要因素。

2)温度分量δT。表1的结果显示，温度因子β3～β6入选模型的只占总测点的1/4，且其标准化系数之和均小于0.06，说明温度变化对滑坡变形影响很小。

3)降雨分量δU。从实测变形的规律看，大部分测点在雨季的变形量要略大于旱季。但通过模型计算，降雨因子的标准化系数之和均小于0.05，说明降雨对变形的影响并不显著。

3.2.2 模型预测

根据实测数据，对滑坡的变形进行预测监控。以TP01为例，代入回归系数，其统计模型如下：

(7)

将时效和温度分量代入式(7)，2019年10月～2020年2月的实测值与模型拟合值见表2。从拟合值与实测值对比可知，模型计算的结果与实测值吻合良好，残差很小。图2为水平位移测点TP01的实测值、拟合值及残差过程线，图3为回归标准化残差。

3.2.3 降雨影响分析

模型显示，降雨因子并非该滑坡体变形的显著因素。根据平面滑动的边坡稳定性计算方法[6]，边坡滑动稳定性的计算公式为

(8)

R=[(G+Gb)cosθ-Qsinθ-Vsinθ-U]tanφ+cL

(9)

T=(G+Gb)sinθ+Qcosθ+Vcosθ

(10)

式中：T为滑体单位宽度重力及其他外力引起的下滑力，kN/m；R为滑体单位宽度重力及其他外力引起的抗滑力，kN/m；c为滑面的黏聚力，kPa；φ为滑面的内摩擦角，°；θ为滑面倾角，°；G为滑体单位宽度自重，kN/m。

表1 滑坡体部分表面变形监测点位移逐步回归模型系数表

表2 滑坡体表面变形监测点TP01实测位移与模型拟合数据 mm

从边坡滑动稳定性计算公式看，影响滑动稳定的主要因素是滑体自重G、滑动面土层的黏聚力c、内摩擦角φ，以及滑面单位宽度总水压力U和后缘陡倾裂隙面上的单位宽度总水压力V。对本滑坡体的相关参数进行分析：

1)为监测滑坡体下地水，共布置14个水位观测孔。监测结果表明，Ⅱ区地下水位总是在滑带之上(见图4)。因为该滑坡体厚达30～70 m，且主要为松散土体，雨水无法大量渗透到滑动面，因此难以对滑动面土层的黏聚力c、内摩擦角φ产生实质影响。即使少量雨水下渗，因为地下水总高于滑带，也不会对c、φ产生实质改变。

图2 滑坡体表面变形监测点TP01实测值、拟合值、残差图

图3 滑坡体表面变形监测点TP01回归标准化残差直方图

图4 滑坡体Ⅱ区地下水位线与滑带位置图

监测成果还表明多数观测孔在降雨前后水位没有明显变化，说明降雨不会导致滑面总水压力U和后缘陡倾裂隙面上的总水压力V明显增大。

2)滑坡体区域自2013年至2019年，最大月降雨量为272 mm，换算单位面积重量为2.72 kN/m2；滑坡体厚度30～70 m，其土层重度按18 kN/m3计算，单位面积土体自重为540～1 260 kN/m2，降雨给滑体单位面积自重G增加的重量仅为0.22%～0.5%，影响甚小。

综上所述，降雨对控制滑动稳定的各个参数影响较小，滑坡体下滑力T及抗滑力R在降雨前后变化很小，说明降雨对变形速率影响不大。这与模型揭示的降雨因子并非该滑坡体变形显著因素的结论一致，说明模型是可信的。

4 结 语

本文选取时效、温度和降雨作为滑坡体变形的影响因子，应用多元逐步回归算法，建立滑坡变形的统计回归模型，较好拟合了时效、温度和降雨量对滑坡变形的影响。

西南某滑坡体变形监测成果的多元逐步回归分析结果表明：该滑坡体的变形主要受时效因子影响。由于该滑坡体地下水位总高于滑带，且滑体厚度大以致降雨入渗困难，降雨对控制滑动稳定的自重G、黏聚力c、内摩擦角等参数影响较小，因此降雨对该滑坡体变形速率影响不大。

应用SPSS软件对滑坡体变形数据进行多元逐步回归分析，计算较为简便，拟合精度较高，可以及时、准确了解变形发展趋势，掌握边坡安全状态，识别影响边坡变形的主要因素，对滑坡体变形趋势进行预测，具有重要的实用价值。