基于多重分形特征和分项组合预测联合响应的滑坡预警预测研究

2022-08-30周晓岚王永强

大地测量与地球动力学 2022年9期

雷恒周晓岚王永强

1 小流域水利河南省高校工程技术研究中心，河南省开封市东京大道1号，475004 2 黄河水利职业技术学院水利工程学院，河南省开封市东京大道1号， 475004 3 长江科学院，武汉市黄浦大街23号，430010

滑坡是常见的地质灾害之一，具有频发性高、危险性大等特点[1-2]。为切实保证滑坡灾害的安全防治，开展滑坡灾害预警预测研究具有重要意义，已成为目前研究热点。

在滑坡预警研究方面，董远峰等[3]在滑坡变形特征分析基础上开展监测预警研究，周辀等[4]结合数值模拟构建滑坡预警判据。上述研究虽涉及滑坡预警，但缺乏滑坡多重分形特征研究，且邓小鹏[5]研究认为，基于滑坡变形监测成果的滑坡预警分级具有可行性。因此，仍可进一步拓展基于滑坡变形数据的多重分形特征研究及预警分级评价。

在滑坡预测研究方面，黄晓虎等[6]构建具有阶跃型变形特征的滑坡预测模型；李秋全等[7]则在滑坡变形趋势分析基础上，利用优化支持向量机进行滑坡变形预测。上述研究在滑坡变形预测方面取得一定成果，但忽略了滑坡变形数据中误差信息对预测精度的影响，且赵淑敏[8]研究认为，通过对误差信息进行分离处理，能有效提高预测精度。因此，进一步开展滑坡变形的分项组合预测具有必要性。

基于上述研究成果认为，滑坡变形预警预测研究十分必要，且仍有进一步拓展的空间。因此，本文基于滑坡变形数据，首先利用多重分形消除趋势波动分析(multifractal detrended fluctuation analysis, MF-DFA)和M-K分析进行滑坡变形的多重分形特征研究及预警分级评价；其次，在利用集成经验模态分解法对滑坡变形数据信息进行分离处理基础上，通过优化递归神经网络和混沌理论实现滑坡变形的分项组合预测；最后，将多重分形特征研究结果和变形预测分析结果进行联合响应，综合评价滑坡变形规律，以便为其防治提供一定的理论指导。

1 基本原理

在前述分析基础上，将论文分析思路进一步细分为：1)先利用MF-DFA模型开展滑坡变形数据的多重分形特征分析，在此基础上构建预警判据及指标，以实现滑坡预警分级，旨在为其防治奠定理论基础。2)先通过集成经验模态分解法(ensemble empirical mode decomposition，EEMD)实现滑坡变形数据的分离处理，并通过优化递归神经网络(recurrent neural networks，RNN)和混沌理论(chaos theory，CT)构建其分项组合预测模型，以实现滑坡变形预测研究。

1.1 多重分形特征分析模型构建

MF-DFA模型属于多重非均匀分形方法，其不仅能揭示变形序列的多重分形特征，还能有效评价变形的发展趋势，优越性较为明显[9]。MF-DFA模型已被广泛应用于岩土领域，因此利用其构建滑坡多重分形特征分析模型具有可行性。

结合MF-DFA模型基本原理，以滑坡变形数据为基础，首先对其进行累积离差求解，并对累积离差序列进行子序列划分，划分依据为：子序列长度为s，个数为Ns=N/s(N为累积离差序列总数)，由于Ns存在非整数可能，难以发挥剩余累积离差节点的作用，因此提出对累积离差序列进行逆序重新划分，则所得子序列个数为2Ns。其次，在q阶波动函数阶次条件下，求解其对应的波动函数F(q,s)为：

(1)

式中，F2(s,v)为第v个子序列的方差。

据前述可知，子序列长度s值与波动函数F(q,s)值具有一一对应关系，因此，通过改变子序列长度s可得到若干(s，F(q,s))，且两者具有如下线性关系：

lnF(q,s)=C+h(q)·lns

(2)

式中，h(q)为q阶波动函数条件下的Hurst指数，C为拟合常数。

将波动函数q值取值范围设定为-8～8之间的偶数值，若h(q)值随波动函数q值变化而变化，则说明滑坡变形序列具有多重分形特征，且通过h(2)值可判断滑坡变形趋势。具体判据为：1)当h(2)值属于[0，0.5)区间时，滑坡变形具有反向持续性，变形会趋于减小，且h(2)值越小，趋势性越强。2)当h(2)=0.5时，滑坡变形具有游离性，无法判断其发展趋势。3)当h(2)值属于(0.5，1.0)区间时，滑坡变形具有正向持续性，变形会趋于增加，且h(2)值越大，趋势性越强。

最后，由于波动函数q值与h(q)值具有对应关系，通过两者求解奇异指数a(q)：

a(q)=h(q)+qh′(q)

(3)

式中，h′(q)为h(q)导函数。

以奇异指数a(q)为基础，进一步计算多重分形谱宽度参数Δa、波形中大小波动所占比例参数Δf(a)：

Δa=amax-amin

(4)

Δf(a)=Δf(amax)-Δf(amin)

(5)

Δa参数主要用于评价滑坡变形序列的多重分形谱宽度，若Δa值越大，则多重分形强度也随之越强，即波动越剧烈；Δf(a)参数主要用于评价滑坡变形序列波形中大小波动所占比例，Δf(a)值越小，大波动波形所占比例越大。依据文献[10]研究成果，通过Δa和Δf(a)参数构建滑坡预警判据，以实现滑坡预警等级划分，具体标准设定如表1所示。

表1 滑坡预警等级划分标准

为实现表1中Δa和Δf(a)参数的趋势判断，再引入M-K分析方法。根据M-K分析方法的基本原理，将其评价参数Z值表示为：

(6)

式中，S为初步统计量，var(S)为初步统计量的方差值。

通过Z值大小即可判断对应评价对象的趋势特征：Z≥Za(Za值为对应显著水平a条件下的临界值，本文将显著水平a设置为0.05，则Z0.05=1.960)说明评价对象具有增大趋势，其越大表明趋势性越强；-Za

综上所述，通过M-K分析可判断Δa和Δf(a)参数的趋势，即可实现滑坡预警分级。

1.2 分项组合预测模型构建

滑坡变形影响因素较多，致使滑坡变形数据具有较大的波动性，且受不确定因素影响，滑坡变形数据会含有一定的误差信息：

y(t)=q(t)+w(t)

(7)

式中，y(t)为滑坡变形值，q(t)为趋势序列，w(t)为误差序列。

由于存在误差信息，在一定程度上会影响预测模型的训练过程，即降低预测精度，这也体现出分形组合预测的优势。经总结分析，将滑坡分项组合预测模型的预测流程分为3个步骤：1)先利用集成经验模态分解将滑坡变形序列分离为趋势序列和误差序列。2)利用GOA-RNN模型实现趋势序列的变形预测，以得到趋势序列的预测值。3)将趋势序列预测结果的预测误差叠加至误差序列中，组成新的残差序列，并通过混沌理论实现其预测处理。将趋势序列的预测值和残差序列的预测值进行叠加，即得到滑坡变形的最终预测值。结合分项组合预测流程，将各阶段的方法原理详述如下。

首先，需对滑坡变形数据进行趋势序列与误差序列的分离处理。考虑到经验模态分解可通过自适应正交基进行时频信号处理，对非线性不平稳信号的处理能力较强，因此，以其为基础构建滑坡变形数据的信息分离模型。需要指出的是，传统经验模态分解存在模态混叠现象，对分离效果具有一定影响，为克服该问题，集成经验模态分解随之产生，其通过增加白噪声来解决频率混叠问题。具体分析步骤为：1)将具有正态分布特征的白噪声添加至滑坡原始序列中，使之重组成具有白噪声信号的新序列。2)利用传统经验模态分解对新序列进行分离，以得到若干本征分量和一个剩余项。3)对本征分量和剩余项进行均值求解，得到EEMD最终分离结果：

(8)

式中，j为本征分量个数，IMF(t)为对应的本征分量，R(t)为剩余分量。

在前述信息分离基础上，需进一步构建信息分离效果评价指标。以往研究多是利用单一指标进行评价，如均方根误差、信噪比及平滑度指标等，由于不同指标的评价原理具有一定差异，单一评价指标已难以满足应用需求，本文提出以3个基础指标为基础，构建综合评价指标k：

k=g1+g2+g3

(9)

式中，g1、g2、g3为3个基础指标对应的归一化值。k值范围为0～3，其值越大，说明信息分离效果越优；反之，分离效果越差。

其次，在信息分离基础上，再构建趋势序列的预测模型。考虑到RNN模型可一次性激活不同神经元，具有较强的即时处理能力，同时RNN模型在拟合过程中具有较快的收敛速度和更优的训练精度，因此通过RNN模型构建趋势序列的预测模型具有可行性[11]。

值得注意的是，RNN模型虽然具有显著的优越性，但也存在一定不足，如连接权值及阈值具有较强的随机性，会对预测精度造成一定影响。因此，为保证趋势序列的预测精度，采用蝗虫优化算法(grasshopper optimization algorithm，GOA)对其进行优化处理。结合GOA算法的基本原理，将其优化过程详述为：1)对粒子位置进行初始化设置，如设定搜索空间上、下界，并设置最大迭代次数等。2)对每个粒子进行适应度计算，并保存最优粒子。3)改变粒子位置，并重新计算粒子适应度值，再对比前期保存的最优粒子，若前者更优，则对最优粒子进行替换；反之，继续迭代寻优。4)当满足最大迭代次数后，终止迭代，并输出最优粒子对应的寻优参数，即完成RNN模型的连接权值及阈值优化处理。

GOA-RNN模型虽具有较优的预测精度，但由于滑坡变形的非线性特征，使得其仍会存在一定预测误差，将其预测误差叠加至误差序列中，组成新的残差序列。

最后，由于残差序列具有较强的波动性和随机性，混沌特性明显，为实现高精度预测，提出利用混沌理论构建其预测模型。首先计算残差序列的混沌指数λmax，λmax>0时表明所得残差序列具有混沌特性，可通过混沌理论进行残差序列的预测处理；然后通过延迟时间参数τ和嵌入维数M进行残差序列的空间重构处理，并以相空间中的第i个相点ψi为预测中心，计算其与最近邻点ψl之间的距离d：

(10)

当d值最小时，通过反推可得到残差序列的预测值。

将趋势序列GOA-RNN模型的预测结果和残差序列CT模型的预测结果叠加，即为滑坡变形的最终预测结果。

2 实例分析

2.1 工程概况

白家包滑坡位于宜昌市秭归县，属长江一级支流香溪河流域，距三峡工程41.2 km。据现场调查成果，滑坡纵向长度约550 m，前宽后窄，均宽约350 m，面积约19.25万m2，平均厚度约45 m(厚度变化差异较大，中部较厚)，总体积约900万m3，属深层大型滑坡。滑坡区地形起伏较大，总体呈东低西高，其东侧相对较平缓，坡度介于0°～25°；西侧相对陡峻，坡度介于35°～50°。

自2003年三峡水库蓄水以来，白家包滑坡历年均会出现不同程度变形，且随着时间持续，局部变形具有持续发展特征。经以往资料统计[12]，近年来白家包滑坡变形特征较为显著。为实时掌握其稳定状态，对其进行变形监测，其中主滑方向上布设ZG324和ZG325监测点，两者监测成果较为完备。因此，本文以其2013～2018年监测成果作为后续分析的数据来源，监测频率为1次/月，共计得到72期变形数据，并对两个监测点的变形-时间曲线进行统计(图1)。由图可知，随着时间持续，两个监测点的变形持续增加，且具有阶梯状特征，这是由于，一方面滑坡月变形具有集中性特点；另一方面，从侧面说明滑坡变形具有显著非线性特征，对其变形进行分项组合预测具有必要性。

图1 滑坡变形-时间序列Fig.1 Time series of landslide deformation

2.2 滑坡预警分级研究

利用MF-DFA模型对滑坡变形数据进行多重分形特征分析及预警分级研究。首先，计算得到两个监测点在不同阶次q条件下的h(q)值(表2)。由表可见，随着阶次q减小，h(q)值也随之减小，说明两个监测点变形均具有多重分形特征；同时，两个监测点h(2)值分别为0.714和0.694，均大于0.5，说明滑坡变形具有正向持续性，变形趋于增加，且ZG324监测点h(2)值相对更大，说明其较ZG325监测点趋势性更强。

表2 多重分形特征分析结果

其次，利用M-K分析对Δa和Δf(a)参数进行趋势判断，以实现滑坡预警分级。结果分析如下。

1)Δa指标判据结果分析。通过计算统计，得到Δa指标判据结果(表3)。由表可见，两个监测点的Z值不仅均大于0，还均大于Z0.05值，表明两者均呈增大趋势，且ZG324较ZG325具有相对更大的Z值，说明前者趋势性更强。结合表1中判据标准可知，两个监测点在Δa指标判据条件下的预警等级为Ⅱ级。

表3 Δa指标判据结果

2)Δf(a)指标判据结果分析。类比前述，再对Δf(a)指标判据结果进行统计(表4)。由表可见，ZG324监测点Z值为-2.304，具有减小趋势；ZG325监测点Z值为-1.862，具有平稳趋势，且其趋势性小于ZG324监测点。结合表1中判据标准可知，两个监测点在Δf(a)指标判据条件下的预警等级为Ⅱ～Ⅲ级。

表4 Δf(a)指标判据结果

3)最终预警结果分析。在Δa指标判据和Δf(a)指标判据结果基础上，对两个监测点的最终预警结果进行分析(表5)。由表可见，在两个判据的预警结果中，仅ZG325监测点在Δf(a)指标判据条件下的预警等级为Ⅲ级，其余均为Ⅱ级。因此按照不利原则综合得出，两个监测点的最终预警等级均为Ⅱ级，即滑坡变形趋向不利方向发展，破坏风险一般，建议加强监测频率及巡视，并做好防灾预案。

表5 滑坡最终预警结果

2.3 滑坡变形预测分析

在滑坡预警分级研究基础上，再利用分项组合预测模型进行滑坡变形预测分析，即将分析过程分为如下两步。

2.3.1 滑坡变形数据分离处理

对两个监测点优化前后的分离效果进行统计(表6)。由表可见，经优化处理，EEMD模型在两个监测点中的综合评价指标k值均有较大提高，且EMD模型的平均综合评价指标k值为2.416，而EEMD模型为2.711，表明后者具有更明显的优势。通过前述分析可知，EEMD模型具有更好的分离效果，说明其优化处理过程具备有效性。

表6 EMD模型优化前后的分离结果

2.3.2 变形分项组合预测分析

在滑坡变形数据分离处理基础上，为充分验证分项组合预测模型的稳定性及滚动预测能力，加之考虑到监测成果的周期较长，因此，将预测过程划分为两个阶段，即前期预测和后期预测，其中，前期预测样本范围为1～36周期，后期预测样本范围包含所有样本(1～72周期)。同时，在预测过程中，将相应预测阶段中的后5个样本作为验证样本，并以ZG324监测点为例，详述不同预测阶段的预测效果。

1)前期预测结果分析

先利用GOA-RNN模型对ZG324监测点趋势项序列进行预测处理，且为验证GOA算法对RNN模型的预测效果，对优化前后的预测结果进行统计(表7)。由表可见，经GOA算法优化处理，5个验证样本的相对误差值均出现不同程度减小，充分说明GOA算法对RNN模型的参数优化具备有效性；在GOA-RNN模型的预测结果中，相对误差变化范围为2.54%～2.91%，预测精度一般，从侧面说明进行后续残差序列预测的必要性。

表7 ZG324监测点前期趋势项预测结果

由表7可知，趋势项序列的预测结果存在一定预测误差，将其叠加至误差项序列中，得到ZG324监测点在前期的残差序列，并通过计算得到λmax=0.057>0，说明残差序列具有混沌特性，可利用混沌理论实现其预测处理。

通过CT模型对残差序列进行预测处理，得到ZG324监测点在前期的最终预测结果(表8)。由表可见，经CT模型处理，ZG324监测点在前期预测结果中的相对误差范围为2.13%～2.27%，平均相对误差为2.21%，相较其趋势项预测结果，预测精度得到一定提高，表明CT模型具备有效性，也初步说明GOA-RNN-CT模型在滑坡变形预测中的适用性。

表8 ZG324监测点前期最终预测结果

类比ZG324监测点前期预测过程，对ZG325监测点进行前期预测，表9为预测结果。由表可见，在ZG325监测点的前期预测结果中，相对误差范围为1.99%～2.18%，平均相对误差为2.12%，其预测效果略优于ZG324监测点。上述分析充分说明，GOA-RNN-CT模型在滑坡变形预测中的稳定性较强。

表9 ZG325监测点前期最终预测结果

2)后期预测结果分析

在前期预测基础上，再通过后期预测来验证组合预测模型的滚动预测能力，并实现外推预测，以评价滑坡变形的发展趋势。经计算统计，两个监测点的后期预测结果见表10。由表可见，两个监测点后期预测效果相当，其中，ZG324监测点后期预测结果的相对误差范围为1.96%～2.13%，平均相对误差为2.05%；ZG325监测点后期预测结果的相对误差范围为1.85%～2.04%，平均相对误差为1.99%，均具有较高的预测精度，表明GOA-RNN-CT模型在滑坡预测中不仅具有较优的预测效果，还具有较强的滚动预测能力。同时通过外推预测可知，滑坡变形仍会进一步增加，与前述趋势判断结果一致。