基于時(shí)間序列與混合核函數(shù)SA-SVR的滑坡位移預(yù)測(cè)模型研究*

尚 敏 熊德兵 張惠強(qiáng) 趙國(guó)飛

(①湖北省地質(zhì)災(zāi)害防治工程技術(shù)研究中心(三峽大學(xué))，宜昌 443002，中國(guó))(②湖北長(zhǎng)江三峽滑坡國(guó)家野外科學(xué)觀測(cè)研究站，宜昌 443002，中國(guó))(③三峽庫(kù)區(qū)地質(zhì)災(zāi)害教育部重點(diǎn)實(shí)驗(yàn)室(三峽大學(xué))，宜昌 443002，中國(guó))(④三峽大學(xué)土木與建筑學(xué)院，宜昌 443002，中國(guó))

0 引 言

滑坡地質(zhì)災(zāi)害是當(dāng)今主要地質(zhì)災(zāi)害之一，對(duì)人民生命財(cái)產(chǎn)造成嚴(yán)重威脅。構(gòu)建合理完善的滑坡位移預(yù)測(cè)模型是實(shí)現(xiàn)滑坡地質(zhì)災(zāi)害監(jiān)測(cè)預(yù)警的關(guān)鍵技術(shù)，能有效減少其給國(guó)家和人民造成的重大損失(尚敏等，2019)。因此，科學(xué)有效地展開滑坡位移預(yù)測(cè)模型研究，具有重要的理論意義和現(xiàn)實(shí)意義。

根據(jù)現(xiàn)有資料，滑坡的預(yù)測(cè)預(yù)報(bào)可追溯到20世紀(jì)60年代日本學(xué)者齋藤(Satio，1969)通過(guò)經(jīng)驗(yàn)公式對(duì)滑坡預(yù)測(cè)的初步探索，其后滑坡預(yù)測(cè)方法經(jīng)歷了從線性預(yù)測(cè)到非線性預(yù)測(cè)的轉(zhuǎn)變(Scheidegger，1973；Murray et al.，1994；韓賀鳴等，2019)。時(shí)間序列分解方法以其完善的理論獲得廣大學(xué)者的青睞，根據(jù)其分解項(xiàng)的復(fù)雜程度不同采用線性或非線性的預(yù)測(cè)方法，大幅度提高了滑坡災(zāi)害預(yù)測(cè)準(zhǔn)確性。該方法主要是基于滑坡位移監(jiān)測(cè)數(shù)據(jù)、降雨量數(shù)據(jù)、庫(kù)水位調(diào)度數(shù)據(jù)等對(duì)滑坡位移進(jìn)行時(shí)間序列分解，然后分別對(duì)分解的子序列應(yīng)用模型進(jìn)行預(yù)測(cè)。目前提取趨勢(shì)項(xiàng)的主要方法有：周超等(2015)應(yīng)用移動(dòng)平均法，將白水河滑坡累積位移分解為周期項(xiàng)和趨勢(shì)項(xiàng)，這種方法操作簡(jiǎn)單，在階躍幅度基本一致的滑坡上能取得不錯(cuò)的效果；黃海峰等(2014)采用平滑先驗(yàn)法經(jīng)過(guò)簡(jiǎn)單計(jì)算快速分離趨勢(shì)項(xiàng)和周期項(xiàng)，對(duì)白家包滑坡進(jìn)行了預(yù)測(cè)，結(jié)果表明，在正則化參數(shù)取值較小時(shí)，其精度高于移動(dòng)平均法但其正則化參數(shù)的取值規(guī)律還有待探究；張凱翔等(2017)應(yīng)用小波變換將滑坡累積位移分解為各個(gè)頻段，低頻分解項(xiàng)即為趨勢(shì)項(xiàng)位移，小波變換具有不需要先驗(yàn)知識(shí)，可以處理具有復(fù)雜變化趨勢(shì)的數(shù)據(jù)等優(yōu)點(diǎn)，但其基函數(shù)和小波階次較難確定；鄧冬梅等(2017)應(yīng)用集合經(jīng)驗(yàn)?zāi)B(tài)分解白水河滑坡位移，具有不需要事先估計(jì)分量的函數(shù)類型且適用性高的特點(diǎn)，但該方法需要對(duì)累積位移序列、降雨量序列和庫(kù)水位分別重構(gòu)，工作量較大，較為復(fù)雜。以上方法各有優(yōu)劣，本文在移動(dòng)平均法的基礎(chǔ)上進(jìn)行改進(jìn)，采用指數(shù)平滑法對(duì)累積位移反復(fù)分解，將最后一次分解的結(jié)果作為趨勢(shì)項(xiàng)，克服了移動(dòng)平均法在累積位移階躍幅度不一致時(shí)分解結(jié)果仍保留明顯階躍性的問(wèn)題。

在周期項(xiàng)的預(yù)測(cè)上，柳青等(2019)將時(shí)間序列和BP神經(jīng)網(wǎng)絡(luò)相結(jié)合，其預(yù)測(cè)位移與實(shí)際位移曲線較一致；湯羅圣等(2012)采用GM(1，1)預(yù)測(cè)八字門滑坡位移，楊背背等(2018)采用長(zhǎng)短時(shí)記憶神經(jīng)網(wǎng)絡(luò)(LSTM)對(duì)三峽庫(kù)區(qū)白水河滑坡位移進(jìn)行預(yù)測(cè)，能更好地反應(yīng)滑坡演化的動(dòng)態(tài)本質(zhì)；然而灰色模型缺乏對(duì)數(shù)據(jù)樣本的學(xué)習(xí)過(guò)程，神經(jīng)網(wǎng)絡(luò)受限于局部最優(yōu)，收斂速度慢等缺點(diǎn)。支持向量機(jī)回歸(SVR)具有小樣本，多維度，非線性的特點(diǎn)，在滑坡預(yù)測(cè)預(yù)報(bào)上獲得廣泛應(yīng)用。針對(duì)SVR參數(shù)尋優(yōu)困難的問(wèn)題，學(xué)者們采用不同的優(yōu)化智能算法進(jìn)行參數(shù)尋優(yōu)，李麟瑋等(2018)采用灰狼優(yōu)化算法(GWO)、楊帆等(2019)采用人工蜂群算法、張俊等(2015)采用粒子群優(yōu)化算法(PSO)、徐峰等(2018)采用自適應(yīng)粒子群算法(AMPSO)；在此基礎(chǔ)上，本文采用模擬退火算法(SA)進(jìn)行參數(shù)尋優(yōu)，通過(guò)線性核與RBF核構(gòu)造混合核函數(shù)，采用構(gòu)造后的混合核函數(shù)對(duì)SA-SVR進(jìn)行優(yōu)化，結(jié)果表明，混合核函數(shù)對(duì)SVR預(yù)測(cè)精度有重要影響。

同時(shí)，滑坡影響因子的選取對(duì)周期項(xiàng)預(yù)測(cè)有重要影響。楊背背等(2018)選擇當(dāng)月降雨量、兩月降雨量、庫(kù)水位高程、當(dāng)月庫(kù)水位、當(dāng)月庫(kù)水位變化、雙月庫(kù)水位變化、當(dāng)月位移增量、前兩月位移增量、前三月位移增量作為周期項(xiàng)預(yù)測(cè)的影響因子，結(jié)合魯芃等(2017)通過(guò)Geo-studio對(duì)滑坡滲流場(chǎng)及應(yīng)力-應(yīng)變進(jìn)行數(shù)值模擬結(jié)果，對(duì)比白家包累積位移、降雨和庫(kù)水位監(jiān)測(cè)曲線，得出白家包滑坡變形對(duì)庫(kù)水位下降更加敏感，變形速率與庫(kù)水位下降速率正相關(guān)的特性。本文在此基礎(chǔ)上，對(duì)庫(kù)水位因子進(jìn)行更加詳細(xì)的劃分，分別為庫(kù)水位、庫(kù)水位變化量、兩月庫(kù)水位變化量、庫(kù)水位下降量、平均庫(kù)水位、平均庫(kù)水位變化量、兩月平均庫(kù)水位變化量和平均庫(kù)水位下降量，經(jīng)過(guò)對(duì)比，前4個(gè)更能反應(yīng)周期項(xiàng)的細(xì)部變化，后4個(gè)更能代表周期項(xiàng)的整體變化趨勢(shì)。

基于上述研究背景，本文以三峽庫(kù)區(qū)白家包滑坡為例，基于時(shí)間序列分解原理，通過(guò)指數(shù)平滑法反復(fù)分解累積位移序列，將最后一次分解所得結(jié)果作為趨勢(shì)項(xiàng)；趨勢(shì)項(xiàng)采用K-flod交叉驗(yàn)證后的3次多項(xiàng)式擬合；周期項(xiàng)采用混合核函數(shù)SA-SVR模型進(jìn)行預(yù)測(cè)；最后，將預(yù)測(cè)的趨勢(shì)項(xiàng)位移和周期項(xiàng)位移疊加，即為滑坡累積位移預(yù)測(cè)結(jié)果，并通過(guò)判定系數(shù)和均方根誤差對(duì)模型進(jìn)行評(píng)價(jià)。

1 滑坡位移預(yù)測(cè)模型

1.1 時(shí)間序列分析

滑坡位移的產(chǎn)生是一個(gè)復(fù)雜的非線性過(guò)程，按時(shí)間序列分解理論可將滑坡位移分解為3個(gè)部分：由滑坡自身地質(zhì)條件(如巖土體類型及成分、地形地貌、地質(zhì)構(gòu)造等)決定的趨勢(shì)項(xiàng)位移，其表現(xiàn)為隨時(shí)間的單調(diào)遞增；由外界誘發(fā)因素(如降雨、庫(kù)水位調(diào)度等)引起的周期項(xiàng)位移，表現(xiàn)為隨時(shí)間呈周期性變化；以及由隨機(jī)事件引起的波動(dòng)項(xiàng)位移(黃顯峰等，2007)?；趯?duì)同類型滑坡的研究，波動(dòng)項(xiàng)的值普遍較小，故本文不對(duì)波動(dòng)項(xiàng)位移進(jìn)行探究，因此時(shí)間序列分解模型可簡(jiǎn)化為：

S(t)=ψ(t)+η(t)

(1)

式中：S(t)為時(shí)間位移序列；ψ(t)為趨勢(shì)項(xiàng)位移；η(t)為周期項(xiàng)位移。

由于白家包滑坡不同年份階躍幅度變化較大，采用移動(dòng)平均法分解所得的趨勢(shì)項(xiàng)位移仍保留明顯的階躍性，且得到的周期項(xiàng)較難預(yù)測(cè)。因此，本文反復(fù)使用指數(shù)平滑法對(duì)累積位移進(jìn)行分解，將第N次分解的趨勢(shì)項(xiàng)作為最終的趨勢(shì)項(xiàng)，來(lái)獲得可以代表滑坡趨勢(shì)的位移曲線。

指數(shù)平滑法是C.C Holt于1958年提出的，在滑坡預(yù)測(cè)中廣泛應(yīng)用(尹光志等，2007)。指數(shù)平滑法是一種特殊的加權(quán)平均法，是一種對(duì)本期觀測(cè)值和上期平滑值賦予不同的權(quán)重，求得本期平滑值的方法。指數(shù)平滑法基本原理如下：

(2)

式中：St為第t期平滑值，也即第t期的趨勢(shì)項(xiàng)；yt為第t期觀測(cè)值；St-1為第t-1期平滑值；a為平滑指數(shù)。由式(2)可知，當(dāng)a越接近1時(shí)，第t期平滑值越接近第t期觀察值。

只使用一次指數(shù)平滑法與移動(dòng)平均法一樣難以較好分離趨勢(shì)項(xiàng)和周期項(xiàng)，因此需要多次使用指數(shù)平滑法，其步驟為：

(1)采用指數(shù)平滑法將滑坡累積位移分解，獲得趨勢(shì)項(xiàng)1。

(2)將趨勢(shì)項(xiàng)1用指數(shù)平滑法繼續(xù)分解，獲得趨勢(shì)項(xiàng)2。

(3)多次重復(fù)上述步驟，可獲得趨勢(shì)項(xiàng)N，此時(shí)周期項(xiàng)即為累積位移與趨勢(shì)項(xiàng)N之差，當(dāng)趨勢(shì)項(xiàng)較平滑且周期項(xiàng)有較好預(yù)測(cè)效果時(shí)，分解結(jié)束。

1.2 混合核函數(shù)支持向量機(jī)回歸

支持向量機(jī)回歸(Support Vector Regression，SVR)是Vapnik于20世紀(jì)90年代提出的一種非線性問(wèn)題回歸預(yù)測(cè)方法(Cortes，1995)。其主要原理是：將樣本數(shù)據(jù)劃分為訓(xùn)練樣本和預(yù)測(cè)樣本，通過(guò)非線性映射函數(shù)φ(x)將訓(xùn)練樣本映射到高維空間，使得低維空間的非線性問(wèn)題轉(zhuǎn)化為高維空間的線性問(wèn)題，然后在高維空間中利用核函數(shù)替換映射函數(shù)的點(diǎn)積運(yùn)算，獲得具有最優(yōu)擬合效果的決策函數(shù)模型。該決策函數(shù)即為SVR回歸預(yù)測(cè)模型，通過(guò)預(yù)測(cè)樣本對(duì)決策函數(shù)進(jìn)行驗(yàn)證。SVR回歸預(yù)測(cè)模型原理如下：

針對(duì)需要訓(xùn)練的樣本數(shù)據(jù)，假定訓(xùn)練樣本集為Ω={(x1，y1)，(x2，y2)，…，(xm，ym)}，其中，xi為輸入向量，yi為輸出向量，則SVR回歸預(yù)測(cè)模型的決策函數(shù)可表示為：

f(x)=ωTφ(x)+b

(3)

式中：ωT為權(quán)重系數(shù)；φ(t)為非線性映射函數(shù)；b為偏置量。

求解決策函數(shù)的過(guò)程可看作求解最小化問(wèn)題：

(4)

(5)

(6)

其中：K(xi，xj)=φ(xi)Tφ(xj)為核函數(shù)，通過(guò)核函數(shù)可以在不確定非線性映射函數(shù)φ(x)形式的情況下求得決策函數(shù)，極大地提高了SVR模型計(jì)算效率，常見的核函數(shù)如下：

其中：線性核泛化能力強(qiáng)、學(xué)習(xí)能力弱，它擁有越是遠(yuǎn)離測(cè)試點(diǎn)對(duì)核函數(shù)越有影響的特性，是一種全局核，RBF核泛化能力弱、學(xué)習(xí)能力強(qiáng)，它擁有越是接近測(cè)試點(diǎn)對(duì)核函數(shù)越有影響的特性。此外，核函數(shù)還具有以下性質(zhì)(劉琰，2012)：

(1)若k1為核函數(shù)，λ為常數(shù)，則λk1也是核函數(shù)；

(2)若k1和k2為核函數(shù)，則對(duì)于任意正數(shù)γ1和γ2，其線性組合γ1k1+γ2k2也是核函數(shù)；

(3)若k1和k2為核函數(shù)，則核函數(shù)直積k1k2也是核函數(shù)；

本文通過(guò)性質(zhì)(2)構(gòu)造混合核函數(shù)，將泛化能力強(qiáng)、學(xué)習(xí)能力弱的線性核和泛化能力弱、學(xué)習(xí)能力強(qiáng)的RBF核函數(shù)組合在一起，構(gòu)造成新的核函數(shù)，其形式為：

Kmix=λKploy+(1-λ)Krbf，λ∈[0，1]

(7)

式中：Kmix為混合核函數(shù)；Kploy為線性核函數(shù)；Krbf為RBF核函數(shù)；λ為線性核比例。

1.3 模擬退火算法

模擬退火算法(Simulated Annealing，SA)是美國(guó)物理學(xué)家S.Kirkpatrick于1983年以Metropolis方法為基礎(chǔ)用來(lái)解決組合優(yōu)化問(wèn)題中最優(yōu)解問(wèn)題的一套算法(康立山等，1994)。模擬退火算法模擬了熱力學(xué)中固體退火過(guò)程中，固體在降溫后，內(nèi)能會(huì)降到最低，達(dá)到一種穩(wěn)定狀態(tài)。模擬退火算法是一個(gè)“產(chǎn)生新解-判斷-接受/舍棄”的過(guò)程，其全局尋優(yōu)能力較強(qiáng)，理論上能以概率1收斂于全局最優(yōu)解。在尋優(yōu)初期，能以很大的概率接受一個(gè)較差的解從而避免過(guò)早尋到局部最優(yōu)解而結(jié)束尋優(yōu)過(guò)程；在尋優(yōu)后期，這個(gè)概率會(huì)逐漸減小，從而保證模擬退火算法在后期時(shí)的收斂能力(李元香等，2019)。以求解f(x)最小值為例，模擬退火算法詳細(xì)步驟如下：

(1)初始溫度為T0，馬爾可夫鏈長(zhǎng)度為L(zhǎng)0，在可行解中隨機(jī)挑選一個(gè)初始解x0，此時(shí)迭代次數(shù)i=0，計(jì)算目標(biāo)函數(shù)f(x0)。

(2)用攝動(dòng)裝置產(chǎn)生一次隨機(jī)擾動(dòng)，在可行解空間中產(chǎn)生一個(gè)新解xi，計(jì)算目標(biāo)函數(shù)f(xi)。

(3)用Metropolis準(zhǔn)則判斷是否接受新解：若f(xi)≤f(x0)時(shí)，接受新解xi，此時(shí)最優(yōu)解xp=xi；若f(xi)>f(x0)，則以概率pt接受新解xi，當(dāng)pt>random(0，1)時(shí)，最優(yōu)解為f(xi)，否則最優(yōu)解仍為f(x0)。

(4)重復(fù)步驟2、3，得到鏈長(zhǎng)為L(zhǎng)0過(guò)程的最優(yōu)解。

(5)判斷是否滿足停止條件，若滿足則終止；若不滿足則開始下一次迭代，溫度為Ti+1=μTi，馬爾可夫鏈長(zhǎng)度為L(zhǎng)i+1，回到步驟2。

混合核函數(shù)SA-SVR預(yù)測(cè)周期項(xiàng)流程圖如圖1所示。

圖1 混合核函數(shù)SA-SVR位移預(yù)測(cè)技術(shù)路線

1.4 評(píng)價(jià)指標(biāo)

判定系數(shù)(coefficient of determination，R2)和均方根誤差(root mean square，RMSE)在滑坡的位移預(yù)測(cè)模型中常作為精度評(píng)價(jià)指標(biāo)，因此，使用它們來(lái)評(píng)價(jià)預(yù)測(cè)模型，其定義為：

(8)

(9)

2 白家包滑坡變形特征分析

2.1 白家包滑坡工程地質(zhì)概況

白家包滑坡位于湖北省秭歸縣歸州鎮(zhèn)，位于香溪河右岸，距香溪河入江口2.5 km。白家包滑坡為土質(zhì)滑坡，坡體為逆向結(jié)構(gòu)，其前緣直抵香溪河，滑坡剪出口位于約高程130 m處，滑坡后緣以基巖為界，高程265 m，滑坡左側(cè)以山脊下部基巖為界，右側(cè)以山梁為界，前緣寬500 m，后緣寬300 m，縱長(zhǎng)約550 m，滑坡面積2.2×105m2。

白家包滑坡物質(zhì)為崩坡積物，坡積物厚度在空間上分布不均，前緣坡積物厚度較厚，滑體物質(zhì)成分為碎塊石土?；矠榛鶐r，成分為長(zhǎng)石石英砂巖、泥巖，產(chǎn)狀285°∠30°，滑坡平面形態(tài)呈短舌狀(尚敏等，2021)。

2.2 白家包滑坡位移分析

白家包滑坡上目前布設(shè)ZG323、ZG324、ZG325和ZG326共4個(gè)GPS監(jiān)測(cè)點(diǎn)，于2006年9月開始監(jiān)測(cè)，各監(jiān)測(cè)點(diǎn)位置如圖2所示，監(jiān)測(cè)位移圖像如圖3所示。ZG325監(jiān)測(cè)點(diǎn)位于滑坡中上部，ZG324監(jiān)測(cè)點(diǎn)位于滑坡中下部，且位移曲線相較于ZG323和ZG326更具有代表性，能夠更好地代表白家包滑坡的整體位移過(guò)程，如圖2和圖3所示。鑒于篇幅有限，因此本文選取ZG325監(jiān)測(cè)點(diǎn)進(jìn)行詳細(xì)分析，以ZG324監(jiān)測(cè)點(diǎn)進(jìn)行對(duì)比分析，ZG325監(jiān)測(cè)點(diǎn)累積位移、降雨量和庫(kù)水位監(jiān)測(cè)曲線如圖4所示。顯然，白家包滑坡位移隨時(shí)間呈階躍式增長(zhǎng)，且每間隔幾年會(huì)有一次較大的階躍增長(zhǎng)。根據(jù)庫(kù)水位的調(diào)度情況，白家包滑坡位移可劃分為兩個(gè)階段：

圖2 白家包滑坡工程地質(zhì)平面圖

圖3 白家包滑坡累積位移曲線圖

圖4 ZG325累積位移-降雨量-庫(kù)水位監(jiān)測(cè)曲線

(2)2009年4月至2020年9月，庫(kù)水位峰值由155 m上升到175 m，波動(dòng)范圍為145～175 m，每年8～12月，庫(kù)水位逐漸從波谷上升到波峰，1～5月，庫(kù)水位從波峰下降到波谷，6～7月庫(kù)水位先小幅上升再下降，滑坡在每年5～9月發(fā)生階躍性位移，其中2009年、2011年、2012年、2015年和2016年階躍較大，分別為178.4 mm、112.1 mm、170.5 mm、182.7 mm、150.4 mm。

綜上所述，白家包滑坡為蠕滑-拉裂型滑坡，其變形主要受庫(kù)水位升降的影響，降雨僅對(duì)滑坡變形起推動(dòng)作用(鄧茂林等，2020)，并且?guī)焖缓徒涤陮?duì)滑坡的影響是不同步的，具有明顯的滯后效應(yīng)。如圖4所示，在8～12月期間，月降雨量逐漸減小，庫(kù)水位處于上升階段，又因白家包滑坡坡體黏土滲透性差，坡體受到與坡面正交且指向坡面的滲透水壓力，減緩坡體下滑。在1～7月期間，月降雨量逐漸增加，庫(kù)水位處于下降狀態(tài)，此時(shí)水向坡體外滲透速率遠(yuǎn)小于庫(kù)水位下降速率，產(chǎn)生指向坡外的滲透水壓力，加劇坡體下滑。值得注意的是，滲透水壓力對(duì)坡體的響應(yīng)在時(shí)間上存在一定的滯后，8～9月庫(kù)水位上升時(shí)產(chǎn)生的階躍位移為庫(kù)水位下降期間的滯后影響造成的。

3 白家包滑坡位移預(yù)測(cè)

本文取2009年1月～2020年9月的GPS監(jiān)測(cè)數(shù)據(jù)作為原始數(shù)據(jù)，但多次使用指數(shù)平滑法提取趨勢(shì)項(xiàng)位移的過(guò)程中，會(huì)損失一部分?jǐn)?shù)據(jù)，因此實(shí)際分析時(shí)間段為2012年1月～2020年9月。取2012年1月～2019年9月的監(jiān)測(cè)數(shù)據(jù)作為訓(xùn)練樣本集，以2019年10月～2020年9月的監(jiān)測(cè)數(shù)據(jù)作為預(yù)測(cè)集。

3.1 趨勢(shì)項(xiàng)位移提取和預(yù)測(cè)

本文采用指數(shù)平滑法進(jìn)行趨勢(shì)項(xiàng)提取，其中平滑指數(shù)a的取值不同分解效果不同，不同a的取值分解效果如圖5所示。當(dāng)a值越趨近于1時(shí)，曲線越接近于原始位移，a值越趨近于0時(shí)，曲線越平滑。經(jīng)過(guò)測(cè)試，當(dāng)a取0.3～0.6時(shí)，周期項(xiàng)的預(yù)測(cè)效果較好，由圖5可知，a在0.3～0.6范圍內(nèi)時(shí)，提取出的趨勢(shì)項(xiàng)仍然呈現(xiàn)明顯的階梯狀，因此本文通過(guò)多次使用指數(shù)平滑法對(duì)原始數(shù)據(jù)進(jìn)行分解，分解次數(shù)越多，提取出的趨勢(shì)項(xiàng)越平滑，但是周期項(xiàng)的預(yù)測(cè)效果也會(huì)逐漸降低，因此綜合考慮預(yù)測(cè)精度和提取效率，取a=0.35，分3次提取趨勢(shì)項(xiàng)，其過(guò)程為：以G325的累積位移作為原始數(shù)據(jù)，應(yīng)用指數(shù)平滑法提取趨勢(shì)項(xiàng)，將所得的數(shù)據(jù)項(xiàng)作為新的原始數(shù)據(jù)，再次應(yīng)用指數(shù)平滑法，循環(huán)此步驟直到第3次提取，其中平滑初值為前12期數(shù)據(jù)的算數(shù)平均值。每次趨勢(shì)項(xiàng)提取結(jié)果如圖6所示。

圖5 不同a值的分解效果圖

圖6 趨勢(shì)項(xiàng)分解結(jié)果

將第3次提取的趨勢(shì)項(xiàng)作為白家包滑坡的最終趨勢(shì)項(xiàng)，根據(jù)圖6結(jié)果將該曲線分為兩段(第1階段：2012年1月～2015年7月，第2階段：2015年8月～2019年9月)，分別用3次多項(xiàng)式進(jìn)行擬合，3次多項(xiàng)式通項(xiàng)公式為：

y=ax3+bx2+cx+d

(10)

直接對(duì)第2階段擬合會(huì)導(dǎo)致第2階段的3次多項(xiàng)式在預(yù)測(cè)集嚴(yán)重不適應(yīng)，這是由于過(guò)擬合造成多項(xiàng)式曲率偏小或偏大造成的。為了解決這個(gè)問(wèn)題，本文在第2階段采用K-flod交叉驗(yàn)證的思想，將第2階段的趨勢(shì)項(xiàng)位移隨機(jī)分為5等份(黃光華等，2020)，不重復(fù)的每次取其中一份作為測(cè)試集，其余4份用來(lái)擬合曲線，保留每次擬合曲線的a、b、c、d參數(shù)并計(jì)算該模型在測(cè)試集上的RMSE，計(jì)算5次RMSE、a、b、c、d的平均值，多次重復(fù)上述步驟，最終取RMSE最小的那一組參數(shù)進(jìn)行擬合。在1000次試驗(yàn)中，獲得最佳K-flod交叉驗(yàn)證RMSE為7.331，3次多項(xiàng)式參數(shù)如表1第2階段①所示，用該3次多項(xiàng)式參數(shù)預(yù)測(cè)2019年10月～2020年9月這12期的趨勢(shì)項(xiàng)位移，獲得其RMSE為7.669，預(yù)測(cè)圖像如圖7a所示。在非階躍段，預(yù)測(cè)絕對(duì)誤差逐漸減小，在2020年2月達(dá)到最小1.120 mm，在階躍段絕對(duì)誤差逐漸增大，最大值為13.834 mm。

圖7 趨勢(shì)項(xiàng)位移預(yù)測(cè)結(jié)果

通過(guò)上述分析可知，在一次預(yù)測(cè)12期時(shí)，由于缺乏最近的數(shù)據(jù)對(duì)多項(xiàng)式參數(shù)進(jìn)行校正，隨著預(yù)測(cè)期數(shù)越多，預(yù)測(cè)值漸漸偏離真實(shí)值，在實(shí)際預(yù)測(cè)中可采用每次預(yù)測(cè)一期的方法，將上一期的真實(shí)值帶入訓(xùn)練集校正多項(xiàng)式參數(shù)，預(yù)測(cè)參數(shù)如表1第2階段②。結(jié)果表明，在階躍階段，預(yù)測(cè)的絕對(duì)誤差呈逐漸減小的趨勢(shì)，其中2020年9月的絕對(duì)誤差由13.834 mm下降到6.491 mm，而在2019年10月～2020年4月絕對(duì)誤差下降不大，這說(shuō)明階躍階段的數(shù)據(jù)對(duì)校正起關(guān)鍵作用。

表1 趨勢(shì)項(xiàng)位移多項(xiàng)式預(yù)測(cè)參數(shù)

3.2 周期項(xiàng)位移提取和預(yù)測(cè)

3.2.1 周期項(xiàng)位移的提取

將滑坡ZG325累積位移減去第3次使用指數(shù)平滑法提取的趨勢(shì)項(xiàng)，即可獲得周期項(xiàng)。如圖8所示，所獲得的周期項(xiàng)位移在2012年7月和2015年7月有較大增幅，其中2012年7月到2015年6月振幅逐漸減小，到2015年7月振幅突然增大，其后振幅逐漸減小，該規(guī)律與ZG325累積位移的階梯性保持一致。

圖8 周期項(xiàng)位移提取值

3.2.2 影響因子選取

影響因子的選取對(duì)周期項(xiàng)預(yù)測(cè)起關(guān)鍵作用，好的影響因子能反應(yīng)預(yù)測(cè)曲線的走勢(shì)。白家包滑坡累積位移序列呈明顯的階躍性，這是由于庫(kù)水位的周期調(diào)動(dòng)和降雨量的聯(lián)合作用造成的。除此之外，由3.2.1分析可知，累積位移的階躍性波動(dòng)也能在一定程度上反應(yīng)周期項(xiàng)的變化趨勢(shì)。

3.2.2.1 庫(kù)水位

庫(kù)水位調(diào)度會(huì)改變坡體地下水滲透壓力。在庫(kù)水位下降期間，由于白家包滑坡滲透系數(shù)遠(yuǎn)小于庫(kù)水位的下降速度，因此產(chǎn)生指向坡外的動(dòng)水壓力，不利于滑坡穩(wěn)定。本文選取庫(kù)水位、庫(kù)水位變化量、兩月庫(kù)水位變化量、庫(kù)水位下降量、平均庫(kù)水位、平均庫(kù)水位變化量、兩月平均庫(kù)水位變化量和平均庫(kù)水位下降量這8個(gè)影響因子。其中平均庫(kù)水位、平均庫(kù)水位變化量、兩月平均庫(kù)水位變化量和平均庫(kù)水位下降量更能總體反應(yīng)周期項(xiàng)的變化趨勢(shì)，而庫(kù)水位、庫(kù)水位變化量、兩月庫(kù)水位變化量和庫(kù)水位下降量則更能從反應(yīng)周期項(xiàng)的細(xì)部變化，周期項(xiàng)和庫(kù)水位關(guān)系如圖9a。

如圖9b，對(duì)比當(dāng)月平均庫(kù)水位下降量、兩月前平均庫(kù)水位下降量與周期項(xiàng)位移曲線可知，兩月前平均庫(kù)水位下降量的峰值與周期項(xiàng)位移峰值基本保持一致，這不僅說(shuō)明庫(kù)水位月下降速率對(duì)周期項(xiàng)位移有重要影響，而且表明庫(kù)水位對(duì)周期項(xiàng)的影響在時(shí)間上滯后兩月。且由圖可知，庫(kù)水位在1～5月逐漸從波峰下降到波谷，6～7月庫(kù)水位先小幅上升再下降，而滑坡在每年5～9月發(fā)生階躍性位移，1～4月累積位移基本保持不變，也即庫(kù)水位下降沒有立即使1～4月累積位移發(fā)生階躍性變化，這說(shuō)明庫(kù)水位對(duì)累積位移的影響還存在一個(gè)閾值，只有低于這個(gè)閾值時(shí)才會(huì)產(chǎn)生階躍性位移。從庫(kù)水位的調(diào)度情況來(lái)看，當(dāng)庫(kù)水位從175 m下降到165 m左右時(shí)，其后約兩個(gè)月，滑坡累積位移將逐漸發(fā)生階躍性躍遷。由于庫(kù)水位影響因子中有兩月因子，考慮到庫(kù)水位對(duì)滑坡位移的綜合影響，因此本文將庫(kù)水位影響因子延后一月，即以一月前的庫(kù)水位影響因子對(duì)當(dāng)前月位移進(jìn)行預(yù)測(cè)。降雨和位移增量對(duì)周期項(xiàng)位移的影響在時(shí)間上和庫(kù)水位保持一致，因此將降雨因子和位移增量因子都延后一月。

3.2.2.2 降 雨

根據(jù)馬銳(2019)對(duì)白家包滑坡變形與降雨、庫(kù)水位相關(guān)性定量化分析研究結(jié)果，白家包滑坡的變形主要是庫(kù)水位快速下降造成的，因此本文僅選取降雨量和降雨變化量作為降雨因子。降雨一定程度上會(huì)使土體軟化，降低巖土體強(qiáng)度，隨著降雨量增加，坡體自重逐漸增大，從而引起坡體位移速率增大(吳仁銑，2013)。如圖9c可知，強(qiáng)降雨和周期項(xiàng)位移峰值在時(shí)間上基本保持同步。

3.2.2.3 位移增量

滑坡累積位移-時(shí)間曲線呈階躍性，并且其階躍的高度與周期項(xiàng)位移的振幅基本保持一致，說(shuō)明滑坡累積位移月增長(zhǎng)速率和周期項(xiàng)位移有一定的聯(lián)系，如圖9c所示。因此將第2次分解后趨勢(shì)項(xiàng)的位移增量作為影響因子。

圖9 周期項(xiàng)位移與影響因子關(guān)系

3.2.3 混合核函數(shù)SA-SVR模型預(yù)測(cè)過(guò)程

3.2.3.1 數(shù)據(jù)預(yù)處理

為了消除不同量綱對(duì)模型的影響，需要將樣本數(shù)據(jù)進(jìn)行歸一化處理。本文采用Min-Max歸一化方法，將影響因子和周期項(xiàng)位移歸一化到[0，1]區(qū)間內(nèi)。

3.2.3.2 構(gòu)造混合核函數(shù)

本文通過(guò)線性核函數(shù)與高斯徑向基核函數(shù)的線性組合構(gòu)造混合核函數(shù)，通過(guò)調(diào)整線性指數(shù)的數(shù)值，以期獲得泛化能力和學(xué)習(xí)能力都較強(qiáng)的SVR模型。試驗(yàn)證明，不同懲罰系數(shù)C和核函數(shù)參數(shù)g下的線性核比例-預(yù)測(cè)結(jié)果RMSE曲線的差別巨大，如圖10所示，在參數(shù)C=8.15，g=0.89時(shí)，預(yù)測(cè)結(jié)果均方根誤差先減小，然后逐漸增大；而在參數(shù)C=6，g=3時(shí)，預(yù)測(cè)結(jié)果均方根誤差隨著線性指數(shù)先減小，再緩慢增大，且最大值與最小值差距較大。這兩種情況充分說(shuō)明線性指數(shù)對(duì)預(yù)測(cè)精確度有重大影響。本文以線性核比例控制混合核函數(shù)作為SVR模型的優(yōu)化參數(shù)之一，通過(guò)模擬退火算法對(duì)線性核比例進(jìn)行參數(shù)尋優(yōu)。混合核函數(shù)實(shí)現(xiàn)采用MATLAB分析軟件，由于LIBSVM工具箱沒有現(xiàn)成的混合核函數(shù)可以調(diào)用，因此需要自己構(gòu)建混合核函數(shù)。

圖10 線性核比例與預(yù)測(cè)結(jié)果RMSE曲線圖

3.2.3.3 模擬退火算法參數(shù)設(shè)定及目標(biāo)函數(shù)構(gòu)造

設(shè)定初始溫度為100，溫度衰減系數(shù)為0.9，最大迭代次數(shù)為200，每個(gè)溫度下迭代次數(shù)為150，將懲罰系數(shù)C、核函數(shù)參數(shù)g和線性核比例λ構(gòu)成一個(gè)解空間，并設(shè)定其上限為[20，10，1]，下限為[1，0，0]，以SVR交叉驗(yàn)證的均方根誤差作為目標(biāo)函數(shù)進(jìn)行參數(shù)尋優(yōu)。

3.2.3.4 混合核函數(shù)SA-SVR模型預(yù)測(cè)

本文取2012年1月～2019年9月樣本作為訓(xùn)練集，2019年10月～2020年9月樣本作為預(yù)測(cè)集，應(yīng)用MATLAB 2018a數(shù)學(xué)分析軟件實(shí)現(xiàn)模型預(yù)測(cè)。將歸一化的11個(gè)影響因子作為SVM輸入向量，周期項(xiàng)位移作為輸出向量，SVM固定參數(shù)采用LIBSVM工具箱中的默認(rèn)值，懲罰系數(shù)C、核函數(shù)參數(shù)g和線性核比例λ使用模擬退火算法進(jìn)行參數(shù)尋優(yōu)，核函數(shù)使用混合核函數(shù)，在訓(xùn)練集中獲得的最佳參數(shù)組合為C=8.154，g=0.891，λ=0.460。將該參數(shù)帶入預(yù)測(cè)集中進(jìn)行預(yù)測(cè)，周期項(xiàng)位移預(yù)測(cè)結(jié)果如圖11所示。混合核函數(shù)SA-SVR預(yù)測(cè)模型的判定系數(shù)R2和均方根誤差RMSE分別為0.791和4.9348，由圖11可知，在2019年10月～2020年4月期間，由于累積位移基本保持不變，因此周期項(xiàng)表現(xiàn)為逐漸下降的趨勢(shì)，在2020年5月～2020年9月，由于庫(kù)水位下降和持續(xù)降雨，滑坡累積位移發(fā)生階躍性躍遷，因此周期項(xiàng)表現(xiàn)為逐漸增大的趨勢(shì)。模型在2019年10月～2020年4月表現(xiàn)良好，平均絕對(duì)誤差為2.648 mm，RMSE為3.075。在2020年7月～2020年8月，模型在峰值處預(yù)測(cè)相較于真實(shí)值略高，絕對(duì)誤差分別為9.99 mm和4.39 mm，表現(xiàn)為預(yù)測(cè)的結(jié)果更加保守，結(jié)合2019年10月～2019年12月預(yù)測(cè)結(jié)果可知，模型在峰值后下降速度較真實(shí)值更快。

圖11 ZG325監(jiān)測(cè)點(diǎn)周期項(xiàng)位移預(yù)測(cè)結(jié)果

3.3 滑坡累積位移預(yù)測(cè)

將趨勢(shì)項(xiàng)位移預(yù)測(cè)和周期項(xiàng)位移預(yù)測(cè)相加即為累積位移預(yù)測(cè)，預(yù)測(cè)結(jié)果如圖12所示。由圖可知，預(yù)測(cè)結(jié)果和真實(shí)值基本保持一致，預(yù)測(cè)階段判定系數(shù)R2為0.790，均方根誤差RMSE為8.5877，在階躍時(shí)預(yù)測(cè)結(jié)果較真實(shí)值略高，其原因是周期項(xiàng)和趨勢(shì)項(xiàng)在階躍階段預(yù)測(cè)結(jié)果比真實(shí)值都偏高。如表2所示，2020年7月出現(xiàn)最大絕對(duì)誤差，其值為15.831 mm，此時(shí)相對(duì)誤差為1.086%，最小絕對(duì)誤差在2020年3月，其值為1.699 mm，此時(shí)相對(duì)誤差為0.119%。從整體上看，預(yù)測(cè)結(jié)果在6月會(huì)提前靠近階躍峰值，有一定的預(yù)警效果，預(yù)測(cè)結(jié)果在階躍段略高于真實(shí)值，然后在非階躍段逐漸回落至真實(shí)值，因此在實(shí)際應(yīng)用中具有一定的安全保障。

圖12 ZG325監(jiān)測(cè)點(diǎn)累積位移預(yù)測(cè)結(jié)果

表2 累積位移預(yù)測(cè)精度

4 在監(jiān)測(cè)點(diǎn)ZG324上的驗(yàn)證

將基于時(shí)間序列分解和SA-SVR混合核函數(shù)預(yù)測(cè)模型應(yīng)用到ZG324監(jiān)測(cè)點(diǎn)上進(jìn)行輔助驗(yàn)證，采用K-flod交叉驗(yàn)證思想的3次多項(xiàng)式預(yù)測(cè)，在2019年10月～2020年9月期間，預(yù)測(cè)結(jié)果的絕對(duì)誤差呈現(xiàn)先減小后增大的趨勢(shì)，最大絕對(duì)誤差發(fā)生在12月，其值為9.014 mm，最小絕對(duì)誤差在4月，其值為0.689 mm，均方根誤差RMSE為6.532。通過(guò)MATLAB分析軟件的LIBSVM工具箱對(duì)樣本進(jìn)行訓(xùn)練并通過(guò)SA算法對(duì)參數(shù)尋優(yōu)，其最佳參數(shù)為C=9.963，g=0.957，λ=0.950，最大絕對(duì)誤差為10.574 mm，最小絕對(duì)誤差為0.038 mm，均方根誤差RMSE為4.559，預(yù)測(cè)結(jié)果在2019年10月～2020年4月效果較好，在2020年5月～2020年8月預(yù)測(cè)值略高于真實(shí)值，并于2020年7月達(dá)到峰值，隨后預(yù)測(cè)值逐漸回落，該趨勢(shì)與ZG325監(jiān)測(cè)點(diǎn)基本保持一致，ZG324監(jiān)測(cè)點(diǎn)趨勢(shì)項(xiàng)位移與周期項(xiàng)位移預(yù)測(cè)結(jié)果如圖13所示。

圖13 ZG324監(jiān)測(cè)點(diǎn)趨勢(shì)項(xiàng)位移與周期項(xiàng)位移預(yù)測(cè)結(jié)果

疊加趨勢(shì)項(xiàng)位移和周期項(xiàng)位移預(yù)測(cè)結(jié)果獲得累積位移預(yù)測(cè)，如圖14所示。預(yù)測(cè)階段最大絕對(duì)誤差為15.598 mm，最小絕對(duì)誤差為0.722 mm，均方根誤差RMSE為8.445，并且預(yù)測(cè)結(jié)果在階躍階段較真實(shí)值略高，有一定的安全保障，該過(guò)程與ZG325監(jiān)測(cè)點(diǎn)基本一致，說(shuō)明本模型具有較強(qiáng)的推廣能力。

圖14 ZG324監(jiān)測(cè)點(diǎn)累積位移預(yù)測(cè)結(jié)果

5 結(jié) 論

本文結(jié)合時(shí)間序列分解和SA-SVR混合核函數(shù)預(yù)測(cè)模型，對(duì)白家包滑坡ZG325監(jiān)測(cè)點(diǎn)累積位移進(jìn)行預(yù)測(cè)，預(yù)測(cè)結(jié)果與真實(shí)值基本一致，其中平均絕對(duì)誤差為4.016 mm，最大絕對(duì)誤差為15.830 mm，最小絕對(duì)誤差為1.699 mm，判定系數(shù)R2為0.790，均方根誤差RMSE為8.5877，隨后采用ZG324監(jiān)測(cè)點(diǎn)進(jìn)行輔助驗(yàn)證，取得了很好的效果，可以作為同類滑坡監(jiān)測(cè)預(yù)警參考模型。

同時(shí)，基于本文研究，時(shí)間序列與混合核函數(shù)SA-SVR的滑坡位移預(yù)測(cè)模型和傳統(tǒng)的時(shí)間序列SVR預(yù)測(cè)模型相比，具有以下優(yōu)點(diǎn)：

(1)在分解方法上摒棄了傳統(tǒng)的移動(dòng)平均法，采用更適應(yīng)白家包滑坡階躍特點(diǎn)的指數(shù)平滑分解法，并反復(fù)使用指數(shù)平滑法使得分解的趨勢(shì)項(xiàng)和周期項(xiàng)更符合客觀規(guī)律。

(2)在趨勢(shì)項(xiàng)預(yù)測(cè)方法上，采用K-fold交叉驗(yàn)證的多項(xiàng)式參數(shù)進(jìn)行預(yù)測(cè)，有效緩解了多項(xiàng)式過(guò)擬合造成預(yù)測(cè)值和實(shí)際值嚴(yán)重偏離的問(wèn)題。

(3)針對(duì)多參數(shù)尋優(yōu)問(wèn)題，采用更加高效的模擬退火算法。

(4)在核函數(shù)的選擇上，不僅僅使用學(xué)習(xí)能力較強(qiáng)的RBF核函數(shù)，而是使用結(jié)合RBF核函數(shù)和線性核函數(shù)的混合核函數(shù)，結(jié)果表明，混合核函數(shù)對(duì)周期項(xiàng)預(yù)測(cè)精度有重要影響。

(5)通過(guò)分析誘發(fā)因子對(duì)周期項(xiàng)位移的滯后關(guān)系，以當(dāng)期影響因子預(yù)測(cè)后一期周期項(xiàng)位移，從而實(shí)現(xiàn)1月內(nèi)短期位移實(shí)時(shí)預(yù)測(cè)，具有較高實(shí)用價(jià)值。

(6)從預(yù)測(cè)結(jié)果上來(lái)看，該模型預(yù)測(cè)結(jié)果在6月會(huì)提前靠近階躍峰值，有一定的預(yù)警效果，在非階躍段預(yù)測(cè)結(jié)果與真實(shí)值基本保持一致，而階躍段的預(yù)測(cè)結(jié)果略高于真實(shí)值，然后逐漸回落至與下一階段非階躍段高度，具有一定的安全保障。

本模型是建立在白家包多年的滑坡累積位移、降雨量、三峽庫(kù)區(qū)庫(kù)水位等監(jiān)測(cè)數(shù)據(jù)上的，由于采用指數(shù)平滑法分解的過(guò)程中不可避免地?fù)p失部分?jǐn)?shù)據(jù)，因此該方法在數(shù)據(jù)較少的滑坡上具有局限性。此外，在參數(shù)尋優(yōu)過(guò)程中，訓(xùn)練參數(shù)過(guò)多會(huì)影響計(jì)算時(shí)間且占用更多內(nèi)存，參數(shù)設(shè)置不合理或訓(xùn)練次數(shù)不足也可能導(dǎo)致模型輸出結(jié)果不穩(wěn)定。