基于多演化階段插值方法的滑坡預測研究

段 功 豪，馬 迅，付 杰

(1.武漢工程大學 計算機科學與工程學院，湖北 武漢 430205； 2.武漢工程大學 智能機器人湖北省重點實驗室，湖北 武漢 430205； 3.中國地質調查局 水文地質環境地質調查中心，河北 保定 071051)

0 引 言

滑坡是一種在山區或丘陵地區十分常見的地質災害。據自然資源部統計，2021年上半年中國共發生地質災害1 150起，造成17人死亡，3人失蹤，直接經濟損失達24 084.9萬元，其中滑坡高達579起，因此對滑坡精準的預測可以極大程度地減少損失。

影響滑坡發生的因素有很多，其中位移是滑坡預測的主要因素之一。從20世紀80年代開始，國內外普遍使用數學模型來擬合監測點的位移-時間曲線進行分析，從而預測滑坡位移的動態趨勢變化。其中影響預測結果準確性的除了模型方法，還有對滑坡監測數據的預處理方法[1-3]。階躍型滑坡具有明顯的演化區間，在實時收集監測點數據時，可能發生儀器故障或受外界影響，產生部分異常的數據，因此數據插值方法的好壞直接影響數據收集的準確性[4]。相關預處理所用的插值算法有很多，Lagrange插值法、均差的牛頓插值法、線性插值法、埃爾米特插值法和三次樣條插值法等[5]。傳統的直接插值方式是將整體的數據帶入插值算法從而算出相應插值，不僅無法保證整體曲線的光滑性，而且會產生龍格現象，極大降低了插值的精確度。所以有時需要將整條曲線進行全局分段插值，縮小數據量可部分提高準確性，但由于通過同一指標分段，降低了分段數據與不同演變階段之間的關聯，尤其對于階躍型滑坡，在拐點附近會出現較大的誤差。

因此，在對曲線分段時使用什么算法、對哪些點進行插值、應在曲線上選取哪些點作為分段點，成為學者們廣泛討論的問題[6-7]。基于上述研究背景，本文選擇目的性強、靈活度高的Lagrange插值法；為保證最后預測的準確性，插值點的選取不僅需解決數據的缺值異常值，有時還要確保數據的等距性，根據需求設置插值點；在插值分段選取上，根據不同的演化階段，以拐點作為分段點，將曲線分為相應的區間，對每個區間進行局部分段插值，既保證了插值的整體關聯性，又能穩定插值與當前運動情況的準確性；最后對白家包滑坡的累積位移-時間曲線進行基于實際運動學變化情況的插值，并對比驗證了該方法的準確性。

1 研究區概況

白家包滑坡屬于典型的階躍型滑坡，位于宜昌市秭歸縣歸州鎮向家店村，處于香溪右岸[8]。通過圖1無人機拍攝可見，周圍地貌部位為低山區凹形斜坡，前緣直抵香溪河，滑坡剪出口位于高程125～135 m處，滑坡后緣以基巖為界，高程265 m，滑坡左側以山脊下部基巖為界，右側以山梁為界，前緣寬500 m，后緣寬300 m，均寬400 m，縱長約550 m，滑坡面積22萬m2。滑坡平面形態呈短舌狀[9-11]。深層滑體前緣厚20～30 m，中部厚47 m，后緣厚10～40 m，滑體平均厚度45 m，滑體體積990萬m3。淺層滑體前緣厚10～20 m，中部厚35 m，后緣厚10～40 m，滑體平均厚度30 m，滑體體積660萬m3[12]。此地氣候類型為亞熱帶季風氣候，且有明顯的峽谷氣候特征，降水豐沛，四季分明，年平均降雨量為1 000～1 400 mm，降雨多集中于夏季，是暴雨高發地區之一。三峽峽谷深邃，風力較強，風向固定，為東西向的順峽方向。勘探資料顯示，此地滑坡體的土壤物質成分由粉質黏土混合碎塊石土組成，土質松散且雜亂。如圖2滑坡體剖面圖所示，滑床物質主要由侏羅系下巴東組的長石石英砂巖及泥巖組成，巖層傾向為260°～285°，傾角為30°～40°[13]。滑坡前緣大部分浸沒在河水中，當水庫水位上升至175.00 m時，部分滑體將完全淹沒于水中。山體植被覆蓋率較好，主要由亞熱帶闊葉林和部分山地灌草叢組成，但在山體上部有部分裸露的土壤巖石，在下部有數十戶居民，以及幾百畝開墾過的農田，其中還有一條秭興公路從中穿過，因此若此地發生滑坡，將造成巨大的損失。

圖1 白家包滑坡全貌Fig.1 Full picture of the Baijiabao landslide

圖2 白家包滑坡體剖面圖Fig.2 Section of Baijiabao landslide mass

為了保障村民的生命財產安全，白家包滑坡在2006年10月開始采用GPS對地表位移進行監測，建設有ZG323、ZG324、ZG325和ZG326共4個監測點，其中ZG324和ZG325位于滑坡主滑坡面，ZG323和ZG326位于滑坡中下部靠近秭興公路附近[14]。這4個監測點的位移方向與坡向一致，各個監測點每月不定期監測1次，夏季6～9月每月監測2次[15]。各儀器均有較高的精度，可精確到小數點后10位。

2 研究方法

2.1 插值算法

本次研究采用的拉格朗日(Lagrange)插值算法結構整齊、實現方便，具體多項式如下：

(1)

其中每個lj(x)為Lagrange基本多項式(或稱插值基函數)，其表達式為

(2)

對于分析監測點的位移-時間的數據，每個時間和對應的監測值可以作為一個坐標點，x為橫軸的時間，y為縱軸的監測值，所有的監測數據都可以通過(x，y)表示。將n個監測點(x0，y0)，(x1，y1)，…，(xn，yn)帶入n項拉格朗日多項式中，通過輸入插入時間x，可算得此刻的插值大小。

在計算插值誤差時，通常計算多項式余項來作為插值的誤差，由公式(1) 可知，y=f(x)的Lagrange插值是Ln(x)，其誤差公式為

(3)

通過公式可以發現，Lagrange插值法根據階數的不同，可以演化出不同的插值算法：當多項式是一階函數時，是線性插值算法；當多項式取二階時，是拋物插值算法；取三階時，相當于三次樣條插值。因此，根據實驗的需要，靈活地設置階數，便于進行理論分析[16]。隨著節點數的增多，階數的上升帶來了高次計算量的增大，插值點的誤差被無限放大，從而導致了龍格現象，當曲線斜率越大時，龍格現象越明顯。為避免龍格現象的發生，在分析數據較多的滑坡位移數據時，往往使用分段插值法，如分段線性插值法、三次樣條插值法等Lagrange低階插值法，但插值的整體關聯性減弱了，因此應根據不同的插值應用背景，科學設置相應的目標函數與動態實現過程。

2.2 指數平滑算法

指數平滑預測法是指通過一個平滑系數，將本期實際值與本期的預測值按比例相加，最終得到下一期數據的預測值的一種加權平均法，其依據時間序列的穩定性和有規律性的態勢，進行合理的延續來達到對近期數據的預測。根據平滑次數不同，分為多種類型，在對滑坡易發性研究時，對斜率明顯的曲線更適用于二次指數平滑法。二次指數平滑法必須與一次指數平滑法配合，即在一次指數平滑法的基礎上再做一遍指數平滑，從而建立預測的數學模型，最后確定預測值。二次指數平滑的計算公式如下：

(4)

(5)

xt+T=AT+BTT

(6)

其中，

(7)

(8)

式中：T為未來預測的期數；AT與BT都是預測模型中的參數[17]。本次滑坡預測使用二次指數平滑模型時，由于數據量較大，靠前數據對模型影響較小，可將實際值代替預測值。

3 數據分析及預測分析

3.1 數據選取

在2006～2012年之間，根據白家包滑坡設置的4個位移監測點ZG323、ZG324、ZG325和ZG326獲取到一共90期數據作為插值方法研究的對象。將獲取的累積位移和對應的時間作為橫縱坐標繪制累積位移-時間曲線圖，如圖3所示。

圖3 白家包滑坡累積位移-時間統計Fig.3 Cumulative displacement-time curve of Baijiabao landslide

觀察圖3曲線，可發現明顯的3大特點：① 具有同步性，4個監測點的累積位移的變化幅度均隨時間同步改變，因此可只選取其中1條曲線作為研究對象；② 具有周期性，累積位移在每年的5～9月份會加速增長，其他時間相對較為平穩；③ 具有明顯的增長趨勢，從整體上看，曲線呈現出階梯狀的增長趨勢，有明顯的拐點，并且隨著時間推移，累積位移的變化越來越明顯，符合模型的應用特點。

根據其同步性，選取累積位移變化最明顯的ZG326曲線作為白家包滑坡的研究對象，同時為方便插值計算，將日期通過做差轉化為天數。在ZG326曲線中，將85期數據作為原始數據，后5期數據作為指數平滑的訓練數據。由于使用了指數平滑法，為了保證每期數據的等距性，需要通過插值算法將數據變為以30 d為間隔的位移數據。

3.2 演化階段劃分

在1963年，日本學者齋藤迪孝提出了齋藤蠕變曲線模型[18]。如圖4所示，通過觀察齋藤蠕變曲線模型，滑坡形變有3個演變階段，分別為初始形變階段、等速形變階段與加速形變階段。研究成果表明：斜坡的穩定性狀況與其變形階段有著直接的聯系，準確地把握斜坡的變形演化階段是進行斜坡穩定性評價和滑坡預測預報的基礎。從圖4斜坡變形三階段演化模式可以分析得出：在斜坡的初始形變階段，當形變在外界因素的作用下突然啟動后，隨著外界因素減弱甚至消失，其變形速率會逐漸降低，其加速度為負值[19]；在斜坡的等速形變階段，由于其速率基本維持在一恒定值，加速度基本為0；而一旦進入加速形變階段，隨著變形速率的不斷增加，其加速度變為正值，并呈逐漸增大的趨勢，超過一界限即表示滑坡進入一臨滑階段[20]。在本實驗中，忽略初始形變階段，將曲線按照等速形變階段與加速形變階段劃分曲線區間，加速度平穩的等速形變劃為平穩區間，加速度突變的加速形變劃為加速區間，在平穩區間使用高階Lagrange插值法，由于曲線斜率較小，不易發生龍格現象；在加速區間，由于斜率較大，使用低階Lagrange插值法可以保證準確度更高。由于在每年5～9月份降雨增大以及白家包滑坡周期性的規律，往往每年5月以及9月為滑坡形變曲線的拐點，將其作為插值區間的劃分。

圖4 齋藤蠕變曲線模型Fig.4 Saivine creep curve model

3.3 數據插值與預測分析

由于插值算法對曲線的加速度敏感，不同區間需要使用合適的插值算法。為了研究不同插值法對預測精度的影響，使用3種插值方法進行對比，分別為直接插值、全局分段插值和局部分段插值。直接插值是將全部85期數據代入Lagrange多項式，構成85階Lagrange插值計算式的一種插值方法；全局分段插值是全程使用同一分段指標，使用一組全程7階Lagrange分段插值的插值方法；最后一種方法是結合上述的滑坡演化區間，將平穩區間使用10階Lagrange插值模擬高階分段插值，加速區間使用5階Lagrange插值模擬低階分段插值，作為局部分段插值。經過直接插值、全局分段插值與局部分段插值后，結果如表1所列。

表1 插值數據對比

由表1可知，由于整體插值的數據過多，發生了龍格現象，造成了極大的誤差，故舍去此方法。將經過插值的數據代入指數平滑算法中，設平滑指數為0.5，使用剩下5期數據訓練，預測出近期的位移變化。為作對比，將原數據、全局插值后的數據與局部插值后的數據分別代入，具體結果如圖5所示。

圖5 滑坡累積位移插值預測對比Fig.5 Comparison of interpolation prediction of landslide cumulative displacement

圖5中左邊為位移插值對比，其中①、③、⑤、⑦、⑨區間為平穩區間，剩余較陡峭的②、④、⑥、⑧區間為加速區間。通過觀察，在整體演變趨勢上，3條曲線基本吻合滑坡的周期性位移形變規律，表現為在平穩區間前后變化較小，加速區間位移變化明顯。通過具體分析，在加速區間上，全局插值曲線的擬合情況與原始數據較為偏離，而局部插值曲線與原數據曲線更為相似。這兩列數據的具體誤差如表2所列。

表2 插值誤差對比

如表2所列，局部分段插值的誤差相比于全局分段插值在各階段的誤差均小一些，在整體上，局部分段插值的累積誤差相較于傳統的全局分段插值降低了近100 mm。

在圖5中，右邊的6個數據為預測對比，可以發現曲線近期呈現加速上升的趨勢。如表3所列，將3組數據放入指數平滑算法中，計算出第1 978天也就是2012年4月20日的預測值，然后進行5組數據的訓練，最終可以預測出2012年9月19日的累積位移值，在誤差方面，使用局部插值預測在訓練上的誤差較小，整體上，全局插值預測累積誤差為1.33 mm，而局部插值預測累積誤差為0.04 mm，提升了近1 mm。

表3 平滑效果對比

4 結 語

本文針對階躍型滑坡的準確位移預測提出了一種改進的Lagrange插值算法。該新算法在傳統插值算法的基礎上，根據Lagrange多項式的特點增加了先驗分階段的步驟，在不同階段使用不同指標的插值方法。由于前后歷史數據的敏感性不同，加速階段和平穩區間采用不同階數的數據擬合插值，能進一步反映滑坡在演化過程中的動態實時特點。最后通過實驗證明了分段插值的方法誤差更低，形成的位移預測值更精確，相較于傳統方法更具有優勢。本文在插值方法上提出了一個新思路，適用于有明顯演化階段的階躍型滑坡類型，但對于形變曲線規律不明顯的滑坡，找到關鍵拐點的方法還有待改進，為此需進一步研究。