复杂型滑坡，如变质岩区大型高位滑坡，规模较大、结构复杂，对其进行变形监测预报难度较大，相关研究也较少。以往的滑坡变形监测预报研究多偏向于现状变形特征分析^[1]，缺乏变形发展规律研究；同时，想要合理评价滑坡后续变形规律，多源信息综合评价是必要的。本文以变质岩区大型高位滑坡为工程背景，利用变形预测和Manner-Kendall分析法^[2](以下简称M-K分析法)联合开展滑坡后续变形的监测预报研究，以充分掌握变质岩区大型高位滑坡的理论基础，旨在为其防治提供更有力的支持。

1 分析模型

以某变质岩区大型高位滑坡变形监测成果为基础，进行既有变形特征分析和变形发展分析，前者主要是利用统计学方法分析滑坡累积变形及变形速率等特征；后者是利用预测模型和M-K分析法对比研究滑坡后续变形规律。

1.1 预测模型

不确定因素(仪器误差、温度变化等)会影响滑坡变形监测，即

$ x(t)=r(t)+\varepsilon(t) $

(1)

式中，x(t)为滑坡变形监测值，r(t)为滑坡变形真实值，ε(t)为不确定变形量。

因此，在滑坡变形预测模型的构建过程中，先利用自适应噪声完备集合经验模态分解(CEEMDAN)模型^[3]进行滑坡变形数据的分解处理，具体步骤如下。

1) 先以集合经验模态分解(EEMD)模型进行I次分解处理，得到平均状态下的第1个模态分量IMF₁^′(t)：

$ \operatorname{IMF}_1^{\prime}(t)=\frac{1}{I} \sum\limits_{i=1}^I \operatorname{IMF}_1^i(t) $

(2)

式中，IMF₁ⁱ(t)为第i次分解的模态分量。

2) 进一步计算第1阶段的余量信号r₁(t)：

$ r_1(t)=\mathrm{IMF}_1^{\prime}(t)+x(t) $

(3)

式中，x(t)为滑坡变形原始信号。

3) 计算第k阶段的余量信号r_k(t)：

$ r_k(t)=\mathrm{IMF}_k^{\prime}(t)+r_{k-1}(t) $

(4)

式中，IMF_k^′(t)为第i次分解的模态分量均值。

4) 重复步骤3)，直至不再满足EMD模型的剩余分量要求，即可将滑坡变形原始信号的分解表示为：

$ x(t)=\sum\limits_{i=1}^K \operatorname{IMF}_i^{\prime}(t)+R(t) $

(5)

式中，R(t)为分解后的最终剩余分量。

以降噪误差比SNR^[4]为降噪性能评价指标：

$ \mathrm{d} n_{\mathrm{SNR}}=10 \lg \left(\frac{P_s}{P_g}\right) $

(6)

式中，P_s为滑坡原始变形数据序列的功率，P_g为滤除误差后滑坡变形数据序列的功率。dn_SNR值越大，对应的分解效果越差；反之，分解效果越好。

通过CEEMDAN模型将滑坡变形数据分解为若干模态分量，考虑到滑坡变形数据的混沌特性^[5]，进一步提出对分解后的各模态分量进行相空间重构(phase space reconstruction，PSR)，以重新还原滑坡原始变形序列的非线性特征。该方法的基本原理见文献[4]，此处不再赘述。该过程的重点是确定嵌入维数和延迟时间参数，结合以往使用经验，使用C-C法确定2个参数。

极限学习机(extreme learning machine, ELM)是通过下降策略来改进预测模型结构的。本文以其构建滑坡变形预测模型，训练过程为：

$ y_j=\sum\limits_{i=1}^l \beta_i g\left(w_i x(i)+b_i\right) $

(7)

式中，y_j为滑坡变形预测值，l为隐层节点数，β_i、w_i和b_i为连接权值，g(x)为激励函数。

据使用经验，并顾及滑坡变形数据的不确定变形量，ELM模型对核函数比较敏感性，这会在一定程度上降低预测精度。引入正则化系数对其进行优化^[6]，以提升预测稳定性及泛化能力，并将优化后的预测模型命名为KELM模型。

值得指出的是，通过KELM模型对各模态分量进行预测处理虽然很大程度上保证了预测精度，但其预测结果也会存在误差。为进一步提高预测精度，再引入ARIMA模型进行误差预测，并将预测值加至前述KELM模型的预测结果中，其预测模型为：

$ R_t=\sum\limits_{m=1}^p \varphi_m r_{t-m}-\sum\limits_{j=1}^q \theta_j a_{t-j}+a_t $

(8)

式中，R_t为误差的补充预测值，φ_m、θ_j为自回归参数和滑动系数，r_t-m为预测误差值，p、q为回归阶次，a_t-j、a_t为对应白噪声值。

综上，本文预测模型的算法流程如图 1所示。

以现有变形成果和外推预测成果为基础，分别构建指标V₁和V₂，V₁为现有监测成果中最后4期变形速率均值，V₂为外推预测结果中的4期变形速率均值，通过V₁、V₂构建出累积变形预测条件下的预报等级标准，如表 1所示。

1.2 M-K分析法

为保证分析结果的准确性，采用M-K分析法对滑坡累积变形序列和速率变形序列进行趋势分析^[7-8]。结合M-K分析法的基本原理，先计算秩系数Z：

$ Z=\left\{\begin{array}{l} (S-1) / \sqrt{\operatorname{var}(S)}, S>0 \\ 0, S=0 \\ (S+1) / \sqrt{\operatorname{var}(S)}, S<0 \end{array}\right. $

(9)

式中，S为初步统计量，var(S)为初步统计量的特征参数，计算公式为[n(n+1)(2n+5)]/18，其中n为分析样本总数

通过Z值即可开展滑坡变形趋势评价，若Z值大于0，其对应评价序列具增加趋势；若Z值小于0，其对应评价序列具减小趋势。同时，Z值越大，其对应的变形趋势性越强。结合检验水平a对应临界值Z_a，可进一步划分变形趋势等级，具体标准如表 2所示。

将M-K分析法结果与前述预测结果进行对比，以佐证分析结果的准确性。

2 实例分析 2.1 工程概况

滑坡地质环境条件如下：

1) 地形地貌。滑坡区具河谷地貌，总体地形起伏较大，河谷呈V型，斜坡坡度较大，多位于30°~65°；滑坡体位于河流左岸，河流流向为110°，滑坡体主滑方向为200°，两者夹角约65°。

2) 地层岩性。滑坡区土层主要为第四系松散堆积层，岩性主要为冲洪积卵石层、崩坡积碎石层和残坡积粉质粘土层。其中，卵石层粒径约10~60 cm，稍密-中密状态，母岩成分以变质砂岩为主；碎石层主要分布于坡脚，由块状碎石为主，结构松散，粒径约60~160 cm，厚度相对较大；粉质粘土呈灰黄色、黄色，可塑，含有一定的碎石或角砾。

3) 地质构造。滑坡区地质构造较为复杂，大构造处于背斜地区，次生构造以断裂为主；同时，新构造运动也较为强烈，以地壳抬升运动和灾害发育为主。

4) 水文地质条件。滑坡区地表水主要为滑坡前缘的河流，降雨条件下坡体表面还可能形成季节性径流。地下水有两类，孔隙水主要赋存于土层孔隙中，接受冰雪融水和大气降水补给；裂隙水主要赋存于下覆基岩裂隙中，富水性受裂隙发育限制，补给与孔隙水一致。

滑坡平面具舌形(图 2)，主滑方向为200°，后缘高程约为2 760 m，前缘高程约为2 550 m，高差为210 m；滑坡体中后部较为平缓，前缘较陡，且前缘临空面较发育，临空高度约490 m；滑坡体纵向长约350 m，横向宽约250 m，平均厚度约15 m，体积约为120×10⁴ m³，属大型土质滑坡体。

滑坡(滑体、滑带、滑床)物质组成特征如下：

1) 滑坡体土层具有显著的分层特征，其上部为含碎石粉质粘土，下部为碎块石层，其中含碎石粉质粘土呈黄褐色，可塑-硬塑，中等强度，碎石含量约20%~35%，磨圆度较差，母岩以变质砂岩为主；碎块石层主要是坡体崩坡积成因，以次棱角状为主，粒径约60~160 cm，结构松散，渗透性强，利于降雨入渗。

2) 滑坡带岩芯主要为粉质粘土，灰色、黄色，厚度5~20 cm，夹杂一定量的碎石，含量一般约5%，碎石磨圆度较好。

3) 滑床主要为三叠系炭质板岩，强-中风化，黑色，岩芯主要呈块状，遇水易软化，岩层产状为205°∠40°。

滑坡体纵剖面I-I′的示意图如图 3所示。

2.2 变形预报 2.2.1 既有变形特征

布设3个监测点进行滑坡变形监测(图 2)。监测频率为2 d一次，监测时段为2020-10-28~2021-01-20，共得到42期监测成果。

1) 累积变形基础特征分析。经统计，得到3个监测点的累积变形曲线如图 4所示。

由图 4可知，3个监测点的累积变形值总体呈持续增加趋势，仅局部期数减小，其中J₁监测点的累积变形值为475.76 mm，J₂监测点的累积变形值为485.39 mm，J₃监测点的累积变形值为373.05 mm，说明滑坡体主轴线上的变形程度相对较大。

2) 变形速率基础特征分析。统计得到3个监测点的速率变化曲线，如图 5所示，可以看出，3个监测点的变形速率具较强的波动起伏特征，总体呈前期速率相对略大、后期速率相对略小的特征。

结合图 5，统计3个监测点的变形速率基础特征参数。J₁监测点速率变化范围为-12.37~45.73 mm/期，速率平均值为11.33 mm/期，速率方差值为135.08 (mm/期)²；J₂监测点速率变化范围为-1.71~43.88 mm/期，速率平均值为11.56 mm/期，速率方差值为98.67(mm/期)²；J₃监测点速率变化范围为0.24~18.01 mm/期，速率平均值为8.88 mm/期，速率方差值为24.95(mm/期)²。J₁监测点具有相对较大的速率变化范围，且变形速率的波动性也较强；J₂监测点的速率平均值最大，但仅略大于J₁监测点；J₃监测点的速率变化范围最小，速率均值也最小。

综上可知，该滑坡体在本次监测时段内变形特征明显，且变形速率具显著波动特征，说明该滑坡体潜在失稳风险较大，有对其展开研究的必要。

2.2.2 预测结果

本节重点对累积变形序列展开预测分析。首先通过CEEMDAN模型开展滑坡累积变形数据的去噪分解处理，且为掌握其优化过程的合理性，对不同过程模型均进行去噪分解结果统计(图 6，去噪模型1为EMD模型，去噪模型2为EEMD模型，去噪模型3为CEEMDAN模型)。由图 6可知，3个去噪分解模型的处理效果存在较大差异，CEEMDAN模型效果最好。

通过CEEMDAN模型将滑坡变形数据分解为若干模态的真实变形量和不确定变形量，并在其基础上开展预测分析。在预测过程中，将前37期数据作为训练集，后5期数据作为验证集，外推预测周期数为4期；同时，为充分验证预测模型构建过程的合理性，以J₁监测点为例，详述不同过程预测模型的有效性。过程预测模型分类为：预测模型1为ELM模型；预测模型2为CEEMDAN-ELM模型；预测模型3为CEEMDAN-KELM模型；预测模型4为CEEMDAN-PSR-KELM模型；预测模型5为CEEMDAN-PSR-KELM-ARIMA模型。

计算得到上述5个预测模型在J₁监测点的预测结果，如表 3所示。可以看出，不同过程预测模型的预测效果存在较大差异，预测模型1相对误差的平均值为3.23%，方差值为0.017 7%²，训练时间为231.71 ms；预测模型2相对误差的平均值为2.91%，方差值为0.015 5%²，训练时间为273.82 ms；预测模型3相对误差的平均值为2.70%，方差值为0.007 9%²，训练时间为253.91 ms；预测模型4相对误差的平均值为2.32%，方差值为0.004 4%²，训练时间为216.52 ms；预测模型5相对误差的平均值为2.02%，方差值为0.001 9%²，训练时间为194.27 ms。预测模型5，即CEEMDAN-PSR-KELM-ARIMA模型在本文滑坡变形预测中具有最好的预测效果。

类比J₁监测点的预测过程，进一步利用CEEMDAN-PSR-KELM-ARIMA模型开展其余监测点的预测研究，结果如表 4所示。可以看出，J₂监测点相对误差的平均值为2.03%，方差值为0.001 6%²，训练时间为197.27 ms；J₃监测点相对误差的平均值为2.00%，方差值为0.002 1%²，训练时间为200.27 ms。

3个监测点的预测效果相当，验证了CEEMDAN-PSR-KELM-ARIMA模型在滑坡变形预测中的普适性，且具有较优的预测效果。

由表 4的外推预测结果可知，3个监测点的后续变形仍会进一步增加。按照表 1的标准进行统计，得到3个监测点在累积变形预测条件下的预报等级结果，如表 5所示。可以看出，3个监测点的预报等级存在一定差异，J₁监测点为Ⅱ级，J₂和J₃监测点为Ⅲ级，总体判断该滑坡体目前的预报等级为Ⅲ级，属较危险状态，存在较大的潜在威胁，应增加监测频率，并尽快执行防灾预案。

2.2.3 M-K分析结果

利用M-K分析法开展滑坡累积变形序列和速率变形序列趋势判断，计算得到累积变形序列，如表 6所示。可以看出，J₁监测点的Z值为0.956，具增加趋势，显著性等级为Ⅰ级；J₂监测点的Z值为1.023，具增加趋势，显著性等级为Ⅰ级；J₃监测点的Z值为1.623，具增加趋势，显著性等级为Ⅰ级。由此可知，滑坡累积变形仍会进一步增加，与前文预测结果一致，相互验证了各自分析结果的准确性。但就显著性等级而言，各监测点的发展趋势相对一般，并不是很强。

利用M-K分析法开展速率变形序列的趋势研究，结果如表 7所示。可以看出，3个监测点的发展趋势分析结果也存在一定差异，J₁监测点的Z值为-0.476，具减小趋势，显著性等级为Ⅰ级；J₂监测点的Z值为-2.011，具减小趋势，显著性等级为Ⅱ级；J₃监测点的Z值为-1.063，具减小趋势，显著性等级为Ⅰ级。由此可知，各监测点速率序列均具减小趋势，这与累积变形的趋势判断和预测结果较为一致。

综上可知，该滑坡体累积变形仍会进一步增加，但变形速率趋于减弱，变形总体趋于稳定。

3 结语

1) 通过计算验证了CEEMDAN-PSR-KELM-ARIMA模型在滑坡变形预测中的普适性，且3个监测点的后续变形仍会进一步增加，该滑坡体目前预报等级为Ⅲ级-橙色预报，属于较危险状态。

2) 通过M-K分析法进行趋势判断，得出滑坡体累积变形具有增加趋势，变形速率序列具有减小趋势，后续变形虽整体趋于稳定，但危险性仍会进一步增加。