滑坡的发生不仅给人类的生命、财产安全带来严重威胁，还给资源、环境、生态等各方面带来巨大破坏。滑坡灾害的预测工作是滑坡灾害综合防治的有效途径之一，但由于滑坡发展演化的影响因素众多，使得滑坡运动具有开放性、复杂性和不确定性等特点^[1]，导致很难采用传统方法对其进行准确预测。以智能算法为主要研究手段开展的滑坡位移预测研究，如BP(back propagation)、RBF(radial basis function)神经网络以及SVM(support vector machine)等得到广泛应用，但实践表明，上述算法存在如下不足^[2-3]：

1) BP算法存在训练时间长、网络收敛速度慢、易陷入局部极小的问题；

2) RBF算法在网络训练过程中需以“聚类算法迭代”为基础选取基函数中心，当建模样本数量较大时，训练耗时较长，也易陷入局部极小；

3) SVM存在如何采取适当惩罚函数来防止过拟合、核参数选取困难等问题。

极限学习机(extreme learning machine，ELM)作为近年来发展起来的一种单隐层前馈神经网络学习算法，具有泛化能力强、学习速度快、网络参数设置简单等优点，逐渐被相关学者引入到滑坡变形预测中来。如周超等^[4]利用小波分析对滑坡位移序列进行分解，进而利用ELM算法对其分解量进行预测，并将其应用于三峡库区八字门滑坡位移预测；李骅锦等^[5]基于ELM与OS-ELM分别对长江三峡库区白水河滑坡的趋势项和周期项进行建模预测，取得了较好的预测效果。然而，由于ELM的输出权重直接由最小二乘估计方法得出，如果训练数据存在粗差干扰，则输出层权值的最小二乘估计结果会很差，从而导致预测结果失真。为此，本文将M估计与ELM相结合，提出一种基于M估计的Robust-ELM滑坡变形预测方法，以减少滑坡监测数据中粗差对ELM预测的干扰。

1 ELM算法

ELM算法^[6]较传统的单隐层前馈神经网络，克服了传统梯度算法易陷入局部极小的局限性，具有网络调节参数少、学习速度快及泛化能力强等优点。ELM网络由输入层、隐含层和输出层组成。设m、L、n分别为输入层、隐含层和输出层的节点数，假定存在N个不同的样本(x_i, t_i)，其中网络输入x_i=[x_i1, x_i2, …, x_im]^T∈Rⁿ，网络的期望输出t_i=[t_i1, t_i2, …, t_im]^T∈R^m，则ELM训练模型可表示为^[7-8]：

$ {y_j} = \sum\limits_{i = 1}^L {{\beta _i}g\left( {{w_i}{x_j} + {b_i}} \right)} ,j = 1,2, \cdots ,N $

(1)

式中，y_j表示ELM期望输出，β_i为隐含层第i个神经元与输出层神经元间的连接权值，w_i表示输入层神经元与隐含层第i个神经元的连接权，b_i为隐含层神经元的阈值，g(x)为激励函数。假定训练样本数量N与隐含层神经元数量L相等，则对于任意给定的w_i、b_i，ELM都可以零误差逼近训练样本，即

$ \sum\limits_{j = 1}^N {\left\| {{t_j} - {y_j}} \right\|} = 0 $

(2)

顾及式(1)，则式(2)可改写为：

$ \sum\limits_{i = 1}^L {{\beta _i}g\left( {{w_i}{x_j} + {b_i}} \right)} = {t_j},j = 1,2, \cdots ,N $

(3)

其矩阵表达式为：

$ \mathit{\boldsymbol{H\beta }} = \mathit{\boldsymbol{T}} $

(4)

式中，H为隐含层输出矩阵，T为网络输出矩阵。ELM的网络训练目标可以归结为非线性优化问题：

$ \begin{array}{l} {\rm{Minimize}}:\left\| {\mathit{\boldsymbol{H\beta }} - \mathit{\boldsymbol{y}}} \right\|\;\;\;\;\;{\rm{and}}\;\;\;\;\;\left\| \mathit{\boldsymbol{\beta }} \right\|\\ {\rm{s}}.\;{\rm{t}}.\;\;\sum\limits_{i = 1}^L {{\beta _i}g\left( {{w_i}{x_j} + {b_i}} \right) = {y_j}} ,j = 1,2, \cdots ,N \end{array} $

(5)

当g(x)无限可微时，连接权值w_i和阈值b_i在网络训练开始时可随机给定，则隐含层输出矩阵H固定不变，那么ELM的训练过程可以看作求解线性系统Hβ = T关于$ \hat \beta $的最小二乘解：

$ \mathit{\boldsymbol{\hat \beta }} = {\mathit{\boldsymbol{H}}^ + }\mathit{\boldsymbol{T}} $

(6)

式中，H⁺为矩阵H的Moore-Penrose广义逆。H⁺的计算方法有多种，其中$ {\mathit{\boldsymbol{\hat \beta }}_{{\text{LS}}}} = {({\mathit{\boldsymbol{H}}^{\text{T}}}\mathit{\boldsymbol{H}})^{ - 1}}{\mathit{\boldsymbol{H}}^{\text{T}}}\mathit{\boldsymbol{T}} $为β的最小二乘估计^[9-10]。

2 ELM输出权值辨识

考虑误差V的存在，将式(4)改写为：

$ \mathop {\mathit{\boldsymbol{H}}}\limits_{N \times L} \mathop {\mathit{\boldsymbol{\beta }}}\limits_{L \times m} {\rm{ = }}\mathop {\mathit{\boldsymbol{T}}}\limits_{N \times m} + \mathop {\mathit{\boldsymbol{V}}}\limits_{N \times 1} $

(7)

其误差方程为：

$ \mathit{\boldsymbol{V}} = \mathit{\boldsymbol{H\hat \beta }} - \mathit{\boldsymbol{T}} $

(8)

式中，$ \mathit{\boldsymbol{\hat \beta }} $为待估参数，且满足V^TPV =min的约束条件(P为观测值权阵)，则

$ \mathit{\boldsymbol{\hat \beta }} = {\left( {{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{PH}}} \right)^{ - 1}}{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{PT}} $

(9)

当权重P为单位矩阵时，式(9)的形式与$ \mathit{\boldsymbol{\hat \beta }} $_LS=(H^TH)^－1H^TT完全相同。

3 基于M估计的Robust-ELM

M估计准则为^[11]：

$ \sum\limits_{i = 1}^n {\rho \left( {{V_i}} \right)} = \min $

(10)

通常残差V为未知参数的函数，将式(10)对参数X求导，并令其等于0，以求出极值点:

$ \partial \sum\limits_{i = 1}^n {\rho \left( {{V_i}} \right)/\partial \mathit{\boldsymbol{X}}} = \sum\limits_{i = 1}^n {\rho '\left( {{V_i}} \right)\frac{{\partial \mathit{\boldsymbol{V}}}}{{\partial \mathit{\boldsymbol{X}}}}} = 0 $

(11)

顾及式(7)、式(8)，则式(11)可改写为:

$ \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{H}}^{\rm{T}}}\frac{{\rho '\left( {{V_i}} \right)}}{{{V_i}}}\left( {\mathit{\boldsymbol{H\beta }} - \mathit{\boldsymbol{T}}} \right)} = 0 $

(12)

令$ \frac{{\rho '({V_i})}}{{{V_i}}} $为权函数P(V_i)，即

$ \mathit{\boldsymbol{P}}\left( {{V_i}} \right) = \frac{{\rho '\left( {{V_i}} \right)}}{{{V_i}}} $

(13)

则式(12)变为：

$ \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}\left( {{V_i}} \right){V_i}} = 0 $

(14)

即

$ \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}\left( {{V_i}} \right)\left( {\mathit{\boldsymbol{H\beta }} - \mathit{\boldsymbol{T}}} \right)} = 0 $

(15)

其矩阵表达式为：

$ {\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}\left( {{V_i}} \right)\mathit{\boldsymbol{H\hat \beta }} - {\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}\left( {{V_i}} \right)\mathit{\boldsymbol{T}} = 0 $

(16)

式(16)与最小二乘估计中的法方程式(11)在形式上完全一致，不同之处在于，式(16)利用权函数矩阵P(V) =diag(P₁(V_i), P₂(V_i), …，P_n(V_i))代替了观测值权阵P，则在M估计下ELM网络输出权值矩阵为：

$ \mathit{\boldsymbol{\hat \beta }} = {\left( {{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}\left( \mathit{\boldsymbol{V}} \right)\mathit{\boldsymbol{H}}} \right)^{ - 1}}{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}\left( \mathit{\boldsymbol{V}} \right)\mathit{\boldsymbol{T}} $

(17)

由式(15)、式(16)可知，式(17)的求解过程类似最小二乘估计的计算过程，不同的是，权函数P (V_i)是残差的函数，计算前未知，通过赋予一定的初始值，由迭代计算方法估计参数$ \mathit{\boldsymbol{\hat \beta }} $。

在迭代计算中，根据残差值的大小，赋予各观测值不同的权重，通过反复计算，迭代完成时，含有粗差的观测值权重趋近于零或某个极小值，从而达到抗粗差的目的。利用上述估计下求解ELM的输出权值$ \mathit{\boldsymbol{\hat \beta }} $，本文称之为Robust-ELM(简记为RELM)。

RELM抗粗差性是通过等价权原理实现的，权函数的选取有多种不同的形式，比较常用的有Huber法、一次范数最小法、P范数最小法、IGG法等^[11]，本文选取一次范数最小法，即

$ \mathit{\boldsymbol{P}}\left( \mathit{\boldsymbol{V}} \right) = \frac{1}{{\left| \mathit{\boldsymbol{V}} \right| + K}} $

(18)

式中，K是一个很小的量。

4 Robust-ELM算法流程

给定任意训练样本(x_i, t_i)，其中x_i=[x_i1, …，x_in]^T∈Rⁿ，t_i=[t_i1, …，t_im]^T∈R^m，激励函数为g(x)，设网络隐含节点为L，L≤N，则RELM训练步骤如下：

1) 随机给定输入权值a_i，隐含层节点的偏差b_i；

2) 根据激励函数g(x)，计算隐含层输出矩阵H；

3) 选取最小二乘估计的$ \mathit{\boldsymbol{\hat \beta }} $ =(H^TH)^－1H^TT为迭代初值，计算初始残差V；

4) 由选取的权函数P(V)，求取每个观测量的初始权值；

5) 利用式(17)求解ELM的输出权值$ \mathit{\boldsymbol{\hat \beta }} $⁽ⁱ⁾；

6) 迭代计算，重复步骤1)~5)，直到max|$ \mathit{\boldsymbol{\hat \beta }} $⁽ⁱ⁾－$ \mathit{\boldsymbol{\hat \beta }} $⁽ⁱ⁾|≤ε或迭代次数达到预设的最大次数时，迭代结束。

5 Robust-ELM滑坡预测中的抗粗差性验证 5.1 数据源 5.1.1 链子崖滑坡

链子崖危岩体是长江流域稳定性最差的大型灾害性崩滑体，与北岸新滩滑坡隔江对峙，上距秭归县城15.5 km，下距三峡大坝26.5 km，其地理位置紧扼长江航道咽喉，一旦崩滑将对长江航运及附近城镇居民带来严重的影响。自1989年以来，我国相关部门针对链子崖的防治工作，采用相关大地测量监测设备，对T₀~T₆缝危岩区、T₇缝危岩区、T₈~T₁₂缝危岩区3个重点区域布设大量监测点进行定期监测。G_A监测点位移曲线如图 1所示。

5.1.1 古树屋滑坡

古树屋滑坡位于湖北省境内，沪-蓉-西高速宜昌至恩施段东起宜昌长江公路大桥，该高速在巴东县境内，地质环境属于典型的喀斯特地貌，加之沟壑纵横，其断面跨越2道峡谷，穿越8座高山，线路走向复杂，在巴东古树屋段施工过程中发现边坡有较大变形，为了保障高速公路施工安全，对古树屋滑坡进行监测，其主滑方向中部3#监测点的变形时序曲线见图 2所示。

5.2 Robust-ELM对单个粗差的抵御性验证 5.2.1 链子崖滑坡

利用链子崖滑坡体G_A点进行如下实验：

1) G_A点共计11期观测数据，在第6期(1983年)观测值44.93 mm中加入2 mm粗差(粗差值的大小根据监测数据趋势和经验确定)，即假定第6期观测数据为46.93 mm。

2) 利用加入粗差后的滑坡数据进行ELM预测：建模样本为1978~1986年9期数据，检验样本为1987~1988年9期数据；网络拓扑结构为：输入节点为2，输出节点为1，隐含层神经元数量为6，激励函数为Sigmoid。

3) Robust-ELM预测：以实验2)中ELM预测结果为基础，构建Robust-ELM预测模型，其中迭代权函数为一次范数法，以最大迭代次数100次为收敛条件。

将2)、3)结果与正常情况下(不含粗差时)的预测残差进行对比，如图 3所示。

5.2.2 古树屋滑坡

利用古树屋滑坡体3#点进行如下实验：

1) 3#点共计33期观测数据，在第28期观测值7.4 mm中加入2 mm粗差(粗差值的大小根据监测数据趋势和经验确定)，即假定第28期观测数据为9.4 mm。

2) 利用加入粗差后的滑坡数据进行ELM预测：建模样本为1~28期数据，检验样本为29~33期数据；网络拓扑结构为：输入节点为3，输出节点为1，隐含层神经元数量7，激励函数Sigmoid。

3) Robust-ELM预测：以实验2)中ELM预测结果为基础，构建Robust-ELM预测模型，其中迭代权函数为一次范数法，以最大迭代次数100次为收敛条件。

将2)、3)结果与正常情况下(不含粗差时)的预测残差进行对比，如图 4所示。

由图 3、图 4可知，当滑坡监测数据中不含粗差时，利用ELM进行滑坡预测所获取的残差值大体符合正态分布。由于粗差的加入，使得预测残差波动起伏较大，破坏了原有残差概率分布。利用RELM对算例1、算例2的检验样本预测结果见表 1(单位mm)、表 2(单位mm)。

从表 1及表 2可以看出，Robust-ELM具有抗粗差的性质，而传统ELM不具备该性质。

5.3 Robust-ELM对多维粗差的抵御性验证

为验证Robust-ELM对多个粗差的抵御性，利用古树屋滑坡体3#点进行如下实验：

1) 3#点共计33期观测数据，在第6期、第18期、第28期观测值中分别加入2 mm粗差。

2) 利用加入粗差后的滑坡数据进行ELM预测：建模样本为1~28期数据，检验样本为29~33期数据；网络拓扑结构为：输入节点为3，输出节点为1，隐含层神经元数量为4，激励函数为Sigmoid。

3) Robust-ELM预测：以实验2)中ELM预测结果为基础，构建Robust-ELM预测模型，其中迭代权函数为一次范数法，以最大迭代次数100次为收敛条件。

将2)、3)的预测结果进行对比，如图 5所示。

由图 5可知，当滑坡数据中含有多个粗差时，利用ELM进行预测时，所得预测结果与真实值存在很大偏差，且预测残差整体偏大；利用Robust-ELM经过100次迭代计算后，在第6期、18期、28期粗差点位上，预测残差明显小于ELM预测结果，且所得预测值与真实值较为逼近，说明Robust-ELM也具有抵抗多个粗差的性质。

6 结语

1) Robust-ELM滑坡预测模型，该方法具有良好的泛化能力，其不仅延续了ELM的快速学习特点，且对存在明显粗差的数据集，根据其每个样本残差的大小赋予不同的权重，从而降低了粗差对预测结果的影响。

2) 基于M估计的Robust-ELM算法能较好地抵御滑坡数据中的单个、多个粗差，且预测精度较高。

3) 在构建Robust-ELM算法过程中，并未讨论不同的迭代权函数及对Robust-ELM预测性能的影响，该工作是后期深入研究的重点。