基于傅里葉函數擬合的滑坡變形監測分析

董耀剛，頡 麗，姚正學，楊 軍，劉興榮，王翔宇

(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#變形趨勢預測曲線