基于主控因子分析与GM-IAGA-WNN联合模型的平推式滑坡位移预测研究*——以垮梁子滑坡为例

黄 健 李 桥 巨能攀 许 强 王昌明

(①地质灾害防治与地质环境保护国家重点实验室(成都理工大学) 成都 610059)(②成都理工大学 环境与土木工程学院 成都 610059)

0 引 言

平推式滑坡是广泛分布在我国西南部三峡库区、川东和川北红层地区的一种特殊结构形式的滑坡(张倬元等， 2016)。此类滑坡多发育在近水平砂、泥岩互层的岩体中，岩层倾角一般仅3～10°、斜坡坡面一般也小于20°，如四川德阳中江县垮梁子滑坡(1949，1981)、重庆巴南麻柳嘴滑坡(1998)、四川达州宣汉县天台乡滑坡(2004)、青宁乡岩门村滑坡(2007)及四川巴中南江县黑山坡滑坡(2012)等。可见，该类滑坡普遍存在，但是针对其形成机制与主控因子的研究尚存在不少争议，对其如何进行预警预报，也成为众多国内学者重点关注的焦点之一(Xu et al.，2016)。

滑坡变形曲线是一条易受内外界因素影响且具有波动性的非平稳、不光滑曲线，笼统性对滑坡历史变形数据进行单因素或单模型分析，无法达到理想预测效果。因此，现阶段的主要研究方向是基于监测数据，提取主控因子建立主控因子-滑坡位移的非线性多因子综合预测模型(李骅锦等， 2017)，相关学者已在该方面获得诸多成果。吴益平等(2007)引入层次分析的思想，将滑坡累计位移划分为趋势项位移和具有不确定性的随机项位移，提出了灰色-神经网络联合模型，取得了较好的预测结果。张俊等(2015)运用滑坡演化规律，分析滑坡主控因子，并基于时间序列将累计位移分解，分别利用不同方法构建预测模型，提高了预测精度(徐峰等， 2011； Du et al.，2013； 李骅锦等， 2016)。但是上述成果大多是基于从外部因子去分析滑坡变形规律，未考虑对滑坡变形更为敏感的坡体内部因子，“平推式滑坡”基于此进行预测结果可靠性有待提高。因此，如何将平推式滑坡各项因子合理整合并纳入滑坡位移预测模型是研究的重难点。

鉴于此，在综合分析各模型优劣的基础上，本文以垮梁子滑坡为例，基于现场多年的监测数据，将滑坡累计位移分成趋势项位移和变化项位移两部分。首先利用GM(1， 1)灰色模型对滑坡趋势项位移进行提取； 再利用灰色关联度分析和相关性分析方法对坡体内外各项因子进行筛选，将优选的主控因子导入基于改进型自适应遗传算法优化的小波神经网络模型对滑坡变化项位移进行预测； 最后，叠加趋势项位移和变化项位移得到滑坡累计位移预测值； 通过与传统BP神经网络模型、小波神经网络模型和未优化遗传算法-小波神经网络联合模型预测值的对比分析，验证了本文所构建模型的合理性和优越性。

1 GM-IAGA-WNN联合模型

滑坡在外部诱发因素(降雨、温度和人类活动等)、坡体内部变化因素(渗压、地下水位和土体含水率等)和自身地质条件(地形地貌、地质构造和地层岩性)的共同影响下，变形位移呈现出不确定性和随机性，且不易预测。影响因素对滑坡变形的影响模式各不相同，自身地质条件决定了滑坡变形的总体趋势； 坡体内部和外部的影响因素导致了滑坡位移的不规则变化(唐天国等， 2005)。因此，可以将滑坡累计位移分解成趋势项和变化项两部分：趋势项位移是随时间变化的近似单调增长的曲线，变化项位移为一个较为复杂的非线性时间序列。

1.1 GM(1， 1)灰色模型

本文运用GM(1， 1)灰色模型对滑坡位移的趋势项进行提取。设原始滑坡累计位移时间序列为X(0)={x(0)(i)，i=1， 2，…，n}，n为数据个数。将x(0)进行一次累加(Aucumulated Generating Operator， 1-AGO)以便弱化随机序列的波动性和随机性，得到新数据序列：X(1)={x(1)(i)，i=1， 2，…，n}，其中

(1)

对x(1)(t)建立下述一阶线性微分方程：

(2)

1.2 小波神经网络模型

小波神经网络(Wavelet Neural Network， WNN)是一种基于小波变换理论和人工神经网络所构建的新型神经网络模型。其基本思想是利用非线性小波基函数取代BP神经网络中通常采用的非线性函数作为隐含层节点的传递函数(Szu et al.，1992； Zhang et al.，1992)，因此小波神经网络较BP神经网络具有更好的学习泛化能力，在滑坡变化项位移的预测中预测精度更高。本文选择3层函数型小波神经网络，模型结构如图 1所示(Zhang et al.，2011)。

图 1 小波神经网络模型结构Fig. 1 Architecture of the wavelet neural network

设系统输入数据集为X={xi，i=1， 2，…，n}，n为输入数据个数，输出数据集为Y={yi，i=1， 2，…，n}，n为输出数据个数，则第i线路的输出计算公式为：

(3)

式中，K为隐含层的神经元个数;N为输出层的神经元个数;w1k，n、w2i，k为各层网络的权值；b1k、b2i为各层网络的阈值；ψa， b为隐含层激活函数为Morlet母小波基函数)，即：

(4)

(5)

式中，ak、bk为各层网络的伸缩平移系数。

σ(t)为输出层激活函数(Sigmoid函数)，即：

(6)

1.3 IAGA-WNN模型

对于WNN的优化问题可以表述为：基于样本数据集P={xi，yi，i=1， 2，…，t}，t为数据个数，通过某种手段寻找最佳网络参数，使得误差函数Error最小且满足精度要求(张松华等，2011)：

(7)

遗传算法是通过模拟自然界进化过程的一种搜索最优解的方法，具有不受空间限制和能同时搜索多个优良点，并在附近继续优化繁殖的特点，在复杂的工程优化问题中得到了广泛引用，但也存在算法收敛速度慢、计算稳定性差和效率底下的问题。为进一步提高模型预测准确度和稳定性，避免参数的盲目试算，本文应用一种改进型自适应遗传算法(Improved Adaptive Genetic Algorithm， IAGA)，考虑利用IAGA优化WNN的参数，参数主要为：伸缩系数、平移系数和各层间的权值和阈值。IAGA-WNN神经网络运算的步骤(Khoa et al.，2004)：

(1)随机初始化。对染色体进行编码，本文采用实数编码法。3层IAGA-WNN神经网络各个体编码构成如下：

akbkw1k,nw2i,kb1kb2i

(2)确定适应度值f。计算公式(Van et al.，1992)为：

(8)

(3)算法运算。由父代群体适应度利用遗传算法的3种算法(选择、交叉和变异)得到子代群体。选择算法基于适应度比例的方法进行操作，各个体被选择概率与其适应度大小成正比。交叉算法中依据一定的交叉概率Pc随机地从被复制群体中选择两个个体进行交叉，Pc自适应调整公式(李延梅， 2012)为：

(9)

最后以变异概率Pm对每一个染色体内的各个基因(个体)进行变异算法计算，Pm自适应调整公式(李延梅， 2012)为：

(10)

式中，Pm1>Pm2>Pm3，取(0，1)之间的值。

(4)循环操作。将新一代群体插入原始种群中，并计算综合体的适应度。

(5)误差判定。计算模型误差值，若达到预设值，则终止循环，得到最佳编码个体，否则回转步骤(3)。

(6)参数转化与导入。将最佳编码个体转化成相应的网络参数，并将其赋值给作为初始值进行仿真预测。

2 垮梁子滑坡变形预测应用

2.1 滑坡工程地质概况及监测设备布置情况

根据坡表变形强烈程度差异，将滑坡分为如图 2的3个区域，坡南半侧区域为滑坡Ⅰ区，滑坡北半侧区域为滑坡Ⅱ区，孤立于北部山梁上发育的小滑坡为滑坡Ⅲ区。受地形地貌控制，滑坡区域地下水渗流场为独立水文地质单元，主要补给来源为大气降水，并且坡表的洼地积水和坡体冲积水也会对地下水进行入渗补给，沿风化裂隙、构造裂隙、层面裂隙从高处向低处运移，以泉水和渗流的形式向坡体前缘沟谷内排泄。

垮梁子滑坡于2011年初开始在滑坡Ⅰ区布设滑坡监测装置，并于2013年6月先后完成了雨量监测装置2个，地表位移监测装置2个，渗透压监测装置8个，土壤含水率监测装置9个，地下水位监测装置7个。考虑数据可靠性和完整性，选取2013年6月至2017年3月作为研究时间区间，相关监测点的监测平面和剖面布置图如图 2和图 3所示。

图 2 垮梁子滑坡监测平面布置图(单位：m)Fig. 2 Floor plan of monitoring arrangement in Kualiangzi landslide(unit： m)

图 3 垮梁子滑坡监测剖面(1-1′)布置图(单位：m)Fig. 3 Section plan of monitoring arrangement in Kualiangzi landslide(unit： m)

2.2 基于GM(1， 1)灰色模型的趋势项位移预测

通过GPS01和GPS02累计位移曲线(图 4)，可见滑坡位移随时间增大而呈上升趋势，并在局部由于内外界因素的影响呈现出不规律波动和阶跃现象。为研究其规律性，利用GM(1， 1)灰色模型从滑坡累计位移中提取出不含随机性变化的趋势项位移，以便挖掘其潜在规律。

考虑到GM(1， 1)灰色模型短期预测精度较好，而长期预测由于不可期的随机扰动效果一般，为此本文采用等维灰数递补动态预测模型，在处理数据时进行“新陈代谢”，其过程为：读入数据时设立变化周期t，到达周期数时在原始数据集中加入最新预测数据，并去掉周期内最旧的数据，形成等维的新数据集，进行下一步计算，如此循环往复，直至完成所有预测。

为避免原始离散序列的光滑度影响预测精度，通过观察两个位移监测点的曲线变化，将其分成2个部分： 2012年9月至2013年9月， 2013年10月至2017年3月(考虑监测数据总体完整程度，本文重点研究2013年6月以后部分)，并设置t=5，利用Matlab R 2016年软件编程分段计算，得到GPS01和GPS02基于GM(1， 1)灰色模型的趋势项位移与实际观测值的对比图(图 4)：

图 4 趋势项位移提取值Fig. 4 Extracted values of periodic displacementsa. 监测点GPS01；b. 监测点GPS02

2.3 基于IAGA-WNN模型的变化项位移预测

滑坡变化项位移，即为滑坡累计位移实际监测值与上节提取的趋势项位移的差值。

由于平推式滑坡的演化机理尚未得到统一的解释(范宣梅等， 2008； 郭晓光等， 2013)，本文以多种监测数据为基础，展开主控因子筛选分析：初步划分为外部主控因子、内部主控因子和坡体自身稳定性因子，尽可能全面地挖掘出各因子与滑坡位移的关系，实现平推式滑坡变形趋势预测，也在一定程度上探索其演化规律。同时为了避免所建模型参数繁杂、信息冗余，每项只选择出最优因子参与计算。

如图 5所示，在每年的雨季阶段，都有向上跃迁的趋势，因此，降雨是影响坡体变形的关键外部因子。同降雨相关的渗压、地下水位及含水率作为影响坡体变形的内部因子进行分析。滑坡变形的大小作为坡体自身稳定性因子参与综合预测模型的构建。各因子的具体量化关系主要是通过利用灰色关联度分析和相关性分析(包含Person相关系数和t检验值)进行计算筛选，获得最终影响坡体变形的关键因子序列。

图 5 变化项位移提取值Fig. 5 Extracted values of Variable displacements

灰色关联度和相关性分析应用十分广泛，其在滑坡敏感性分析中也有实现(王贝莉， 2013； 邓冬梅等， 2017)。本文给定灰色关联度计算中分辨系数取ρ=0.5，灰色关联度系数CG取值为0～1之间，且当CG越接近1则表明该因素对滑坡位移的影响越敏感，反之亦然。在计算Person相关性系数时，CP绝对值小于0.4时，可认定该因素与滑坡位移弱相关；CP绝对值介于0.4～0.6之间，可认定为中等相关；CP绝对值介于0.6～0.8时，可认定为强相关； 当CP绝对值大于0.8时，可认定为极强相关。在进行t检验时设置显著性水平为∂=0.05，即当t检验概率值Ct小于∂时，则拒绝原假设，认为该因素与滑坡位移总体是相关的，且Ct越小说明相关性越显著，并约定：当Ct小于0.01时，认定该因素与滑坡位移总体具有显著相关性； 当Ct介于0.01～0.05之间，认定该因素与滑坡位移具有一般相关性。

在进行主控因子分析时，先基于灰色关联度系数大小进行排列，再判断Person相关系数，最后基于t检验值确定相关系数是否满足显著性检验，以杜绝虚假相关现象。

2.3.1 降雨主控因子

降雨事件对滑坡的影响方式主要分为持续性降雨和突发性降雨两种形式。为了确定何种降雨形式对滑坡的影响程度更大，并避免相似影响因子的冗余，本文提出4种降雨类因子作为备选： ①单月降雨量：当月降雨总量； ②连续2月降雨量： 2个相邻月降雨总量； ③截至年内累计降雨量：从年初至当月降雨总量； ④当月最大日降雨量：当月内记录的最大日降雨量。

为分析降雨影响因子与滑坡位移之间的关系，以实际监测数据处理得到滑坡月位移、单月降雨量、连续2月降雨量、截止年内累计降雨量和当月最大日降雨量随时间变化曲线，如图 6所示。

图 6 滑坡单月位移及多型降雨量监测曲线Fig. 6 Monitored monthly displacement of landslide and multiple types of rainfall

由图 6可见，降雨主要集中在每年的雨季(6～10月)，在雨季期间内滑坡位移速率明显增大，各降雨备选因子曲线的波动与滑坡位移速率波动曲线有较好的一致性。

基于灰色关联度系数CG、Person相关性系数CP和t检验值Ct进行综合筛选，最终计算结果如表 1所示：

表 1 降雨影响因子相关性分析结果统计Table 1 Correlation analysis of rainfall influencing factors

结合灰色关联度分析和相关性分析，确定GPS01和GPS02的降雨影响因子为连续2月降雨量。

2.3.2 坡体内部主控因子

通过在坡体内部多种类型的监测装置而得到的监测曲线可以发现，其变化与滑坡变化项位移有一定关联性，但关联程度各不相同，结合数据的完整性提出备选因子： ①渗透压：坡体后缘渗透压(SY01)、坡体中部渗透压(SY03、SY04)、坡体前缘渗透压(SY05、SY08)； ②地下水位：坡体中部地下水位(DX05)、坡体前缘地下水位(DX06、DX07)； ③土体含水率：坡体中部土体含水率(HH02)、坡体前缘土体含水率(HH03、HH05)。

基于原始监测数据得到渗透压-时间变化曲线(图 7)，土体含水率-时间变化曲线(图 8)和地下水位-时间变化曲线(图 9)：

图 7 滑坡单月位移及多个渗透压监测曲线Fig. 7 Monitored monthly displacement of landslide and multiple types of osmotic pressure

图 8 滑坡单月位移及多个地下水位监测曲线Fig. 8 Monitored monthly displacement of landslide and multiple types of groundwater table

图 9 滑坡单月位移及多个土体含水率监测曲线Fig. 9 Monitored monthly displacement of landslide and multiple types of soil moisture contenta. HH02； b. HH03 & HH05

从图 7可以得出SY01监测曲线与滑坡月位移量曲线的变化规律较为一致，但具有一定的“滞后效应”，如2013年7月滑坡月位移达到最高值274.29imm和289.36imm，而渗透压值则在1个月之后2013年8月达到最高值37.69ikPa和86.23ikPa。SY03、SY04和SY05位于坡体中部，该部分渗透压一直处于平稳状态与时间变化关系不大。SY08呈现出明显的季节性变化，每年的6月份随着滑坡区域雨季的来临呈现上扬趋势，而后在9月份到达峰值，再缓慢回落。除SY01外其他监测点与月位移量曲线协同性不强。

GPS01和GPS02的变化项位移与渗透压备选因子的灰色关联度分析和相关性分析，最终计算结果如表 2所示：

表 2 渗透压备选因子相关性分析结果统计Table 2 Correlation analysis of osmotic pressure alternative factors

GPS01和GPS02的变化项位移与地下水位备选因子的灰色关联度分析和相关性分析，最终计算结果如表 3所示：

表 3 地下水位备选因子相关性分析结果统计Table 3 Correlation analysis of underwater alternative factors

由图 9a可见，HH02对外部降雨有明显的响应变化，呈现“下降-上升-下降”的总体趋势，且随时间推移，土体含水率变化范围不断变小，趋于平稳。HH02在2013年、2014年、2015年和2016年的极值差(最大含水率-最小含水率)为1.56， 1.24， 1.18和0.45。HH03和HH05的两条监测曲线，如图 9b所示，自监测装置安装后即开始走高至2014年3月保持平稳，后小幅回落，于2015年6月再次抬升，结合降雨监测数据未发现明显规律，与滑坡位移变化曲线无明确响应关系。

GPS01和GPS02的变化项位移与土体含水率备选因子的相关性分析，最终计算结果如表 4所示：

表 4 土体含水率备选因子相关性分析结果统计Table 4 Correlation analysis of alternative factors of soil moisture content

综合灰色关联度分析和相关性分析，确定GPS01和GPS02的坡体内部影响因子为SY01、DX06和HH02。

2.3.3 坡体自身稳定性因子

滑坡的变化项位移的大小不仅与坡体内外影响因子的作用有关，也要考虑滑坡自身的稳定性。滑坡越稳定，对坡体内外影响因子的响应就越弱，反之亦然。结合前述对垮梁子滑坡变形趋势分析，可以明显地发现滑坡在雨季位移速率增大，而在非雨季期间位移速率较小且平稳，周期内的总体变化趋势大致一致，但位移量值存在差异。因此，本文将滑坡的变形趋势过程以年(当年6月至次年5月)为周期进行划分，周期内滑坡位移量越小，滑坡越稳定。统计周期内单月滑坡位移量，并逐月叠加，得到监测点周期内位移累计增量，以此作为衡量坡体自身稳定性因子。经计算该因子与变化项位移的CG为0.570，CP为0.601，Ct为0.001，因而可以作为滑坡变形的主控因子之一。

2.3.4 IAGA-WNN模型预测结果

基于滑坡累计位移与主控因子的确定，进一步采用IAGA-WNN模型对滑坡变化项位移进行预测，模型预测过程如下：

(1)框架搭建。以2013年6月至2016年5月共计36组监测数据作为训练学习样本，以2016年6月至2017年3月共计10组为测试样本。模型采用3层网络设置，输入降雨影响因子、渗透压影响因子、地下水位影响因子、土体含水率影响因子、坡体自身稳定性因子共计5个节点，并将所有影响因子及其变化项位移进行归一化处理。

(2)网络参数的确定。依据国际上较为认可的确定隐含层个数的公式(杨涛， 2012)：

(14)

式中，m、n分别为输入层和输出层的个数。

取m=5，n=1，计算可得S取在5～11之间，最终通过试算确定网络结构为5-7-1，即输入层节点为5个，隐含层节点为7个，输出层节点为1个，学习率Ir1=0.01，Ir2=0.001，遗传优化算法中种群规模P=20，最大进化代数为100代，交叉率Pc=0.3，变异率Pm=0.05，误差精度为0.001。

(3)模型预测与对比。用已确定的网络参数对训练样本进行学习，得到预测结果，并将结果与传统BP神经网络模型和小波神经网络模型进行对比，结果如图 10所示。

图 10 变化项位移预测及对比Fig. 10 Perdiction and comparison of variable displacementsa. 监测点GPS01；b. 监测点GPS02

表 5 各模型预测精度及误差Table 5 Contrast of accuracy and error for prediction of each model

各模型预测精度及误差对比如表 5所示，结果表明：IAGA-WNN模型在GPS01和GPS02两个监测点的确定系数R2分别达到了0.976和0.943，均方根误差RMSE分别为4.297和6.916，最大相对误差分别MaxE分别为5.931imm和5.975imm，总体精度较高，与传统BP神经网络模型和小波神经网络相比有明显优势。且与GA-WNN模型进行效率对比，耗时相差无几，R2值分别提高了2.5%和2.1%，RMSE和MaxE也相应降低。

2.4 滑坡累计位移预测结果与误差分析

通过叠加GM(1， 1)灰色模型得到的趋势项位移与IAGA-WNN模型得到的变化项位移，最终获得联合模型的滑坡累计位移预测值，如图 11所示。模型预测值变化趋势与实际监测值较为一致，总体预测精度较好。GPS01模型预测值的确定系数R2值和均方根误差RMSE分别为0.968， 27.689，GPS02模型预测值的确定系数R2值和均方根误差RMSE分别为0.971， 26.383。

图 11 滑坡实际监测值与模型预测值的比较Fig. 11 Comparison of the observing displacements and the predicted valuesa. 监测点GPS01；b. 监测点GPS02

由图 11可见，误差主要发生在2013年9月和2013年10月，分析其原因可能是由于坡体内影响因子与滑坡位移变化协同具有一定的滞后性，由外部强降雨引起的坡体位移速率较大幅度增快，经过一段时间后才会反应在坡体内部监测点的数据变化上，因此模型的最终预测结果上产生了误差较大的点。

3 结 论

本文以四川德阳中江县垮梁子滑坡为例，从滑坡监测数据着手，优化筛选出滑坡位移关键控制因子，应用GM-IAGA-WNN联合模型对滑坡进行位移预测研究。取得以下主要结论：

(1)通过对监测数据的分析，发现了坡体各部位变化大致同步，但坡体中部变形较坡体前缘更剧烈。坡体变形主要受双月累计降雨量、SY01、DX06和HH02的控制。因此对于此类“平推式滑坡”的监测重点应放在外部降雨量、坡体后缘渗透压、中部土体含水率和前缘地下水位上。

(2)本文提出的联合模型，与传统BP神经网络模型、未添加遗传算法的小波神经网络模型和未优化遗传算法-小波神经网络联合模型预测结果对比，具有稳定性强、精度高的优势。尽管由于外部强降雨引起坡体内部影响因子变化具有一定的滞后现象，个别时间点的位移预测值未达到预期预测，但是模型有效降低了预测误差，总体效果较好，在实际运用中具有较高的推广应用价值。

(3)通过对监测数据的挖掘，揭示了滑坡变形与坡体内外影响因子的相互关系。为进一步开展平推式滑坡成因机理深入分析和滑坡预警预报提供了一定的技术支撑。