基于傅里叶函数拟合的滑坡变形监测分析

董耀刚，颉 丽，姚正学，杨 军，刘兴荣，王翔宇

(1.甘肃省科学院 地质自然灾害防治研究所，甘肃 兰州 730000；2.中国石油天然气股份有限公司 西南管道贵阳输油气分公司，贵州 贵阳 550000)

滑坡深部位移监测是滑坡变形监测的重要内容,对准确确定滑动面的位置、评价滑坡活动现状及未来发展趋势具有重要意义[1]。目前主要是根据位移-时间曲线确定滑坡的稳定性，根据累积位移-深度曲线和钻孔柱状图确定滑动带的位置[2]。具体应用中，前期主要是根据滑坡深部位移变形曲线形态来定性分析滑动面的位置和滑坡变形发展过程，但由于监测结果是时间序列的随机变量[3]，一些学者尝试用非线性的方法，如Apriori算法[4]、Kalman滤波数据融合[5]、小波去噪后回归[6]、信息维[7]、多重分形[8]、Verhulst-ARMA组合预测模型[9]、半参数模式的 GM(1,1)[10]等，进行监测数据挖掘，从而了解滑坡形变及其原因。这些方法的应用取得了一定效果，但也存在诸多局限，如预测精度不高、适应性不强、部分算法较为复杂等。研究计算简便、精度较高、运算快捷的模型来分析滑坡监测数据是研究者一直追求的目标。本研究以2015年G316线稍子坡段(K2556+345～K2563+000)10#滑坡深部位移监测数据为例，采用MATLAB软件中的curve fitting tool工具对监测数据进行处理，具有计算简单方便、结果通俗易懂且精度较高的特点，可为更加准确地分析滑坡的发展趋势提供基础数据。

1 滑坡特征

稍子坡段10#滑坡位于G316线K2560+526～K2560+777 km处，属于老滑坡。滑坡长137 m、宽222 m、平均厚度8.5 m、面积25 400 m2、体积约29.4万m3，滑体由扰动破碎的棕红色含砾泥岩及少量砂岩组成，主滑方向为NE43°。滑坡后缘坡陡，上方山体地形坡度较缓，滑体前缘滑坡舌推覆于河床之上，逼迫河道，前后缘相对高差180 m。滑坡左右两侧的变形破坏有明显差异：滑体左侧的中后部地势相对较高，前缘堆积在泥盆系褶皱地层上，并已被河流侵蚀呈陡坡；滑体右侧在主滑之后至少又发生过3次小范围的滑动，现已完全滑落，后部地势较低、坡面较缓，地下水埋深仅1.0 m左右，有湿地和3处泉水出露，有几处小范围浅层坡体因含水量高而明显变形[11]，同时滑坡舌受河流侵蚀经常发生坍塌。国道316线从滑坡中部通过，受滑坡左右两侧差异变形影响，通过滑体右侧的公路损坏严重。2015年初相关单位邀请专家对该路段滑坡进行考察会诊时，考虑到该滑坡规模大、投资大和治理相对困难等，专家一致建议暂不治理，先进行深部变形监测，根据监测结论再作决定。

2 滑坡深部位移监测方案

滑坡深部位移监测是通过仪器测量地下岩土体相对于稳定地层的位移量，据此确定移动岩土体的滑移面和变形变化的规律[12]。其监测原理示意和监测点照片见图1、2。

滑坡深部监测孔的位置一般选在滑坡体主滑方向和具有代表性的主滑轴线上[12]。根据该滑坡的发育特征，在公路上方，沿其主滑断面自上而下布设ZK10-9#、ZK10-8#、ZK10-5#三个监测孔，以确定滑动面深度。受人为活动影响，ZK10-8#监测孔孔口严重堵塞，因此在与ZK10-8#监测孔斜距2.5 m处补设BZ8#监测孔，以代替ZK10-8#监测孔。此外，为便于继续观测该滑坡变形特征，在公路外侧增设1个监测孔BZ1#。各监测孔基本参数如表1所示。监测仪器采用成都世纪安岩仪器有限公司生产的AY-1B型固定式测斜仪，可直接监测和读取位移数据，并有微电脑时控开关实时控制将数据定时存入电脑。

图1 滑坡深部监测原理示意

图2 滑坡深部监测孔照片

表1 滑坡深部监测孔基本参数

参数各监测孔基本参数ZK10-9#BZ8#ZK10-5#BZ1#孔口标高/m1 711.061 700.901 690.231 672.23孔深/m25.823.525.023.0测深/m25.523.024.522.5总测深/m95.5

3 监测数据成果

滑坡深部位移监测可以提供大量的观测数据，根据这些数据可以整理出各种曲线，从不同角度了解滑坡的性状和发展趋势。2015年稍子坡段10#滑坡各监测孔监测数据见表2。

表2 2015年各监测孔累计位移量

注：其中“—”表示因施工覆盖等原因导致变形不明显，未测得数据。

4 监测数据的处理及分析

4.1 累计位移-深度曲线

利用累计位移-深度曲线可以进行滑动面判断，滑动面位置一般在累计位移-深度曲线相邻一对正负曲率最大点之间，若为多滑动面或滑动带则将出现几对正负曲率极大值点。埋设在滑坡后部的监测孔ZK10-9#无明显的位移，说明滑坡后缘处于相对稳定状态；BZ8#监测孔在6.5～8.0 m处出现位移突变，说明在这个深度范围内可能存在1个滑动面(带)；ZK10-5#监测孔在7.0～11.5、16.5～17.5、21.0～23.5 m位置均出现位移突变，其中7.0～11.5 m处仅在1月28日出现位移，可认为是仪器误差等偶然事件导致，而在16.5～17.5、21.0～23.5 m处存在滑动面的可能性较大。综合3个监测孔的监测结果，可以在滑坡纵剖面上确定出滑动面位置，如图3所示。

4.2 位移-时间曲线

利用位移-时间曲线可以确定滑坡的稳定性，但想通过简单的位移变化规律来判断滑坡状态是困难的，这就需要采用一定的数学分析方法来预测滑坡深部位移的变形趋势，其实质就是一组数据的曲线拟合问题。曲线拟合的方法很多，仅在MATLAB的curve fitting tool中就提供了12大类近40种拟合模型。早在19世纪早期，法国数学家傅里叶就提出“任何函数都能用无穷多个正弦和余弦函数的和表示”，说明傅里叶函数具有广泛的适应性，特别适合对周期性函数进行描述。考虑到滑坡位移数据具有周期性的特点，本研究选择傅里叶函数进行拟合，为简便起见，选择一阶傅里叶函数进行拟合，其基本公式为

y(t)=a0+a1coswt+b1sinwt

(1)

式中：t为时间，d；y(t)为位移，mm；a0、a1、b1分别为待定系数；w为基频，hz。

为检验曲线拟合效果，一方面利用均方根误差、误差平方和,以及拟合度等统计量来检验，同时考虑到本次监测数据有限，故仅留最后一组数据(ZK10-9#监测孔)用作预测效果检验，其余数据均参与建模，另一方面选择三次多项式拟合和高斯拟合两种模型作为效果对比分析模型。三次多项式拟合模型和高斯拟合模型分别为

y(t)=p1t3+p2t2+p3t+p4

(2)

(3)

上二式中：p1、p2、p3、p4、a1、a2、b1、b2、c1、c2均为待定系数；其余参数意义同上。

采用MATLAB软件中的curve fitting tool拟合工具箱，对上述拟合函数进行最小二乘法求解，3种拟合方法的拟合结果见表3。由表3可知，各拟合函数的拟合度均大于0.9，拟合效果均较好，其中：BZ8#监测孔高斯拟合效果最好，拟合度最高，误差平方和、均方根误差最小，但预留数据组预测误差率最小的是傅里叶拟合；ZK10-5#监测孔拟合效果最佳的是三次多项式，拟合度最高，均方根误差最小，但就误差平方和、预留数据预测误差率来说，效果最佳的是高斯拟合；BZ1#监测孔误差平方和及均方根误差最小、拟合度最高、预留数据预测误差率最小的都是傅里叶拟合。总体来说傅里叶拟合效果较好，因此本研究采用傅里叶拟合函数对滑坡深部位移变形进行预测，预测结果见图4—6。由图4—6知，该滑坡深部位移仍处于变形扩展阶段，急需采取有效措施以遏制滑坡变形的进一步扩展。

5 结果与讨论

以2015年G316线稍子坡段10#滑坡深部位移监测数据为例，利用MATLAB软件中的curve fitting tool曲线拟合工具箱中的傅里叶函数对3组实测数据进行拟合分析，得出如下几点认识：

(1)滑坡深部位移监测是研究滑坡深部位移特征行之有效的手段，利用累计位移-深度曲线上发生突变的位置，结合滑坡纵剖面图，可以较准确地确定出滑动面的位置。

(2)利用曲线拟合的方法，如三次多项式、高斯函数、傅里叶函数等，可以对滑坡位移-时间曲线进行深层次处理，进而对滑坡变形发展趋势进行预测，方法简便且精度较高。预测结果显示，该滑坡目前仍处于变形扩展阶段，急需采取有效措施来遏制变形的发展。

(3)滑坡位移-时间曲线拟合效果可以通过拟合度、误差平方和、均方根误差等统计量来检测，拟合结果显示：3种拟合函数均具有较佳的拟合效果，拟合度均大于0.9，但综合来看以傅里叶函数拟合效果最佳，拟合度最高、误差平方和、均方根误差最小。

(4)鉴于傅里叶函数拟合的普适性及对周期性函数描述效果最佳的特点，可以利用傅里叶函数拟合的方法建立滑坡预测预报模型，但有关滑坡阶段划分、临滑判据等需要进一步研究。

表3 3种拟合模型拟合结果分析

注：*表示实测值。

图4 BZ8#变形趋势预测曲线 图5 ZK10-5#变形趋势预测曲线 图6 BZ1#变形趋势预测曲线