以速率为参量的GM(1,1)滑坡时间预报模型研究

，谢金

(厦门大学 a.建筑与土木工程学院, 福建 厦门 361005； b.深圳研究院， 广东 深圳 518057)

1 研究背景

滑坡是世界上比较活跃的自然灾害之一。一旦达到一定规模的滑体沿着软弱结构面下滑时，就有可能对其影响范围内的建筑物(或构筑物)、人员或财产造成严重的损害或损失。因此，为了科学地防治滑坡，减少伤亡和损失，滑坡的预测预报显得更为重要。

自20世纪60年代以来，国内外学者先后提出了30多种滑坡预测预报的方法和模型。这些方法和模型大致可分为5类：①确定性预报模型,如斋藤模型[1]、福囿斜坡时间预报法[2]等；②统计预报模型,其代表性的模型有灰色GM(1,1)模型[3]、Verhulst模型[4-5]、等；③非线性预报模型,其代表性的模型有突变理论模型[6]、协同预报模型[7]、BP神经网络模型[8]等；④宏观预报模型,该模型是根据滑坡演化过程中表现出的各种征兆、迹象，如裂缝宽度，滑带土的孔隙率等建立各类滑坡特征相适应的预报模型，模型因滑坡个体差异而不同。⑤空间预报模型,如逻辑回归模型[9-11]、随机森林模型[12]。

从摩擦学的角度看，滑坡演化的本质是构成滑坡的颗粒或者颗粒团的内摩擦行为。当用电子显微镜观察滑带土的特征时，可以发现滑带土中有许多酥裂状结构，切砾微断裂，以及放射状或不规则状裂纹；同时也以观察到碎砾定向排列和颗粒压扁拉长等变形现象，宏观上表现为位移监测曲线呈现波动上升的趋势，这正是滑坡做黏滑运动和蠕滑运动的标志[13]。目前已有的预报模型大都没有考虑滑坡的运动速率对摩擦力的影响。与此同时，现有的预报模型大多数以滑坡的位移为参量建立预报模型。

基于上述分析，有必要考虑滑坡运动速率对摩擦力的影响。该文以黏滑运动和灰色系统为理论基础，建立以速率为参量的GM(1,1)预报模型，并结合黄茨滑坡、新滩滑坡进行实例分析，取得较好的预报效果。

2 滑坡孕育特征和变形运动特性

2.1 滑坡孕育特征

图1 不同类型的滑坡位移随时间的变化Fig.1 Displacement of different types of landslides against time

许强[14]根据斜坡的变形-时间曲线特征将滑坡分为：稳定型滑坡、渐进型滑坡、突发型滑坡(图1)。该文用残余抗剪强度[15]σr、长期抗剪强度σl、峰值抗剪强度σp、流变下限强度σcl以及滑坡的下滑应力强度σ，这5个强度指标来区分滑坡的类型。当滑坡的下滑应力强度σ小于滑坡的流变下限强度σcl时，滑坡仅发生弹性变形，不发生宏观的位移；当σcl<σ<σl时，滑坡存在流变变形，有宏观的位移，但经过一段时间后，流变停止。这主要是因为在没有其他干扰(地震、突降暴雨等)的情况下，滑坡的残余强度σr足以抵抗滑坡的下滑应力强度σ，一般把这种类型的滑坡称为稳定型滑坡(图1)；当σl<σ<σp时，滑坡存在流变变形，经过一段时间的发展，滑坡的残余强度σr不能抵抗滑坡的下滑力σ，变形逐渐增大，渐渐进入加速破坏阶段，直至整体失稳破坏，这种具有滑动时效的滑坡称为渐变型滑坡(图1)；当σp<σ时，滑坡不发生流变变形，滑动面上的各点的强度均达到或超过滑坡的破坏强度，滑坡会突然失稳，发生突发型滑坡(图1)。由上述分析可见，不同的受力特征，滑坡表现出的变形特征也不同。因此，从分析滑坡的受力状态出发，分析滑坡的变形以及运动特性，或许可以找到预报滑坡的方法。

图2 岩土应力-应变关系的2种基本模式Fig.2 Two basic modes of stress-strain relationship of rock and soil

大量研究表明，滑坡不是瞬间发生的，通常具有一定的变形时效。滑坡主要发在坡度为30°～50°的斜坡上[16]。突变型一般具有较深的滑动面，往往发展为液化流动，初始蠕动阶段、等速蠕变阶段经历的时间很短，发生时往往没有预兆；渐进型滑坡规模较大，演化时间较长，具有明显的预兆，比如地表裂缝、隆起、下陷等。已有的研究表明，大多数的自然斜坡和人工边坡，随着时空的演进，多发生渐进型滑坡。渐进型滑坡的滑带土一般夹杂着具有剪胀性的超固结黏土或密实的砂土，对于剪胀性的土而言，随着下滑力和抗滑力的变化，土的密实性和黏聚力逐渐下降，表现出应变软化的现象(图2)。随着变形的发展，摩擦系数逐渐降低，岩层间的错动和弯曲越来越明显。从时间上看，渐进式滑坡一般会经历3个阶段：初始变形阶段、等速变形阶段、加速变形阶段。从空间上也可以将滑坡的变形分为3个阶段[14]：第1阶段，当斜坡的下滑应力强度σ超过滑坡的流变下限σcl时，岩层向坡下弯曲，后缘产生拉力，坡体表层产生蠕滑。后缘的拉力使得坡体的潜在滑移面产生裂缝，并不断延伸发展，逐渐形成主裂缝；与此同时，拉力的竖直分量使得坡体下错，造成面向坡下的台阶；第2阶段，随着蠕变的不断发展，中部剪应力集中的部位被扰动，使得斜坡的下半部分出现隆起或者鼓丘，裂缝继续发展。第3阶段，随着裂缝的发展，潜在的滑移面逐渐贯通，斜坡破坏，形成滑坡。

2.2 滑坡变形运动特性

无论是稳定型滑坡、渐进型滑坡还是突发型滑坡，滑坡的运动方式可以归纳为2种：一种是稳定的滑动，即滑动过程中没有应力降，如滑坡的等速滑动阶段，这种滑动形式被称为蠕滑；另一种是不稳定的滑动，滑动过程中伴随有周期性的应力降。前者滑动速度较慢，后者滑动速度较快。滑坡的孕育和发生过程实际上是这两种运动交替发生的过程，显微镜下更为常见的情况是黏滑和蠕滑标志往往共存一处[17]。

在滑动体发生蠕动之前，潜在滑动面附近的颗粒或者颗粒团在地质作用下随机分布，具有高度的对称性，处于无序状态。当斜坡的下滑应力强度σ超过滑坡的流变下限σcl时，滑坡体开始蠕动，滑带土内的颗粒开始发生明显的有序运动。通过显微镜观察可以发现，滑带土内的软弱结构部分有从内向外扩展的微裂缝。随着微裂缝的不断扩展，平行于滑动面的颗粒或颗粒团在摩擦力或者剪应力的作用下不断被挤压或者拖拽。少量的硬度较低的颗粒或者颗粒团被压碎、拉长或者形态逐渐向扁平化方向发展[18]，而硬度较大的颗粒或者颗粒团则产生滑移；垂直于滑动面的硬度比较低的颗粒或者颗粒团在正压力的作用下长轴方向开始向滑动面倾斜，而硬度较高的颗粒或者颗粒团开始被剪断；随着蠕变变形的继续发展，平行于滑动面的大部分颗粒或者颗粒团在剪应力作用下沿着滑动面发生定向运动[18]。

大量硬度比较低的颗粒或者颗粒团被压碎或者呈现扁平状，而硬度比较高的颗粒或者颗粒团逐渐填充到孔隙率比较高的区域。垂直于滑动面的大部分硬度比较低的颗粒或者颗粒团在正压力的作用下长轴方向不向滑动面倾斜的程度不断加大，而硬度较高的颗粒或者颗粒团由于被剪断成三角楔，在摩擦力的作用下发生转动[18]。通过显微镜观察可以发现，滑带土内的软弱结构部分有放射状或者不规则的裂缝形成。随着剪切位移的不断增大，裂缝不断发展，颗粒或者颗粒团被压碎、啃断的现象越来越严重，潜在的滑移面逐渐贯通斜坡破坏，形成滑坡。

通过以上的分析发现，滑坡的孕育、发展和发生过程实际上是蠕滑和黏滑运动交替发生的过程。

3 滑坡演化的预测预报模型

3.1 传统灰色预测预报模型及以速率为参量的GM(1,1)模型

3.1.1 以位移为参量的GM(1,1)预报模型

1986年，我国学者陈明东等[3]首次以位移切线角为参量建立灰色系统模型，对新滩滑坡进行了预报验证，并取得成功。由于位移切线角获取较困难，此小节将以位移为参量进行建模预报。

灰色系统模型的灰微分方程为

(1)

式中a1，b1为系数。对于一组等时距位移观测序列

x(0)=[x(0)(1),x(0)(2),…,x(0)(n-1),x(0)(n)],

x(0)一次累加后得到

用最小二乘法求解系数a1，b1，计算表达式为

(2)

其中：

(3)

Z(1)(tk)=0.5(x(1)(tk-1)+x(1)(tk)),

k=2,3,…,n;

(4)

Y=[x(0)(t2),x(0)(t3),…,x(0)(tn)]T。

(5)

求解灰微分方程得到位移的时间响应公式为

(6)

式(6)两边对时间求导，得到速度随时间的变化率为

3.1.2 以位移为参量的Verhulst预报模型

Verhulst是1987年德国物理学家发现的一种生物生长模型，由于滑坡的位移特征与生物的生长规律类似，很多学者用Verhulst模型对滑坡进行预报研究[4-5]。

Verhulst模型的灰微分方程为

(8)

式中a2，b2为系数。

Verhulst模型系数的求解方法和过程同GM(1,1)模型，唯一不同的是需要将式(3)中的矩阵B替换成如下形式，即

模型位移随时间的演化方程为

3.1.3 以速率为参量的的GM(1,1)预报模型

绝大多数机械运动中存在摩擦振动现象(黏滑运动)，滑坡运动也不例外。为了深入认识这种运动，学者们对这种现象进行大量的研究。我国学者李春波[19]进行了一系列试验，并提出用下述方程来描述摩擦力与滑动速度之间的依赖关系，即

(11)

(12)

取式 (12)前2项代入式(1)得

(13)

F=uNcosα-kv。

(14)

式中：u为滑动摩擦系数；N为垂直于坡面的支持力。

根据摩尔-库伦准则，可以将抗滑力计算式表达为

F抗滑=uNcosα-kv+c。

(15)

式中c为黏聚力。

滑坡受自重力及内力作用，且内力平行于坡面(图3)。

图3 滑坡的简化滑块模型Fig.3 Simplified sliding model of landslide

则由牛顿第二定律可知：

(17)

式(17)即为描述边坡演化的非线性动力学模型。

该模型是以速率为参量的灰色GM(1,1)预测模型。模型的求解方法同以位移为参量的灰色GM(1,1)预测模型相比，只需将位移参量替换成速率参量。

以速率为参量的GM(1,1)预测模型的速率计算式为

式 (18)对时间积分，可以得到以速率为参量的灰色GM(1,1)预测模型的位移时间公式为

(19)

3.2 模型预测预报判据和时间

许强等[20]认为改进的切线角>85°，斜坡变形进入加速变形阶段。此方法的切线角85°判据是在对一个滑坡位移曲线进行分析的基础上得出的，缺乏一般性。王家鼎[21]通过对大量滑坡实例的统计分析，提出将位移切线角89°～89.5°作为滑坡的预警判据。此后，李聪等[22]建立了滑坡数据库，对不同类型的滑坡的临界速率进行了统计，发现滑坡进入急剧变形阶段位移临界速率一般≤50 mm/d。根据切线角与速率之间的转换关系v=tanθ，计算得到滑坡剧烈滑动时的位移切线角θ一般≤88.855°。由统计规律可知，滑坡实例数据样本量越大，越能反映滑坡的变形规律。故到目前为止，根据切线角变量来确定预警判据的研究中，李聪的临界速率判据最为科学和严密；贺小黑等[5]认为Verhulst模型以速率最大值为预报判据缺乏合理的依据，以加速度最大值作为预报判据更为合理。综合上述分析，对于GM(1,1)模型，本文取临界速率vcr=50 mm/d为滑坡发生的判据；对于Verhulst模型，本文以加速度最大值为判据。将vcr分别代入式(7)、式(18)，反解出滑坡发生时间。

以位移为参量的GM(1,1)模型预报时间的公式为

以速度为参量的GM(1,1)模型预报时间的公式为

以位移为参量的Verhulst模型预报时间的公式[6]为

(22)

式(20)、式(21)、式(22)中的Δt为位移监测的等时间间隔天数。

4 实例说明

为了验证以速率为参量的GM(1,1)模型的预测预报效果，分别对黄茨滑坡、新滩滑坡进行理论计算和分析。值得注意的是，实际监测过程中，滑坡位移监测数据往往会出现波动，速率会出现负值。为了减弱观测误差对数据分析的影响，可先用Vondrak平滑法[23]对原始数据进行预处理，处理后，如果速率没有负值，则可直接进行建模计算。若速率依然有负值，可选用滑坡进入加速变形阶段后的数据进行计算和分析。 为了方便进行对比分析，将以位移为参量的GM(1,1)模型编号为模型①，以位移为参量的Verhulst模型编号为模型②，以速率为参量的GM(1，1)模型编号为模型③。

4.1 实例1——黄茨滑坡

黄茨滑坡位于甘肃省永靖县,该滑坡起始监测时刻为1994年8月1日，监测时间间隔为15 d，剧滑时间为1995年1月30日凌晨2点30分，监测资料[5]如表1。

表1 黄茨滑坡B2监测点的监测数据Table 1 Monitored displacements at monitoring point B2

建模前，先对原始数据进行平滑处理，平滑方法采用Vondrak平滑法[23]，该法对等时距观测数据和不等时间距离都能进行平滑处理。平滑效果如图4。

图4 黄茨滑坡原位移和平滑位移与测次曲线Fig.4 Original and fitted cumulative displacement vs.time curves of Huangci landslide

对于以位移为参量的GM(1,1)和Verhulst模型，用于建模的原始序列为

x(0)=[1.05，3.26，6.07，8.2，10.7，14.99，18.05，21.65，27.4，33.14，38.56，49.43，63.54]。

设时间序列为t=[1,2,3,4,5,6,7，8，9，10，11，12，13]，对x(0)进行一次累加后有：x(1)=[1，4.31，10.38，18.58，29.28，44.27，62.32，83.97，111.37，144.51，183.07，232.5，296.04]。

对于以速率为参量的GM(1,1)预报模型，由于滑坡监测的物理量是位移，故需要将位移转化为速率，才能用速率模型参量进行建模预报。位移转化为速率的计算公式如下,即

计算得到建模的速率序列为

v(0)=[2.21，2.81，2.13，2.5，4.29，3.06，3.6，5.75，5.42，10.87，14.11]。

以1994年8月16日为速率模型建模的初始时间，v(0)对应的时间序列为

vt=[1,2,3,4,5,6,7,8,9,10,11,12]。

对v(0)进行一次累加后有

v(1)=[2.3，5.1，7.3，9.5，14.3，17.1，20.6，26.5，32.6,38,47.7,63.3],用最小二乘法求解，得到模型参数a1=-0.217 61，b1=5.155；a2=-0.223 98，b2=0.424 55；a3=0.308 06，b3=0.000 29。将求得的系数以及v(1)(t1)，vcr分别代入式(20)、式(21)、式(22)得模型的预报时间值tr1=108.779 26，tr2=115.715 41，tr3=271.671 06。需要说明的是，位移监测单位是cm,vcr取5 cm/d；由于以速率为参量的GM(1,1)模型的计算时间从1994年8月16日开始，故代入式(20)时需要将式(20)的t替换成vt计算。把预报时间值转换成日期，得到滑坡预报结果见表2。表中，Δ是预报绝对误差(d)，是预报时间与滑坡实际发生时间的差值，即Δ=tcr-ts；δ为相对误差(%)，相对误差是绝对误差Δ换算成天数后除以滑坡实际发生时间ts与时间序列t1之间的间隔天数再乘以100%。当滑坡预报时间与实际发生时间相同时，Δ=0，δ=0；当滑坡预报时间早于实际发生时间时，Δ>0，δ>0；当滑坡预报时间晚于实际发生时间时，Δ<0，δ<0。下文中的计算规则与此处相同。将预报结果列于表2中。

表2 黄茨滑坡模型预报结果Table 2 Forecasted occurrence time of HuangciLandslide by different models

图5 黄茨滑坡模型位移预测值-时间曲线Fig.5 Comparison of predicted displacement vs.time curve among three models

从表2中可以看出，3个预测预报模型的预报误差均在50%以内，GM(1,1)模型的预报时间比Verhulst模型的预报时间早，能够提前预报滑坡的发生。以速率为参量的GM(1,1)模型的预报时间更加接近滑坡实际发生时间，精度较高。

将求得的模型参数分别代入式(6)、式(10)、式(19)可以求得模型位移预测值。图5是3个预测模型的位移时间曲线与原位移时间曲线和平滑位移曲线的对比图。从图5中可以看出，以速率为参量的GM(1,1)预测模型的位移预测效果明显优于Verhulst模型，基本上能够反映这段时间内位移随时间的变化趋势。

4.2 实例2——新滩滑坡

新滩滑坡位于湖北省秭归县境内的新滩镇长江北岸,其下游27 km 处是三峡大坝. 滑坡发生于1985年6月12日。将新滩滑坡B3监测点的位移数据[7]列于表3中。为了分析以速率为参量的GM(1,1)模型的预测预报效果,首先对监测点的位移数据进行平滑处理。平滑效果见图6。

图6 新滩滑坡原位移和平滑位移随时间的变化曲线Fig.6 Original and fitted cumulative displacement vs.time curves of Xintan landslide

表3 新滩滑坡B3监测点位移Table 3 Monitored displacements of point B3 of Xintan landslide

用平滑后的位移全部数据建立GM(1,1)模型。从图中可以看出，第36个月以后，位移速率不断增大，故认为从第37个月开始，滑坡进入加速变形阶段。

对于Verhulst模型，选取第36个月以后的平滑数据进行预测预报。分别用最小二乘法求解得模型系数a1=-0.048 46,b1=110.605 22，a2=-0.046 45,b2=4.413 26，a3=0.066 64,b3=1.204 5×10-7。将求得的系数以及v(1)(t1)，vcr分别代入式(20)、式(21)、式(22)得模型的预报时间值tr1=47.423 99，tr2=54.002 07，tr3=79.789 39。将预报结果列于表4中。

表4 新滩滑坡模型预报结果Table 4 Forecasted occurrence time of XintanLandslide by different models

注:amax表示滑坡运动的加速度达到最大值

从表4中可知，3个预测预报模型的误差均在50%以内，GM(1,1)模型的预报时间比Verhulst模型的预报时间早，能够提前预报滑坡的发生。以速率为参量的GM(1,1)模型的预报时间能做到提前预警。

将求得的模型参数分别代入式(6)、式(10)、式(19)可以求得模型位移预测值。图7是3个预测模型的位移时间曲线与原位移时间曲线和平滑位移曲线的对比图。

图7 新滩滑坡模型位移预测值-时间曲线对比Fig.7 Comparison of predicted displacement vs. timecurve among three models

从图7中可以看出，以速率为参量的GM(1,1)预测模型的位移预测效果明显优于Verhulst模型，基本上能够反映位移随时间的变化趋势。

5 结 论

该文以摩擦学中的黏滑运动以及灰色系统为理论基础，建立了以速率为参量的GM(1,1)预报模型，并结合黄茨滑坡、新滩头滑坡实例进行了预测预报分析，主要得到以下结论：

(1)基于监测位移建立模型对进行滑坡预测预报时，所建立的模型应该考虑滑坡的变形速率对摩擦力的影响。

(2)该文建立的以速率为参量的GM(1，1)预报模型的预报效果比传统的以位移为参量的GM(1，1)预报模型和Verhulst预报模型更好。

(3)以速率为参量的GM(1,1)滑坡预测预报模型比以位移为参量的Verhulst预测预报模型更能够反映位移随时间的变化规律。