基于改进灰色系统模型的库岸边坡变形分析

刘冠豪++杨品杰++叶险峰

摘 要：边坡变形监测与预测预报是边坡灾害防治的重要技术措施，也是目前研究的热点和难点。该文针对边坡变形监测时间间隔不定、数据样本少、预测精度要求高等问题，以某库岸边坡上的D02监测点为例，采用拉格朗日插值算法进行非等时距序列转换，结合灰色系统原理，建立非等时距灰色预测模型进行变形分析与预测。研究结果表明：（1）库岸边坡在该点的垂直位移变化经历了加速变形-减速变形-缓慢变形三个阶段，其变形速率分别为-1.599mm/月、-0.103mm/月和-0.039mm/月，总体呈暂时趋稳态势；（2）建立的非等时距灰色预测模型有效解决了因天气、施工和人为等因素影响，导致观测周期难以保持一致，出现原始数据时距不等的问题，且预测结果具有很高的精度，有一定的推广应用价值。

关键词：库岸边坡 拉格朗日插值 灰色系统模型 变形分析与预测

中图分类号：G64 文献标识码：A 文章编号：1674-098X（2015）01（a）-0051-04

随着社会经济发展和能源需求增大，我国先后开工或建成了许多大型水利水电工程，如三峡大坝、小浪底水库、南水北调工程、金沙江梯级电站等。大型水利水电工程库区的水位与水文变化，导致库岸山体稳定性降低，诱发大量不稳定边坡、崩塌、滑坡、泥石流等地质灾害，常造成山体变形破坏、河道堵塞、泥沙淤积、工程施工建设与运营使用受损，甚至引起工程失效和溃坝等灾难，已成为社会关注的热点。1963年意大利262 m高的瓦依昂拱坝因左岸山体滑坡而溃坝，造成1925人死亡。1976年美国93m高的提堂土坝溃决造成11人死亡、25000人无家可归[1]。2008年汶川地震造成震区山体破坏，诱发大量次生地质灾害，导致近2000座水库受损[2]。因此，对库岸边坡进行变形监测，并根据监测数据进行边坡变形分析与预测预报是一项重要的技术工作，对保障水利水电工程施工与运营安全具有重要意义。

库岸边坡变形破坏受众多不确定性因素影响，是一个复杂的非线性过程。其变形分析涉及变形数据处理与分析、变形物理解释和变形预测预报等方面[3]。目前，国内外学者对边坡变形监测数据处理与分析做了深入研究，取得了一定的理论研究成果，并发挥了实用效益[4]。早期基于岩土力学原理，结合工程经验，采用有限元（1967）、离散元[5]、边界元[6]和非连续变形分析[7]等方法，建立边坡变形值与相关环境参数值的函数关系式进行边坡变形分析。但因受其力学机理粗浅或假设不合理、以及参数获取等限制，致使计算结果与实际相差较大，预测精度受到很大影响。随着概率与数理统计、时间序列分析、非线性理论和人工智能等现代数学与计算理论的诞生和发展，国内外学者在此基础上建立了众多边坡变形分析模型与方法[8-11]。每种变形分析模型都有一定的适用条件和局限性，单一的模型或方法已难以准确预测出库岸边坡的变形情况，需要改进或建立已有的模型方法，使其更好地分析边坡实际变形状态与趋势[1]。库岸边坡变形的影响因素众多，往往不能用单一变量表示，同时监测时间常为不等间隔，因此灰色预测理论需在多变量非等步长研究上取得突破。该文针对边坡变形监测时间间隔不等，且影响变形因素众多、数据样本少、预测精度要求高等问题，引入拉格朗日插值算法进行非等时距数据序列转换，结合灰色系统原理，建立非等时距灰色预测模型，进行边坡变形分析与趋势预测，并采用后验差方法检验预测精度。

1 研究方法

1.1 拉格朗日插值法

1.1.1 插值算法原理

在库岸边坡变形研究中，可以用函数y=f（x）来表示其变形值在某点随时间变化的数量关系。但这种函数关系没有明显解析式，通常通过实验观测得到y=f（x）在某区间[a，b]上的一系列点的函数值yi=f（xi），即一张表（xi，yi），（i=0，1，2，…，n）。根据该函数表“插出”（构造）一个既能反映函数f（x）特性、又计算简单的函数（x）来近似f（x），并使得（xi）=f（xi），（i=0，1，2，…，n）。这种求函数近似式的方法称为插值法，其数学表达式为：

设函数y=f（x）定义在区间[a，b]上，x0，x1，…，xn是[a， b]上取定的n+1个互异节点，且在这些点处的函数值f（x0），f（x1）， …，f（xn）为已知，即yi=f（xi），（i=0，1，2，…， n）。若存在一个f（x）的近似函数（x），满足

（1）

则称（x）为f（x）的一个插值函数。f（x）为被插值函数，点xi为插值节点，式（1）为插值条件，而误差函数R（x）=f（x）-（x）称为插值余项，区间[a，b]为插值区间，插值点在插值区间内的称为内插，否则称外插。

1.1.2 拉格朗日插值多项式

通过n+1个不同的节点（xi，yi），构造一个特殊n次多项式li（x），使其在各节点上满足：

（2）

即 （3）

由条件lk（xi）=0 （i≠k）可知：x0，x1，…，xk-1，xk+1，…，xn都是n次多项式lk（x）的零点，故可设

（4）

其中Ak为待定常数。由lk（xk）=1，可求得

（5）

将公式（5）代入公式（4），则

（6）

因此，以n+1个n次基本插值多项式（公式（6））为基础的n次代数插值多项式为：

（7）

公式（7）称为n次拉格朗日插值多项式，用于解决该文灰色系统模型中的非等时距问题。

1.2 灰色模型

1.2.1 灰色原理与方法

灰色系统理论以“部分信息已知、部分信息未知”的“小样本、贫信息”不确定型系统为研究对象，用灰色数学模型将观测数据序列看做随时间变化的灰色量（过程）。通过累加生成和累减生成逐步使灰色量白化，从而建立相应于微分方程解的模型并作出预测预报。其主要模型方法与计算步骤如下。endprint

1）生成原始数列：

（8）

2）累加生成序列：

（9）

并且X（0）与X（1）满足：

（k=a，…，n） （10）

若生成的数列仍不满足条件，则可进行多次累加：

（11）

3）建立微分方程：

（12）

4）平滑差分：

（13）

5）最小二乘法求解：

（14）

公式（11）和公式（12）即为灰色预测的两个基本模型。当kn时，称为模型预测值。

1.2.2 精度评定

采用后验差方法检验灰色模型的精度，它由后验差比值C和小误差概率P共同描述。

设由灰色模型得到：

（15）

计算残差：

（16）

记原始数列及残差数列的方差分别为和，则

（17）

式中

（18）

计算后验差比值C和小误差概率P：

（19）

模型精度等级=max{所在的级别，所在的级别}，其判别式如下（表1）。

1.3 改进的灰色模型

1.3.1 等时距变换

灰色模型要求原始数据必须是等时距，但由于沉降观测持续时间长，易受到天气、施工及人为因素影响，常导致观测周期难以保持一致，从而出现原始数据时距不等。因此，在按等时距灰色理论建模型之前，需要对实测数据进行等时距变换，其具体过程如下。

取库岸边坡沉降观测断面在观测时刻的累计沉降作为原始数据序列，记为

（20）

各时段时间间隔为ti=ti+1-ti，⊿tj=tj+1-tj，i，j{1，2，…n}。当i≠j时，⊿ti≠tj，即表示各时段的时间间隔不相等。

平均时间间隔：

（21）

计算等时间间隔点的累计灰沉降值XL（1）（k），（k=1，2，…，n）：

（22）

式中，i的取值应满足，从而得到等时距数据序列：

（23）

1.3.2 非等时距建模

将拟合成一阶线性微分方程，即

（24）

式中a，u为待定常数。a为发展系数，其大小反映了序列的增长速度，u为灰色作用量。按照最小二乘法原理求解可得：

（k为时间间隔点） （25）

取时所对应的时间t=0，则k时间间隔点所对应的时间t=（k-1）。故t时刻所对应的预测沉降值为：

（26）

联合公式25和公式26，得到非等时距预测模型：

（27）

由于非等时距灰色模型是经过拉格朗日插值与等时距变换而成，故其精度评定方法与等时距灰色模型相同。

2 实例分析

2.1 形变分析

以湖南省张家界市慈利县境内的某水库坝址右岸山体一边坡为研究对象，该边坡存在变形失稳风险，直接影响水库大坝及其下游5km处的城镇安全。自2000年以来，在该边坡体上布设了沉降监测点，并采用GPS、水准测量和变形监测方法进行各监测点的形变监测，以分析边坡的变形规律与趋势。现以该边坡上的D02号监测点的垂直位移变化值为例进行边坡变形分析与预测（表2）。

利用公式（2）～（7）和公式（20）～ （23），进行等时距变化处理（表3）。

为了便于分析计算结果，则将D02监测点垂直位移变形的等时距变换值以时间为横坐标、变形值为纵坐标，绘制成变形过程线图（图1），以直观反映出观测点垂直位移变形的规律、趋势和幅度。

从图1和表3可知：D02监测点垂直位移在前5个月内沉降显著，呈加速变形趋势，变形速率为-1.599 mm/月；在第5～49月期间，D02监测点沉降速度减少，呈减速变形趋势，变形速率为-0.103 mm/月；在第49～62月期间，D02监测点沉降速度减少，呈缓慢变形趋势，变形速率为-0.039 mm/月。总体而言，边坡体上的D02监测点在监测后的5年内经历了加速变形、减速变形和缓慢变形等阶段，并趋于稳定。

2.2 模型预测

由于D02监测时间间隔各不相等，而等时距变换后的时间为非整月或半月的月份，即非规则月份，不便于模型预测和应用分析。现以0.5个月为单位，经拉格朗日插值，计算各规则监测时间间隔，则可以建立插值后的时间序列：{t}={0， 0.5， 1， 4， 5， 7.5， 9， 11.5， 14.5， 17， 18.5， 19.5， 22.5， 28， 32.5， 35.5，37.5， 40.5， 44， 47.5， 50.5，53，56，60.5， 62}。

由于该时间序列观测值仍不等时距，则采用灰色非等时距建模方法，进行累加生成和累减生成变换处理，且保证处理用于灰色非等时距预测模型的原始序列值满足非负性和递增性要求。由表3可知，经拉格朗日插值后的观测数据仍不能满足非负性要求，因此，需要在建模前进行如下处理。

注：平均时间间隔=62/24=2.58（月）

①移轴。

寻找数列中绝对值最大的负数（-13.3），然后取该负数的绝对值，并将该绝对值设为min，再对数列中每个数值加上min，得到数列{x（0）（t）} ={11.00， 6.56， 5.03， 5.27， 4.26， 4.31， 4.08， 3.47， 2.71， 2.40， 2.40， 2.38， 2.27， 2.06， 1.83， 1.69， 1.50， 1.79， 1.16， 0.51， 0.65， 1.72， 0.40， 0.40， 0.00}。

由于x（0）（2）< x（0）（1），表明移轴后数列仍然不满足建模所需的递增条件，不能建立一阶灰色非等距模型，因此，需要对移轴后的数列进行累加生成。

②累加生成。

一次累加生成，结果为：{x（1）（t）} ={11，17.56， 22.59， 27.86， 32.12， 36.42， 40.50， 43.97， 46.68， 49.08， 51.48， 53.86， 56.12， 58.18， 60.01， 61.70， 63.20， 64.99， 66.15， 66.66， 67.30， 69.02， 69.42， 69.82， 69.82}。该数列x（1）（23）=x（1）（24），仍不满足建模所需的递增条件，则需对其进行二次累加生成。

二次累加生成，结果为：{x（2）（t）} ={11， 28.56， 51.15， 79.01， 111.12， 147.55， 188.05， 232.02， 278.70， 327.77， 379.25， 433.11， 489.23， 547.41， 607.42， 669.13， 732.33， 797.32， 863.46， 930.12， 997.42， 1066.44， 1135.87， 1205.69， 1275.51}。该数列已满足非负、递增的建模要求。

③二阶建模。

利用二次累加生成的数列{x（2）（t）}，进行二阶灰色非等距建模，可得：

最小二乘原理求解，可得：a=-0.0206，u=-8.3530，则对应的非等距灰色预测模型为：

用公式（28）计算D02监测点各非等时距灰色预测值（表4、图2）。

注：表中“实际值”为监测点实际监测时间经等时距插值处理后得到的垂直位移形变值；“预测值”为按监测点时间序列进行非等时距灰色模型计算的垂直位移形变值；“-”表示没有实际观测值；残差值=实际值-观测值；相对误差=|差值/实际值|。

由表4和图2知：D02监测点的实际值和预测值比较吻合，残差值围绕y=0上下波动，初步预测模型预测良好。

2.3 精度检验

由公式（16）～（18），计算可得：

=-0.064，=-10.66，=5.14，=0.24

计算后验差比值：，即C<0.35。

计算小误差概率：1，即0.95

模型预测精度等级=max{P所在的级别，C所在的级别}=1级（好）。因此，非等距灰色模型在边坡变形预测中具有很高精度，可以满足实际工作需要。例如，该模型预测D02监测点在2007年7月9日的变形值为-15.4mm，而该点在2007年8月16日采集的实际数据计算出的变形值为-14.6mm，残差值为0.8mm、相对误差为5.48%，满足边坡变形分析的精度要求。

3 结论与讨论

1）D02监测点近5年的观测结果表明，库岸边坡在该点的垂直位移变化经历了加速变形-减速变形-缓慢变形三个阶段，在监测至第5个月、第5～49个月和第49～62个月的变形速率分别为-1.599mm/月、-0.103mm/月和-0.039mm/月，故该边坡总体呈暂时趋稳态势。

2）采用拉格朗日插值算法进行非等时距序列转换，结合灰色系统原理，建立非等时距灰色预测模型进行边坡变形分析与预测，有效解决了因天气、施工和人为等影响因素而导致观测周期难以保持一致出现原始数据时距不等的问题，且预测结果具有很高的精度，能够满足实际需要。

3）随着预测时间的延长，非等时距灰色模型的预测精度会有所降低，需要在实践应用中对模型反复进行试验分析，厘定出满足预测精度要求的数据序列的有效时间长度，以提高模型的精度与实用性。

参考文献

[1] 高玉云.基于小波理论的变形监测数据处理模型的研究[D].西南石油大学，2011.

[2] 孔亚平，韩用顺，朱颖彦.震区公路走廊次生地质灾害监测与评估地图集[M]. 人民交通出版社，2013.

[3] 王孟穹.模糊分析方法及时序分析模型在变形监测中应用[D].西南交通大学，2009.

[4] 王忠学，梁东伟.基于小波神经网络在滑坡变形监测预报中的应用[J].测绘与空间地理信息，2012，35（11）：190-192.

[5] Cundall P. A. A computer model for simulating progressive large scale movements in blocky rock systems[C]//ISRM Symp.， Nancy， France， Proc.2.1971：129-136.

[6] Crouch S.L. Solution of plane elasticity problems by the displacement discontinuity method. I. Infinite body solution[J]. International Journal for Numerical Methods in Engineering， 1976， 10（2）：301-343.

[7] Shi G.H. & Goodman R.E. Two dimensional discontinuous deformation analysis[J]. International Journal for Numerical and Analysis Methods in Geomechanics， 1985，9（6）：541-556.

[8] 石双忠，岳东杰，梅红.时间序列分析在变形监测数据处理中的应用[J].工程勘察，2004（3）：59-61.

[9] 孙树海，曹兰柱，张立新.露天矿边坡稳定性的模糊综合评判[J].辽宁工程技术大学学报，2007，26（2）：177-179.

[10] Prakash N. & Manikandan S.A.， et al. Prediction of biosorption efficiency for the removal of copper using artificial neural networks[J].Hazard Mater，2008：1268-1275.

[11] 付建军，邱山鸣，赵海斌.基于灰色关联度的边坡稳定影响因素分析[J].长江科学院院报，2011，28（1）：53-58.