基于MACD指标的渐变型滑坡临滑预报模型研究

2022-11-23王家柱巴仁基铁永波高延超

水文地质工程地质 2022年6期

王家柱，巴仁基，葛华，铁永波，高延超

（1.中国地质调查局成都地质调查中心, 四川成都 610081；2.成都理工大学地质灾害防治与地质环境保护国家重点实验室, 四川成都 610059）

滑坡是世界上最为常见的地质灾害之一，每年造成大量的人员伤亡和巨大的经济损失[1-4]。据统计，2020年我国共发生地质灾害7 840 起，造成139 人死亡（失踪）、58 人受伤，直接经济损失50.2 亿元，其中60%以上由滑坡灾害引起。由于滑坡诱发因素及其破坏机制的复杂性，滑坡的预警预报既是热点也是难点[5-8]。近年来，虽有不少学者从加强监测设备硬件能力和提升监测平台数据解算能力方面进行研究，但尚未形成统一认识，不同地区适用性较差[9-11]。自20世纪60年代，日本学者Saito[12]利用蠕变试验成果，首次提出滑坡预报的经验公式以后，国内外学者先后提出了数十种模型[13-15]，其中日本学者Fukuzono[16]在Saito 模型的基础上，通过分析实验室获得的滑坡位移数据，得出速度倒数与时间的曲线有凹、凸和线性3 种关系，而大量经验数据表明，速度倒数与时间的拟合曲线以线性关系居多，因此速度倒数法成为一种广泛使用的临滑预报方法[17]，目前在矿山边坡[18,19]、自然边坡[20]和室内物理模拟[21]上，均取得了较好的应用。速度倒数与时间的拟合曲线多呈线性关系，主要与滑坡中新生成的剪切面和裂纹发展有关[22]，在临滑阶段，速度倒数与时间关系更加接近线性。选择合适的拟合曲线起始点，即临滑加速变形阶段开始加速（Onset of Acceleration，OOA）点，将显著提高速度倒数模型的预警精度[23-25]。

大量蠕变试验和斜坡变形监测资料表明[26-28]，在常规重力作用下，根据斜坡的变形-时间曲线特征可将滑坡分为渐变型、突发型和稳定型3 种类型。渐变型滑坡需要经过相对较长时间的孕育和演化才能逐步形成，一般发生于松散的土质斜坡或临空条件不好、具有时效变形的岩质滑坡中。渐变型滑坡从变形到启动具有明显的3 阶段蠕变特征，为监测预警预报的实现提供了有利条件，目前滑坡预警成功的案例均以渐变型滑坡居多[29-31]。

平滑异同平均线（moving average convergence divergence，MACD）是阿佩尔于1979年提出的，其利用长期和短期平均线之间差离值，反应股票整体趋势的变化情况，对股票买进和卖出时机作出研判[32]。由于该指标能较好的反映数据趋势突变的时间点，本文引进该指标，寻找滑坡变形趋势的突变点，即滑坡OOA 点。因此，本文以渐变型滑坡为研究对象，利用经济学上已非常成熟的MACD指标快速寻找滑坡OOA 点，结合速度倒数模型，提出渐变型滑坡临滑时间预测模型，并结合典型实例，对该方法准确性进行评估，期望研究成果为滑坡预警预报提供一种快速而有效的手段。

1 基于MACD 指标的临滑时间预测模型

1.1 速度倒数模型

如前所述，日本学者Fukuzono[16]通过滑坡变形监测数据，得到速度倒数与时间的关系，如式（1）所示。当α＞2 时，拟合曲线呈凸形；α＜2 时，拟合曲线呈凹形；α=2 时，拟合曲线呈直线型（图1）。

图1 速度倒数值与时间的关系曲线图Fig.1 Relationship between inverse velocity and time

式中：v—变形速率/（m·s-1）；

α—常量系数，是材料属性的函数，本文取α=2；

A—常量系数，是材料属性的函数；

tf—滑坡实际发生时间/s；

t—时间变量/s。

1.2 MACD 模型指标构建技术方法

MACD利用短期均线与长期均线之间的聚合与分离状况，对预测事物的高点及低点进行研判，其原理是利用不同速度（短期与长期）的平滑移动平均线（exponential moving average，EMA）计算两者的正负差，又称离差值（differential value，DIF），再由DIF的多日平均计算离差平均值（discrete exponential average，DEA）。经过双重平滑运算，便可较好的算出趋势的转折点位；将离差值与离差平均值画在以时间为横轴，MACD值为纵轴的坐标轴上，通过观察DIF和DEA的方向、绝对位置和相对位置关系，把它们的同向、异向和交叉作为判断趋势翻转的时机。具体步骤如下：

（1）计算短期平均线（SMA）与长期平均线（LMA）：

vt—t时刻变形速率/（m·s-1）；

n—平均数周期。

借鉴Carla 等[18]关于平均数滤波的方法，SMA中n=3，LMA中n=10。

（2）计算DIF和DEA：

式中：m—离差平均值指定期限/d。

在金融领域，m一般取为9 d，因此本模型中m=9 d。

（3）计算MACD值：

当MACD指标在-1～1 区间内震荡时，为滑坡正常变形区间；MACD＞1 且后续数值在短时间内急速增加，代表速率在短期内急剧上升，即认为滑坡进入加速变形阶段，该点便为滑坡OOA 点。

2 区布嘎滑坡预报模型实例研究

2.1 滑坡发育特征与监测设备部署

2.1.1 滑坡地质背景

区布嘎滑坡位于云南省会泽县大井镇，平面上呈近似“舌形”，横向宽约120 m，纵向长约280 m，平均厚约10 m，体积约33.6×104m3，主滑方向297°（图2），剖面上整体坡度约25°，前缘沟道切割较陡约50°，中-后缘为宽缓台地，坡度约15°，为中型土质滑坡（图2）。滑坡边界较为明显，后缘以斜坡陡缓交界处为界，左右两侧以深切冲沟为界，具有明显的双沟同源特征，滑体三面临空。滑体北侧由于前缘临空条件良好，局部复活，为滑坡强变形区，该区纵向长约130 m，宽约45 m，厚约10 m，体积5.85×104m3；滑体南侧为弱变形区，目前滑坡变形特征不明显。

图2 区布嘎滑坡平面分布图与剖面图Fig.2 Geological map and profile of the Qubuga landslide

根据现场调查，老滑坡体主要以残坡积粉质黏土夹角砾为主，基岩主要为下寒武统沧浪铺组红井哨段泥质粉砂岩，新近滑坡体（强变形区）变形迹象明显，中部发育横向延伸长度约40 m 的拉张裂缝，基本贯穿整个滑坡体。滑体两侧受滑坡变形影响，形成错落台坎，高约1 m。

2.1.2 滑坡监测设备部署

在强变形区前缘和后缘各布设1 台GNSS，弱变形区中部布设1 台GNSS，雨量计和GNSS 基站均布设于滑坡影响区以外的稳定区域（图3）。

图3 区布嘎滑坡监测设备布置图Fig.3 Location and distribution of the installed monitoring instruments

2.1.3 滑坡变形失稳过程

2020年9月3日，雨量计显示，最大小时降雨量达34 mm，单日降雨量达68 mm（图4）。2020年9月9日，滑坡体开始出现变形，至2020年9月17日此轮强降雨结束，滑坡强变形区变形特征明显，中部形成宽约50～80 cm 的贯通拉裂缝，后缘新增多级拉张裂缝，宽约20～80 cm，延伸长度5～20 m，滑坡两侧陡坎下错加剧，深度约5～8 m，如图5 所示。

图4 滑坡区单日降雨量与累计降雨量图Fig.4 Daily and cumulative rainfall of landslide area

图5 强降雨后滑坡区变形破坏特征Fig.5 Deformation characteristics of landslide after heavy rainfall

2.2 滑坡临滑时间预报

强变形区GNSS-1 和GNSS-2 位移和速率数据显示，2020年8月30日，累计位移和变化速率开始急剧上升，滑坡进入临滑阶段（图6）。累计位移-时间曲线具有明显的初始变形阶段、等速变形阶段和加速变形阶段，基本符合渐变型滑坡的3 阶段特征。弱变形区GNSS-3 累计位移平稳增加，变化速率在一定范围内波动，可见滑动变形主要发生在强变形区，后期现场复核也证实此观点。因此着重选取GNSS-1 和GNSS-2监测点进行临滑时间预报。

图6 GNSS 监测累计位移与速率时效变形曲线Fig.6 Curve of accumulative displacement and velocity-time of GNSS

2.2.1 滑坡位移速率MACD指标

以监测设备安装时间即2020年5月26日为起点，记为1，并以天为单位对位移速率进行解算，按前文所述步骤计算GNSS-1 和GNSS-2 位移速率的MACD指标，并以柱状图形式表达，见图7、图8。

图7 GNSS-1 位移速率MACD 柱状图Fig.7 MACD histogram of deformation rate of GNSS-1

图8 GNSS-2 位移速率MACD 柱状图Fig.8 MACD histogram of deformation rate of GNSS-2

2.2.2 确定滑坡OOA 点

根据MACD值以及短期平均线与长期平均线交叉情况综合分析。GNSS-1 监测数据显示，在监测开始后第98 天（即2020年8月31日），MACD指标由负变为正，且随后急剧增加，最大值达3.304；短期平均线同样在8月31日由下至上穿越长期平均线，随后与长期平均线形成较大偏离，这意味着速率在短期内急剧增加，因此判断2020年8月31日为GNSS-1 监测区域（即滑坡前缘）OOA 点，也是速度倒数模型拟合曲线的开始点。

GNSS-2 监测数据表明，在监测开始后第100 天（即2020年9月2日），MACD值由负变正，且随后急剧增加至4.486；短期平均线同样在9月2日由下至上穿越长期平均线，随后偏离长期平均线，因此9月2日为GNSS-2 监测区域（即滑坡后缘）OOA 点。

由2 处GNSS 监测数据可以看出，当短期平均线向上穿过长期平均线时，代表滑坡在短时间内速度增加较快，如果MACD值仍在-1～1 的区间，代表速度虽然有所增加，但变形趋势并未整体改变；一旦MACD值超过-1～1 的区间，代表变形趋势整体发生改变，便可将短期平均线上穿长期平均线的点，同时也是MACD值由负转正的点，作为滑坡OOA 点。

可以看出，利用GNSS-1 和GNSS-2 监测数据计算出的滑坡OOA 点存在一定的差异，其中GNSS-1 所在区域OOA 点为2020年8月31日，而GNSS-2 所在区域OOA 点为2020年9月2日。究其原因，主要是由于滑坡体中部发育的完全贯通的长大拉张裂缝，将强变形区切割成较为独立的两部分，造成两部分滑体OOA 点存在一定差异。

2.2.3 利用速度倒数模型进行滑坡临滑预报

为降低GNSS 监测数据的系统误差和偶然误差，需对采集的速率数据进行预处理。Carla 等[25]认为，指数平滑法（ESF）兼容了全期平均和移动平均所长，不舍弃过去的数据，但是随着数据的远离，仅给予逐渐减弱的影响程度，即赋予逐渐收敛为0 的权数，对于短期趋势预测有很好的效果，计算公式见式（5）：

式中： β—平滑系数，取 β=0.5。

为选取最优的数据处理方法，利用短期平均（n=3）、长期平均（n=9）以及指数平滑法，分别对GNSS-1和GNSS-2 监测数据进行线性拟合。结果显示GNSS-1 监测点（图9）中，原始数据和长期平均线数据（LMA）的拟合效果较差，其确定系数（R2）低于0.9，而利用指数平滑函数优化后的数据，其R2最大为0.981；同时就预测精度而言，LMA预测结果最差，误差率达8.23%，ESF 预测结果最优，误差率仅为0.04%。GNSS-2 监测数据拟合效果（图10）与GNSS-1 监测数据基本一致，LMA预测结果最差，误差率达9.06%，其R2最低，仅为0.507，ESF 预测结果最优，误差率仅为0.45%。综上，本文与Carla 等[25]研究结果基本一致，ESF 对短期趋势预测具有较好的效果。

图9 GNSS-1 速度倒数拟合曲线图Fig.9 Fitted curve of inverse velocity of GNSS-1

图10 GNSS-2 速度倒数拟合曲线图Fig.10 Fitted curve of inverse velocity of GNSS-2

鉴于指数平滑函数在拟合程度和预测结果上表现出的优势，本文采用指数平滑函数对GNSS-1 和GNSS-2 采集的位移速率数据进行处理后，利用速度倒数模型进行临滑时间预报。预测结果见表1 与图11。由于完成线性拟合需至少3 个数据点，因此GNSS-1的第98 与99 天无预测数据，预测数据由第100 天开始；GNSS-2 的第100 与101 天无预测数据，预测数据由第102 天开始。

图11 预测滑坡时间与真实滑坡时间Fig.11 Predicted time and real landslide time

表1 预测时间与真实滑动时间误差表Table 1 Error rates of predicted time and real landslide time of GNSS-1 and GNSS-2

预测结果显示，通过指数平滑处理后，即便在数据量较少的情况下，预测误差率也不超过2%。GNSS-1 预测数据显示，在OOA 点后第4 天，预测误差率便低于1%，滑坡发生滑动时间在第106.56～107.22 天之间，与实际发生滑动的第107 天（即2020年9月9日）接近，模型显示出较好的预测结果。同时，随着时间推移，监测数据不断更新，预测时间在真实时间上下波动且不断逼近真实滑动时间。GNSS-2 预测数据显示，在OOA 点后第6 天，预测误差率低于1%，滑坡发生滑动时间在第106.08～106.52 天之间，预测结果在前期数据量较少的情况下，误差较大，随着监测数据不断更新，预测结果不断逼近真实滑动时间。