大华滑坡位移预测模型研究

王 伟，邹丽芳，周倩瑶，姜宇航，陈鸿杰，徐卫亚

(1.河海大学 岩土工程科学研究所，南京 210024； 2.河海大学 岩土力学与堤坝工程教育部重点实验室，南京 210024； 3.河海大学 地球科学与工程学院，南京 211100； 4.华能澜沧江水电股份有限公司，昆明 650214)

1 研究背景

水动力型滑坡指斜坡在冰川融雪、降雨、水位变动、地表径流及地下水活动等水动力因素驱动下发生的岩土体失稳灾害[1]。我国西南地区山体众多，降雨丰富，高山峡谷地貌发育，水动力型滑坡分布较多，且建设有大型水电站，库水位落差变化很大，水库蓄水导致库区许多滑坡复活，安全隐患巨大。为了降低滑坡灾害造成的一系列损失，众多国内外学者都对滑坡预测预报展开了研究，如基于降雨阈值方法以及现场监测数据(如位移、变形、孔隙水压力等)等进行预报，均取得了良好效果[2]，其中针对滑坡位移的预测仍是该领域的重点之一。

滑坡位移曲线一般是非平稳的时间序列，直接处理相对困难，许多学者将时间序列中“分解与集成”的思想运用到滑坡位移预测中[3-6]。如徐峰等[4]、杨背背等[5]、鄢好等[6]运用移动平均法、二次移动平均法，将滑坡位移分解为反映滑坡位移长期变化趋势的趋势项位移和受外界降雨、库水位变化等诱发因素控制的周期项位移，分别用不同的模型加以预测，但移动平均法容易受现象复杂性影响；文静[7]、Zhou等[8]应用小波分析法先将滑坡位移时间序列进行分解再对位移进行预测，但小波分析在预估小波分解阶次以及确定基函数方面有一定难度；姚琦等[9]应用经验模态分解法(Empirical Mode Decomposition，EMD)分解三峡库区木鱼包滑坡位移，但会出现模态混叠现象的问题，而集合经验模态分解(Ensemble Empirical Mode Decomposition，EEMD)通过加入白噪声解决了这一问题，且不需要设定分量个数，并根据时间序列特征进行高效自适应分解[10]。

滑坡变形的演化过程受到外界影响因素的作用，位移时程曲线上也呈现出与影响因素对应的变化特征，故在做滑坡位移预测研究时，学者往往将滑坡位移与外界影响因素进行综合响应分析[11]，并将外界影响因素作为预测输入量，如彭令等[12]计算了与降雨量、库水位变化及位移增量有关的6个影响因子与滑坡累计位移的关联度，并将它们作为预测模型输入量；Cai 等[13]通过研究表明同样的遗传算法-最小二乘向量机(Genetic Algorithm-Least Square Support Vector Machine，GA-LSSVM)模型，考虑外界影响因素的多因素模型预测效果明显优于单因素模型。近年来随着计算机硬件、软件技术与人工智能的发展，使用机器学习算法建立滑坡预测模型已经成为主流研究方向，其中以反向传播神经网络(Back Propagation Neural Network，BPNN)等人工神经网络为代表被广泛应用[14-15]，但这类人工神经网络往往需要大量的样本数据，在有限的学习样本下进行精确预测存在困难，且在高维空间中容易陷入局部最优。相比而言，支持向量机(Support Vector Machine，SVM)理论基础完善，在较少的样本学习样本下仍能表现出较好的预测能力，具有运行速度快、预测精度高等优点。Suykens和Vandewalle[16]对支持向量机进行了改进，提出最小二乘支持向量机(Least Square Support Vector Machine，LSSVM)，使得模型具有更快的运行速度和更好的适应性，并在滑坡位移预测中得到广泛运用[17]。而LSSVM模型的预测精度依赖于两个内置参数的选择，快速准确地选择最佳参数组合对预测精度的提高有重要意义。粒子群算法是一种群体优化算法，具有算法简洁、收敛速度快等优点，有着良好的全局寻优能力[18]。粒子群优化算法-支持向量机模型(Particle Swarm Optimization-Support Vector Machine，PSO-SVM)应用于滑坡易发性定量评价中[19]。本文采用粒子群算法对LSSVM的参数进行寻优以避免主观因素的干扰，高效地实现最佳参数的选择。

基于以上研究背景，本文以某流域典型水动力型滑坡——大华滑坡为例，采用集合经验模态分解法(EEMD)将滑坡位移分解为多个本征模态分量，提取残余项作为趋势项位移，并将剩余分量累加作为波动项位移。对分解得到的趋势项位移用多项式方程进行拟合预测；使用粒子群算法(PSO)对最小二乘支持向量机(LSSVM)进行参数寻优，并用于波动项位移的建模预测。同时对滑坡监测资料进行分析，选取当周降雨量、两周降雨量、一周最大降雨量、库水位、库水位变化、两周库水位变化、地下水位、地下水位变化、两周地下水位变化共9个因子作为滑坡波动项位移的初始影响因子，采用灰色关联度分析对输入变量进行选择，对筛选得到的因子进行核主成分分析(Kernel based Principal Component Analysis，KPCA)以降低输入因子的维度。最后将各部分预测位移叠加，实现滑坡累计位移的预测。

2 理论与方法

2.1 集合经验模态分解

EMD可以对非线性、非平稳的原始序列数据进行线性和平稳化处理，将原始的时间序列分解为一个残余函数和一组振荡函数，然后利用希尔伯特变换(Hilber Transform，HT)获得时频谱图，得到具有物理意义的频率。残余函数反映了原始时间序列的主要趋势，振荡函数又被称为本征模函数(Intrinsic Mode Function，IMF)，其包含了原信号的不同时间尺度的局部特征信息。所有复杂的时间序列都能够利用EMD算法分解为有限个本征模函数和一个反应原始序列趋势的残余函数。Wu和Huang[20]在EMD的基础上进行了改进，提出了一个新的信号分解方法，称为集合经验模态分解(EEMD)，有效地解决了EMD的模态混叠问题。EEMD过程如下：

(1)将白噪声加入原始的信号序列。

(2)利用EMD法，将加入白噪声后的序列分解为一系列IMF函数。

(3)加入不同幅值的白噪声，重复上述步骤(1)和步骤(2)。

(4)把每次分解得到的各个IMF的均值作为最后的IMF分解结果。

2.2 核主成分分析

主成分分析(Principle Component Analysis，PCA)是一种实现大量数据降维的方法，本质上讲它是一种线性变化，只能处理线性数据的降维，遇到非线性问题时，处理效果往往不尽如人意。因此在主成分分析中引入了核函数(kernel)的概念，即核主成分分析(KPCA)。核主成分分析的基本思路是通过核函数将在低维空间中线性不可分的数据通过某种映射函数隐式映射到更高维的空间中去，而无需知道映射函数的具体形式，使其在该高维空间中线性可分，之后再使用适用于线性可分数据的相关算法进行后续处理。

设Xk∈Rd(k=1，2，…，n，d为维数)是原始空间的一组原始样本数据，利用非线性映射函数Φ(x)将样本数据从原始空间映射到高维特征空间F中。假设映射数据的均值为0，即

(1)

(2)

(3)

对于上式的特征解V，存在系数T1，T2，…，Tn，使得

(4)

定义n×n的矩阵K，Kij=Φ(xi)Φ(xj)，根据式(3)、式(4)可以得到nλT=KT，其中矩阵T=(T1，T2，…，Tn)。对于主成分的选取，只需要计算函数Φ(x)在特征向量VK上的投影，即

(5)

由于上述推导都基于一个特定前提，即假设映射数据的均值为0。在现实情况中，这种假设往往是不成立的，此时引入矩阵K′代替矩阵K，则

K′=K-InK-KIn+InKIn。

(6)

式中In为n阶单位矩阵。

2.3 PSO-LSSVM耦合模型

2.3.1 粒子群算法

粒子群算法(PSO)是一种基于种群的随机优化技术，主要用于全局寻优。在PSO中，每个优化问题的解都称之为“粒子”，所有的粒子都具有一个位置向量(粒子在解空间的位置)和速度向量(决定下次飞行的方向和速度)，并可以根据目标函数来计算当前的所在位置的适应值。在每次的迭代中，种群中的粒子除了根据自身的经验(历史位置)进行学习以外，还可以根据种群中最优粒子的“经验”来学习，从而确定下一次迭代时需要如何调整和改变飞行的方向和速度，就这样逐步迭代，最终整个种群的粒子就会逐步趋于最优解。本文用粒子群算法对LSSVM模型的两个参数进行寻优。

2.3.2 最小二乘支持向量机

支持向量机(SVM)是Vapnik等1995年提出的一种基于统计学原理与结构风险最小原则的机器学习新算法[21]。它的基本思想是经由一个非线性映射，将原输入空间数据映射至一个高维度特征空间中，这样在原空间的非线性回归问题就转化成了特征空间中的线性回归问题。支持向量机对解决小样本、非线性、高维度问题有很大优势，且能避免陷入局部最优。Suykens等[16]提出最小二乘支持向量机(LSSVM)改进了支持向量机算法，原SVM训练中的不等式约束条件在LSSVM中变为等式约束。在LSSVM模型中，设训练集为：(xi，yi)∈Rn×R(i=1，2，…，l)，其中xi表示输入向量，yi表示输出向量，l表示样本数量。在特征空间中，回归函数f(x)可以表示为

f(x)=ωTφ(x)+b。

(7)

式中：ω表示权重向量；b表示偏差量；φ(x)表示将输入向量投射至高维度特征空间的核函数。

基于算法结构风险最小化原理，函数优化问题可以表示为

(8)

式中：J(ω，ξ)为构造的函数；γ表示惩罚参数；ξi表示第i个误差变量；ξ为误差变量。

用拉格朗日法求解优化问题，其结果为

(9)

式中：φ(·)：Rn→Rnh表示映射到一个更高维度的向量空间；αi为拉格朗日乘子；L(ω，b，ξ，α)为拉格朗日函数，对ω、b、ξi、α的偏微分，需满足下列要求：

(10)

式(10)可转化为

(11)

根据Karush-Kuhn-Tucker最优条件，可以得到一个线性方程，即

(12)

式中：核函数K(x，xi)=φ(xi)Tφ(x)为满足Mercer条件的对称函数。

则最小二乘支持向量机回归估计y(x)可以构造为

(13)

本文预测方法的流程见图1。

图1 预测方法流程Fig.1 Flow chart of the prediction method

3 滑坡监测数据分析

3.1 大华滑坡工程地质概况

大华滑坡地貌形态如图2所示，地质剖面图与监测点布置分别如图3、图4所示。滑坡分布高程在1 410～1 870 m之间，后缘至前缘长度约1 000 m，顺河向宽度1 060 m 左右，周边为基岩陡壁呈“圈椅”状，形成后缘及上下游侧缘被基岩陡坡围限、前缘临空的态势，属于典型的纵横等长式滑坡。整个坡体形似扇状，前宽后窄，近NS向展布，前缘顺河宽度约1.2 km。整个坡体区内地势西(后缘)高东(前缘)低，局部呈阶梯状斜坡地貌，上陡下缓。滑坡前缘已达江边，形成坡度50°～60°的岸坡。滑坡堆积物体大约为4 840×104m3，从滑动规模上属于特大型坡。水库正常蓄水位1 477 m时，滑坡体前缘约有67 m的深度位于水面以下。

图2 大华滑坡全貌Fig.2 Panorama of Dahua landslide

滑坡物质主要由表层10～50 m左右崩积土夹碎块石及下部20～50 m左右倾倒变形的全、强风化紫红色板岩组成。滑带埋深40～50 m，最大埋深80 m。由于河谷下切、侧蚀作用前缘临空，在重力长期作用下，岩层产生倾倒折断变形破坏，并多次蠕变、错落、滑动形成多期古滑坡堆积体。从力学角度分析，该坡属于牵引式滑坡，前缘发生失稳滑动，牵引后部坡体产生拉裂缝，导致整个滑坡抗滑力进一步减小。根据钻孔揭露情况，滑坡体无统一的底滑面，但局部发育次级滑动面，滑带物多由紫红色、灰绿色砾石、碎屑及泥质物组成。

图3 大华滑坡地质剖面图Fig.3 Geological profile of Dahua landslide

图4 大华滑坡监测点布置Fig.4 Layout of monitoring points for Dahua Landslide

3.2 监测数据分析

大华滑坡上布置了19套全球导航卫星系统(Global Navigation Satellite System，GNSS)接收器用于监测滑坡表面变形，GNSS接收器安装以来，对部分测点进行了连续观测，观测周期为1次/周。在这些测点中，监测点DHLD1-1的观测数据较全，所跨周期长，且变形量最大，具有较好的代表性。故本文对大华滑坡监测点DHLD1-1从2016年9月4日至2019年8月21日的变形数据、降雨数据以及2018年2月2日至2019年8月21日间的库水位调控、地下水位数据进行分析，如图5所示。由图5可知：

(1)位移速率变化趋势与降雨量变化趋势有明显的正相关性，两个快速变形阶段均发生大幅度降雨，并且位移速率的增加一般稍滞后于大幅度降雨，滞后的时间与降雨特征、滑坡内部土层渗水能力等因素有关。

(2)大华滑坡为老滑坡，在高水位蓄水影响下复活，并受库岸再造活动的影响而处于蠕滑状态，最终前缘失稳引起坍塌，因此库水位变动对大华滑坡变形及稳定性具有较大影响。库区自然江水位为1 413 m，2018年2月第1阶段蓄水升至1 437 m，2018年4月第2阶段蓄水，水位升至库区正常蓄水位1 477 m。第1阶段蓄水期间，位移变化规律不明显，但在第2阶段蓄水期间，可以明显看到位移增量有明显的减小，接近为0，即位移速率为0，且维持一段时间。这表明，第2阶段的蓄水过程增加了滑坡的稳定性。原因是随着库水位上升，滑坡浸没体积变大，有效重量减小。

(3)地下水位与滑坡位移速率之间存在较明显的正相关关系，且地下水位变化与滑坡位移速率变化之间也存在延迟现象。如图5(c)所示，地下水位第28周出现第3次峰值1 489.48 mm，随后第29周位移速率达到最大值2.37 mm/d。

图5 观测数据时程曲线Fig.5 Time-histories of monitoring data

综上所述，库水、降雨以及地下水位的作用是导致大华滑坡发生位移运动的主要因素，在后续的滑坡因子选择中应当予以充分考虑。

4 大华滑坡位移预测

以典型水动力型滑坡大华滑坡为例，以大华滑坡表面位移监测点DHLD1-1的累计位移为预测对象，取2018年1月1日至2019年8月21日以周为单位共85组累计位移数据为原始时间序列，其中2018年1月1日至2019年5月9日70组各类监测数据作为训练样本集，2019年5月16日至2019年8月21日15组各类监测数据作为测试样本集。

4.1 位移原始序列分解

首先采用EEMD对位移时间序列进行分解。在EEMD中加入振幅为0.2的白噪声，集成总数N选常用值200，分解结果如图6所示。

图6 滑坡累计位移时间序列EEMD结果Fig.6 Accumulated displacement of the landslide with time by using EEMD

由图6可以看出，原始位移时间序列被分解为5个IMF分量和一个平稳的残余项分量r；IMF频率依次减小，残余项分量r接近一条直线。很多研究常采用分别预测各个IMF分量以及残余项最后将各个预测结果相加得到总位移的方法[22]，但是由于分量较多，导致预测过程较为繁琐，大大增加了预测的所需时间，并且分离出的频率最高的IMF1分量为随机项，有很大的随机性，导致IMF1无法被准确预测。故本文借鉴邓东梅等[23]的方法，将分离出的残余项作为滑坡趋势项位移，IMF1—IMF5相加作为滑坡波动项位移，如图7所示。

图7 大华滑坡累计位移、趋势项位移与波动项位移Fig.7 Displacement of accumulated term，trend term and fluctuation term of Dahua landslide

4.2 趋势项位移预测

经EEMD后分离出的滑坡趋势项位移逼近于一条平稳增长的直线，故采用二次多项式对滑坡趋势项位移进行拟合(图7(b))，拟合得到如式(14)所示，精度指标R2达到1，预测值和真实值的绝对误差最大仅为0.08 mm，均方根误差(RMSE)为0.04 mm。

x(t)=0.002 7t2+5.334 9t+273.238 9 。

(14)

4.3 波动项位移预测

前文分析表明，降雨量、库水位和地下水位为影响滑坡位移的三大重要因素。由于降雨量和库水位对滑坡位移的影响存在一定的滞后效应，故本文选取当周降雨量X1、两周降雨量X2、一周最大降雨量X3、库水位X4、库水位变化X5、两周库水位变化X6、地下水位X7、地下水位变化X8、两周地下水位变化X9共9个因子作为滑坡波动项位移的初始影响因子备选数据库。

分别计算9个备选因子与波动性位移的灰色关联度，将计算结果列入表1中。根据灰色关联度≥0.6的选取原则，取周降雨量X1、两周降雨量X2、一周最大降雨量X3、库水位X4、库水位变化X5、地下水位X7、地下水位变化X8共7个因子作为波动项位移初始影响因子，并将它们按顺序重新编号为Y1、Y2、Y3、Y4、Y5、Y6和Y7。这7个指标之间并不是互相独立的，例如3个降雨量指标间相互影响，降雨、库水位与地下水位之间存在因果关系，如果直接将这7个因子作为输出指标，会造成数据库冗余并影响预测结果，故应该先对初始影响因子进行降维以消除因子间可能存在的多重共线性影响。

表1 波动项位移影响因子灰色关联度计算结果

为了比较KPCA方法与传统PCA方法的优劣，对7个滑坡位移影响因子分别进行PCA和KPCA处理，其中KPCA核函数设置为高斯径向基核函数，取核函数参数σ=1，得到各指标的特征值、方差贡献率以及累计贡献率，计算结果如表2和表3所示，并对两种方法的主成分累计贡献率进行对比，如图8所示。

表2 主成分分析总方差

表3 核主成分分析总方差

图8 PCA和KPCA累计贡献率对比Fig.8 Comparison of accumulated contribution between PCA and KPCA

表4 核主成分分析的成分矩阵

从表2、表3和图8可知，KPCA获取的第1个主成分的贡献率达到了85.25%，即得到的第1个主成分就涵盖了原始数据85.25%的信息，而PCA第一个主成分贡献为45.24%，几乎只有KPCA法的一半。并且，在累计贡献率超过95%的前提下，核主成分分析只需要提取2个主成分，累计贡献率为96.52%，而主成分分析则需要提取5个主成分，与原始指标相差不大。这表明核主成分分析比主成分分析有更好的数据降维能力，更适合非线性数据的降维，所提取的主成分也能保留足够的原始数据信息。所以，本文选取KPCA方法，将得到的前两个主成分用于滑坡位移的预测。

根据表4中的KPCA成分矩阵，可以计算得到KPCA两个主成分表达式为

(12)

将根据上述两个主成分计算公式以及各个降雨量、库水位、地下水位数据计算得到的Z1、Z2作为波动项训练集的输入变量，前次位移作为训练集的输出变量，建立起LSSVM预测模型，并用PSO算法对LSSVM参数进行寻优。LSSVM模型的核函数选为径向基函数(Radial Basis Function，RBF)；在用粒子群算法对LSSVM的两个参数寻优时，初始种群规模为100，进化次数为500，速度更新公式中，系数c1、c2、ω分别取值为：c1=1.7，c2=1.5，ω=0.9；得到的波动项位移预测结果及预测精度分别如图9和表5所示。

图9 不同优化方式下LSSVM模型波动项位移预测值Fig.9 Predicted displacement of the fluctuation term in LSSVM model by using different optimized methods

表5 不同模型波动项位移预测精度

由图9与表5可知：EEMD-KPCA-PSO-LSSVM模型能较好地预测滑坡波动项位移，其预测值曲线趋势与监测值曲线基本保持一致，且相对传统常用的支持向量机SVM模型及BP神经网络模型，EEMD-KPCA-PSO-LSSVM模型预测精度最好。

4.4 累计位移预测

将趋势项位移预测值与波动项位移预测值相加即可得到滑坡累计位移的预测值。 为了探究本文提出的EEMD方法和KPCA方法提高预测位移准确性的能力， 将不采用EEMD方法、 不采用KPCA方法和既不采用EEMD也不采用KPCA方法的PSO-LSSVM模型用于比较预测能力， 各个模型的累计预测结果和精度如图10和表6所示， 可知：

(1)在以上4个预测模型中，EEMD-KPCA-PSO-LSSVM的预测精度最高，PSO-LSSVM模型预测精度最差，EEMD-KPCA-PSO-LSSVM模型3种预测精度值对比PSO-LSSVM模型分别提升了78.94%、82.44%、82.61%，证明了本文采用的EEMD及KPCA优化方法的有效性。

(2)同样对比PSO-LSSVM模型，EEMD-PSO-LSSVM模型的3种预测精度指标值分别提升了66.34%、67.96%、67.39%；KPCA-PSO-LSSVM的3种预测预测精度指标值分别提升了36.74%、33.27%、33.33%，这说明EEMD方法提升预测精度的能力要比KPCA方法强。

图10 不同优化方式下LSSVM模型累计位移预测值Fig.10 Predicted accumulated displacement in LSSVM model by using different optimized methods

表6 不同优化方式下LSSVM模型预测精度

5 结 论

本文建立了大华滑坡位移预测的EEMD-KPCA-PSO-LSSVM组合模型，给出了模型的预测流程，得出以下结论：

(1)EEMD算法在时间序列分解中具有较强的适应性，本文基于因素响应分析与EEMD算法将滑坡位移序列分解为趋势项位移与波动项位移，并分别采用多项式方程与PSO-LSSVM进行拟合预测，预测结果较好；KPCA可以依据输入因子进行有效降维，避免数据集间多重共线性的影响，提高模型预测精度。

(2)通过对大华滑坡的分析可得，库水、降雨以及地下水位是影响大华滑坡位移运动的主要因素。本文选取降雨、库水位等9个因子，结合EEMD-KPCA-PSO-LSSVM组合预测模型对滑坡位移进行预测，得到的滑坡预测值与实际监测值显示出良好的一致性，均方根误差为2.59 mm，平均绝对误差为1.72 mm，平均相对误差为0.24%。该组合预测模型在滑坡位移预测研究中具有较高的应用价值。