段 功 豪，馬 迅，付 杰

(1.武漢工程大學(xué) 計算機(jī)科學(xué)與工程學(xué)院，湖北 武漢 430205； 2.武漢工程大學(xué) 智能機(jī)器人湖北省重點(diǎn)實驗室，湖北 武漢 430205； 3.中國地質(zhì)調(diào)查局 水文地質(zhì)環(huán)境地質(zhì)調(diào)查中心，河北 保定 071051)

0 引 言

滑坡是一種在山區(qū)或丘陵地區(qū)十分常見的地質(zhì)災(zāi)害。據(jù)自然資源部統(tǒng)計，2021年上半年中國共發(fā)生地質(zhì)災(zāi)害1 150起，造成17人死亡，3人失蹤，直接經(jīng)濟(jì)損失達(dá)24 084.9萬元，其中滑坡高達(dá)579起，因此對滑坡精準(zhǔn)的預(yù)測可以極大程度地減少損失。

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

因此，在對曲線分段時使用什么算法、對哪些點(diǎn)進(jìn)行插值、應(yīng)在曲線上選取哪些點(diǎn)作為分段點(diǎn)，成為學(xué)者們廣泛討論的問題[6-7]?；谏鲜鲅芯勘尘埃疚倪x擇目的性強(qiáng)、靈活度高的Lagrange插值法；為保證最后預(yù)測的準(zhǔn)確性，插值點(diǎn)的選取不僅需解決數(shù)據(jù)的缺值異常值，有時還要確保數(shù)據(jù)的等距性，根據(jù)需求設(shè)置插值點(diǎn)；在插值分段選取上，根據(jù)不同的演化階段，以拐點(diǎn)作為分段點(diǎn)，將曲線分為相應(yīng)的區(qū)間，對每個區(qū)間進(jìn)行局部分段插值，既保證了插值的整體關(guān)聯(lián)性，又能穩(wěn)定插值與當(dāng)前運(yùn)動情況的準(zhǔn)確性；最后對白家包滑坡的累積位移-時間曲線進(jìn)行基于實際運(yùn)動學(xué)變化情況的插值，并對比驗證了該方法的準(zhǔn)確性。

1 研究區(qū)概況

白家包滑坡屬于典型的階躍型滑坡，位于宜昌市秭歸縣歸州鎮(zhèn)向家店村，處于香溪右岸[8]。通過圖1無人機(jī)拍攝可見，周圍地貌部位為低山區(qū)凹形斜坡，前緣直抵香溪河，滑坡剪出口位于高程125～135 m處，滑坡后緣以基巖為界，高程265 m，滑坡左側(cè)以山脊下部基巖為界，右側(cè)以山梁為界，前緣寬500 m，后緣寬300 m，均寬400 m，縱長約550 m，滑坡面積22萬m2?；缕矫嫘螒B(tài)呈短舌狀[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]。此地氣候類型為亞熱帶季風(fēng)氣候，且有明顯的峽谷氣候特征，降水豐沛，四季分明，年平均降雨量為1 000～1 400 mm，降雨多集中于夏季，是暴雨高發(fā)地區(qū)之一。三峽峽谷深邃，風(fēng)力較強(qiáng)，風(fēng)向固定，為東西向的順峽方向?？碧劫Y料顯示，此地滑坡體的土壤物質(zhì)成分由粉質(zhì)黏土混合碎塊石土組成，土質(zhì)松散且雜亂。如圖2滑坡體剖面圖所示，滑床物質(zhì)主要由侏羅系下巴東組的長石石英砂巖及泥巖組成，巖層傾向為260°～285°，傾角為30°～40°[13]?；虑熬壌蟛糠纸]在河水中，當(dāng)水庫水位上升至175.00 m時，部分滑體將完全淹沒于水中。山體植被覆蓋率較好，主要由亞熱帶闊葉林和部分山地灌草叢組成，但在山體上部有部分裸露的土壤巖石，在下部有數(shù)十戶居民，以及幾百畝開墾過的農(nóng)田，其中還有一條秭興公路從中穿過，因此若此地發(fā)生滑坡，將造成巨大的損失。

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

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

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

2 研究方法

2.1 插值算法

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

(1)

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

(2)

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

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

(3)

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

2.2 指數(shù)平滑算法

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

(4)

(5)

xt+T=AT+BTT

(6)

其中，

(7)

(8)

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

3 數(shù)據(jù)分析及預(yù)測分析

3.1 數(shù)據(jù)選取

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

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

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

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

3.2 演化階段劃分

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

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

3.3 數(shù)據(jù)插值與預(yù)測分析

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

表1 插值數(shù)據(jù)對比

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

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

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

表2 插值誤差對比

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

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

表3 平滑效果對比

4 結(jié) 語

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