在滑坡的发育-失稳演化过程中，缓慢变形往往先于滑动破坏。对滑坡区域地表的变形进行监测可以揭示滑坡详细的运动学信息，进而通过分析滑坡变形的时空特征，为滑坡风险评价、触发分析和发生机理等研究提供科学支撑。大地测量新技术凭借其大范围、小尺度、高精度的技术优势已成为滑坡时空变形特征研究的重要途径。现有的大地测量技术已能完全满足在时间域和空间域精确监测滑坡变形演化过程的需求^[1-5]，但这些研究对监测数据的利用仍不充分，缺少对时空形变场特征的深入挖掘。不同类型的滑坡在其演化过程中的变形模式可以是拉伸、压缩、剪切、旋转等多种基础变形的综合^[6]，滑坡变形过程中的地表时空形变场是揭示滑坡运动学机理和实现量化风险评价的重要依据^[7]。本文以山西省晋城市沁水县里必煤矿附近一处滑坡体为例，通过建立小变形和有限变形条件下的应变模型，整合分布滑坡体各处离散测站上的位移变化信息，以解算里必滑坡体不同演化阶段的应变分布，揭示该滑坡变形整体与局部的时空特征，为滑坡运动学研究及定量风险评价提供支撑。

1 应变模型 1.1 小变形条件下的应变模型

地球曲率产生的系统误差在小区域内可以忽略不计，因此针对滑坡的应变计算可基于平面基准^[8]。当测站间距不大时，变形体内部的位移分量可以写成其位置坐标的函数，并可在其邻域内的任意一点展开为泰勒级数：

$ \left\{\begin{array}{l} u_x^0=u_x+\frac{\partial u_x}{\partial x} \Delta x+\frac{\partial u_x}{\partial y} \Delta y \\ v_y^0=v_y+\frac{\partial v_x}{\partial x} \Delta x+\frac{\partial v_x}{\partial y} \Delta y \end{array}\right. $

(1)

式中，Δx=x₀－x，Δy=y₀－y。

据计算，里必滑坡体在其滑前和滑后阶段，测区内两两测站之间距离变化量的最大值仍小于1/34 000。在这种小变形的条件下，物体受力发生变形的过程中其内部任意点上的绝对位置变化也可以表达为某一点及其邻域的整体平移(u_x, u_y)、刚性旋转ω和应变(ε_x, ε_y, γ_xy)3种基本变化的总和，具体模型与求解可见参考文献[9]。

1.2 有限变形条件下的应变模型

设二维平面上的变形区域内存在距离相近的两点O和P，在变形前有：

$ \overrightarrow {OP} = \mathit{\boldsymbol{r}}(\Delta x, \Delta y) $

(2)

式中，(Δx, Δy)为P点关于O点的坐标差。设$\overrightarrow {OP}$关于x轴和y轴的方向余弦分别为l、m，则有Δx=l|r|，Δy=m|r|，代入式(1)则有：

$ \left\{\begin{array}{l} u(P)=u(O)+\frac{\partial u_x}{\partial x} l|\boldsymbol{r}|+\frac{\partial u_x}{\partial y} m|\boldsymbol{r}| \\ v(P)=v(O)+\frac{\partial v_x}{\partial x} l|\boldsymbol{r}|+\frac{\partial v_x}{\partial y} m|\boldsymbol{r}| \end{array}\right. $

(3)

令$\overrightarrow {OP'} $为变形后的$\overrightarrow {OP} $，则有：

$ \begin{gathered} \left|\overrightarrow{O P^{\prime}}\right|=[\Delta x+u(P)-u(O)]^2+ \\ {[\Delta y+v(P)-v(O)]^2} \end{gathered} $

(4)

由伸长度的定义可得：$\left|\overrightarrow{O P^{\prime}}\right|=(1+e)|\boldsymbol{r}|$，代入式(4)整理可得：

$ \begin{gathered} (1+e)^2=\left[\left(1+\frac{\partial u}{\partial x}\right) l+\frac{\partial u}{\partial y} m\right]^2+ \\ {\left[\frac{\partial v}{\partial x} l+\left(1+\frac{\partial v}{\partial y}\right) m\right]^2} \end{gathered} $

(5)

当变形不再满足微小变形假设时，变形过程中某点位移分量的偏微分和材料伸长度的高次项不能再视为微小量而被忽略，因此将式(5)展开可得：

$ \begin{gathered} e+\frac{1}{2} e^2=\left\{\frac{\partial u}{\partial x}+\frac{1}{2}\left[\left(\frac{\partial u}{\partial x}\right)^2+\left(\frac{\partial v}{\partial x}\right)^2\right]\right\} l^2+ \\ {\left[\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}+\left(\frac{\partial u}{\partial x} \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \frac{\partial v}{\partial y}\right)\right] l m+} \\ \left\{\frac{\partial v}{\partial y}+\frac{1}{2}\left[\left(\frac{\partial u}{\partial y}\right)^2+\left(\frac{\partial v}{\partial y}\right)^2\right]\right\} m^2 \end{gathered} $

(6)

根据文献[10]给出的有限变形拉格朗日方法，平面上二维基本应变张量参数与位移分量偏微分的关系可写为：

$ \left\{\begin{array}{l} \varepsilon_x=\frac{\partial u}{\partial x}+\frac{1}{2}\left[\left(\frac{\partial u}{\partial x}\right)^2+\left(\frac{\partial v}{\partial x}\right)^2\right] \\ \varepsilon_y=\frac{\partial v}{\partial y}+\frac{1}{2}\left[\left(\frac{\partial u}{\partial y}\right)^2+\left(\frac{\partial v}{\partial y}\right)^2\right] \\ \gamma_{x y}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}+\left(\frac{\partial u}{\partial x} \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \frac{\partial v}{\partial y}\right) \end{array}\right. $

(7)

代入式(6)得：

$ e+\frac{1}{2} e^2=\varepsilon_x l^2+\gamma_{x y} l m+\varepsilon_y m^2 $

(8)

2 滑坡体时空变形特征分析 2.1 区域概况

里必滑坡体海拔为756.64~796.52 m，区域地形起伏较大但总体面积较小，滑坡体最大宽度约550 m，轴向长度约230 m。区内地层主要由第四系全新统人工填土、黄土状土、上更新统粉质粘土、含碎石粉质粘土、块石及二叠系上石盒子组砂质泥岩等组成，其中各类土层中孔隙发育程度较高且夹杂有大量砾石。土层中间夹有多层薄层强风化泥岩，风化程度、裂隙发育程度及岩体破碎程度都较高。滑坡变形监测采用高精度全站仪，精度等级为二级。如图 1所示，变形监测网沿主滑动轴方向布设了5个近似平行的观测剖面，测站间距约50 m，所有测站观测标志均安装在强制观测桩上，观测桩埋深于大于厚度0.8 m的冻土层。

在小变形条件和有限变形条件下求解应变张量参数，需要利用各测站在相邻两期观测之间的位移速率矢量构建观测方程。本文以单位时间内各测站的坐标变化量为待定参数，通过以速率为参数的动态平差方法^[11]求得各测站相邻两期之间的位移速率，结果如图 2所示。从图中可以看出：1)2020-10~12是滑坡变形的显著阶段，x轴分量波动强于y轴分量，表明滑坡体的优势变形方向接近NS向；2)2021-04-16~22滑坡体中后部发生滑动失稳，各观测剖面03、04号测站变形显著，x轴分量大于y轴分量，表明滑坡体的失稳滑动方向更接近x轴方向；3)滑坡体发生失稳破坏之前经历了一段相对稳定的准静态阶段，各测站位移速率几乎为0，这与2020-12对滑坡体前缘支挡墙进行锚索加固有关。

已知各测站不同时段的位移速率可以计算得到相应时段内任意两测站之间平面距离的相对变化量，这是在有限变形条件下应变模型求解时所需的基本观测量。求解应变张量参数时需要至少3个不同的实测点提供起算数据，本文以各测站为顶点建立Delaunay三角形应变单元，并采用最小二乘模型求解基本应变张量参数，结果如图 3~8所示。

2.2 滑前演化阶段

目标滑坡体滑前演化阶段的时空变形过程并不均匀。时域上，存在随时间推移而趋于稳定的总体趋势；空间上，不同部位在不同时期的变形强度、变形特征均有明显的差异变化。监测初期(2020-10~11)，滑坡体的最显著变形特征为前缘中部沿主滑动轴方向的局部滑动变形，其中最大位移发生在PC01测站，日平均位移速率达0.6 mm/d，位移速率自前缘中心向两侧和后缘衰减。滑坡体两侧边界测站的位移速率及方向都与中心测站有显著差异，随时间推移，滑坡体各部份位移均发生明显的衰减。2020-12前缘位移几乎不可见，转而向后缘发生显著扩散，这一变化特征与12月初完成支挡墙的锚索加固相对应。2021年以来，滑坡体仅在边缘位置的少数测站发生滑动，滑坡体整体处于准静止状态。图 3中绿色扇形为非对称旋转张量，表示各单元的转动变形趋势，可以看出，滑坡体各部分之间的相对差异旋转并不显著，异常主要发生在滑坡体的边界位置，且总体表现为逆时针的差异转动。

图 4为当前阶段滑动区域的最大、最小主应变分布，相较于位移矢量的时空分布，最大、最小主应变分布可以更加明确地突出监测区域局部变形的空间差异。以2020-12-03为节点，滑坡体的演化变形可以分为2个阶段，前期滑坡变形相对更强且主要发生在滑坡体前缘与支挡墙相邻的部位，最大压应变值达－49.3 μstrain。2020-10两期滑坡前缘分布的两处强压应变方向不一致，表明滑坡体的局部变形受到支挡墙走向的显著影响，同时滑坡体中后部以不均匀的弱张应变为主要特征，表明其处于不稳定状态。2020-11滑坡体前缘主压应变轴与主滑动轴方向一致程度较高，但是11月两期应变分布有明显变化。支挡墙被锚索加固前后(2020-11-19~12-03)，滑坡体前缘PD01位置的应变强度大幅增强，主压应变值由－49.2 μstrain突变为－82.3 μstrain，由此可以估算出，在锚索加固之前支挡墙大约只吸收了40%滑坡体前缘滑动的冲击。值得注意的是，2020-12-13滑坡体后缘830 m等高线位置突然出现强拉张变形，主张应变轴方向与监测初期前缘主压应变方向一致。2020-12-10野外勘测发现该张应变位置新出现一条10~20 cm宽的裂缝，裂缝走向NE~NNE，大致与主张应变轴方向垂直。随后的监测结果表明，该裂缝位置恰好是滑坡体后缘发生滑动失稳的边界。

2.3 滑动失稳阶段

滑坡体在2021-04-16~22发生滑动失稳，日均位移速率最大值达到4.5 mm/d。大致以830 m等高线为界，滑坡体滑动可以划分为前后2个变形区，2个分区水平位移趋势相反，本次滑坡为后缘滑动。加固后的支挡墙对滑坡体前缘滑动有十分显著的限制作用，滑动后滑坡体前缘在支挡墙的阻碍作用下发生反向逆推，形成如图 5所示的位移分布。

图 6中最大、最小主应变强度的空间分布与图 5给出的位移空间分布相似，即以830 m等高线为界，可以将滑坡体划分为2个不同强度的变形区，且分界存在明显的强度突变。滑坡体后缘应变强度更高，以压缩变形为主要特征，横向上存在明显的空间趋势变化。滑坡体南北两侧主压应变近似垂直于邻近的滑动边界，中部以近似滑动轴方向的压缩变形为主要特征，同时具有显著的横向拉张变形，表明滑坡体滑动过程中伴随有显著的横向扩张。应变极值出现在滑坡体的南部边界PA02测站后方，变形以垂直于滑坡体边界的压缩变形为主，最大应变率超过1 450 μstrain。滑坡体前缘各部应变强度十分微弱，表明期间该区域内部相对变形较为微弱，滑坡体前缘与支挡墙之间的相互挤压作用不显著。这是由于本文计算所得的应变和位移速率实际是2021-04-16~22的平均值，其中必然包含部分滑前和滑后分量，其结果不能实时反映滑坡体滑动失稳的真实动态形变过程。

2.4 滑后调整阶段

如图 7所示，滑坡发生后3个月水平位移速率与最大、最小主应变的时空变化差异巨大，表明期间滑坡体一直处于震荡调整过程，尚未达到稳定状态，不具备明显的统一趋势性滑动。当前阶段滑坡变形最为复杂，变形模式也同时包含平移、旋转和应变3种基本类型，特别是2021-05-28滑坡体中部广泛存在顺时针转动变形。该转动趋势随时间衰减，其中前缘和中部衰减快、后缘衰减慢，至2021-07后缘仍然存在明显的转动变形。结合图 8中最大、最小主应变的时空分布特征可知，滑坡前中部受2021-04失稳滑动的影响有限，在滑后阶段调整至平稳的速度更快。后缘由于裂缝发育的阻隔及失稳滑动带来的强变形效应，其调整至平稳的过程相对缓慢。2021-07-24滑坡体后缘中部水平位移出现局部统一的趋势，同时滑坡体前缘相对变形基本稳定，但后缘仍然具有显著的相对变形。主要变形特征以近轴向的挤压和近横向的拉张变形为主，相对变形剧烈区域主要集中在中南部，北部边缘相对平静。

3 讨论

应变张量参数求解方面，式(1)和式(7)均至少需要3个已知量构建观测值向量，才可按照最小二乘模型求解。小变形条件下的应变张量参数求解需要变形监测点的位移矢量，而有限变形条件下求解应变张量参数则需要进一步计算两测站之间距离的相对变化量。因此，应变计算时宜以测区内各测站为顶点组建基本图形单元进行参数求解，在该前提下求得的应变张量参数应该是相应时段内代表各图形单元变形特征的平均值。显然，应变单元面积越小，变形监测的频率越高，所得的结果越精确。

滑动与失稳的触发方面，根据沁水县气象站观测资料，每年6月至9月为雨季，且近年来9月多现极端降雨天气，本文监测到滑坡体在2020-10初出现较强的局部滑动。降水对滑坡的触发作用通常存在一定时间的滞后性，因此里必滑坡体在2020-10出现的显著滑动应该与当年雨季，特别是9月份的极端降雨天气密切相关。滑坡体的后缘在2021-04-16~22发生了失稳滑动，滑坡发生前，滑坡体基本处于准静止状态，其整体位移和内部相对变形几乎为0。考虑到研究区域每年10月中旬至次年3月上旬为霜冻期，最大冻土厚度可超过80 cm，因此本次的后缘滑坡与冻土解冻效应有较高的时间关联性。

应变张量能够定量描述滑坡体不同部位、不同时期的内部相对变形强度，进而有效揭示滑坡体内部的相对变形细节。在滑坡体的滑前演化阶段和滑动失稳阶段，水平位移矢量的局部差异不明显，滑坡体各部多具有相近的位移趋势，而通过应变计算可以深入挖掘这类区域之间的相对变形差异。特别是锚索加固之前支挡墙仍然对滑坡体前缘有一定的限制作用，可以吸收近40%的前缘滑动，这是位移矢量的时空分布特征难以体现的。另一方面，应变的时空分布特征对于滑坡体的局部相对突变高度敏感，在变形特征的定性分析上也优于传统的位移矢量分析，如830 m等高线附近在2020-12-13出现沿滑动轴的强局部拉张，该异常特征难以由位移矢量或时间序列变化获取。经过实地勘察对比发现，该局部拉张与一条10~20 cm宽的横向裂缝发育关系密切，是滑坡体结构性破坏的有力证据，这点亦被后期滑动失稳阶段和滑后调整阶段的形变场空间分布特征所验证。

4 结语

本文构建了小变形条件和有限变形条件下的应变模型，并与位移矢量分布相结合，研究了里必滑坡体滑前演化阶段、滑动失稳阶段和滑后调整阶段的相对形变时空分布特征。结果显示，应变场能够更加直观地定量描述滑坡体不同部位、不同演化阶段形变场的时空分布特征。

本文应变分析结果显示，滑前演化阶段滑坡体前缘与支挡墙之间存在较强的挤压冲击，并且其作用强度和趋势的时空分布并不均匀，不同时段强压缩变形存在于滑坡体与支挡墙之间的不同位置。此外，尽管位移矢量自后向前存在显著的空间加速特征，但应变场显示支挡墙对于滑坡的限制作用一直存在。通过对比锚索加固后的滑坡变形强度差异可知，前期支挡墙对滑坡大概起到了40%左右的限制作用。相较于位移速率分析，应变分析可直观且定量地描述滑坡体局部异常变形的属性特征，如2020-12-13出现在滑坡体830 m等高线位置的轴向强拉张，该应变张量很好地对应了区域内一条10~20 cm宽的裂缝。因此，局部异常的强应变集中对监测对象的结构性破坏有很高的灵敏度。张性裂缝的发育一直都是滑坡发育-失稳过程中的重要标识，在本文滑坡体准稳定、小变形背景下出现这样的局部异常，表明滑坡体的破坏性发育仍在持续。