羅袆沅，蔣亞楠,，許 強(qiáng)，唐 斌

(1.成都理工大學(xué) 地球科學(xué)學(xué)院,四川 成都 610059；2.成都理工大學(xué) 地質(zhì) 災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點(diǎn)實(shí)驗(yàn)室,四川 成都 610059)

0 引言

GPS的靜態(tài)觀測作為滑坡變形監(jiān)測的常用方法，其變形監(jiān)測數(shù)據(jù)是滑坡預(yù)警預(yù)報(bào)的重要依據(jù)，在滑坡實(shí)時(shí)監(jiān)測過程中起指導(dǎo)作用[1-2]?；谧冃伪O(jiān)測數(shù)據(jù)建立有效的預(yù)測模型，不僅能減少監(jiān)測工作量與成本，同時(shí)又能降低變形發(fā)育致災(zāi)發(fā)生的機(jī)率。

滑坡作為一種常見的地質(zhì)災(zāi)害，通常由多種不確定因素，如地質(zhì)條件、地貌、水文地質(zhì)和物理因素以及人類活動(dòng)等共同作用導(dǎo)致。因此，需結(jié)合多種監(jiān)測手段對(duì)滑坡變形狀態(tài)進(jìn)行定量分析。而滑坡變形預(yù)測一直是滑坡預(yù)警預(yù)報(bào)的重要研究方向和研究熱點(diǎn)[3-4]。由于滑坡變形演化受季節(jié)性因素(如庫水位調(diào)度、周期性強(qiáng)降雨)影響，滑坡累積位移曲線通常表現(xiàn)出很強(qiáng)的非線性動(dòng)力學(xué)特征，如階躍性變化等。常規(guī)方法對(duì)此類滑坡進(jìn)行預(yù)測預(yù)報(bào)時(shí)，極易將階躍性的變形特征誤認(rèn)為滑坡已進(jìn)入臨滑階段，造成誤判[5]。因此傳統(tǒng)的分析方法很難對(duì)滑坡的穩(wěn)定性做出準(zhǔn)確評(píng)價(jià)。在分析此類滑坡時(shí)，需充分考慮滑坡位移與其外部影響因素間的關(guān)系，只有在深入研究兩者的響應(yīng)關(guān)系的基礎(chǔ)上才能準(zhǔn)確地實(shí)現(xiàn)滑坡的預(yù)測預(yù)報(bào)[6]。

由于滑坡累積位移時(shí)序包含滑坡內(nèi)部影響因素(地質(zhì)結(jié)構(gòu)、地形地貌)加持的趨勢部分，以及外部影響因素(降雨、庫水位變化)波動(dòng)的周期部分，將位移時(shí)序分解為趨勢項(xiàng)和周期項(xiàng)并分別建模預(yù)測變得尤為重要。支持向量回歸機(jī)(Support Vector Regression, SVR)模型是基于統(tǒng)計(jì)學(xué)習(xí)理論所建立[7]，目前已廣泛用于此類滑坡的周期項(xiàng)位移預(yù)測。然而，其預(yù)測性能易受模型參數(shù)選擇的影響，過去通常采用多種參數(shù)尋優(yōu)法解決這一問題[8-11]，并未考慮影響因子之間的線性關(guān)系，造成樣本數(shù)據(jù)冗余；同時(shí)趨勢項(xiàng)位移預(yù)測多采用多項(xiàng)式擬合的方法，容易造成擬合結(jié)果不穩(wěn)定、收斂效果差。因此，本文通過采用核主成分分析(Kernel Principal Components Analysis, KPCA)實(shí)現(xiàn)影響因子的特征提取[12]，消除影響因子間的相關(guān)性，提高影響因子的質(zhì)量，降低樣本數(shù)據(jù)冗余；并通過二次指數(shù)平滑(second exponential smoothing, DES)擬合趨勢項(xiàng)位移[13]，以消除多項(xiàng)式擬合結(jié)果的不穩(wěn)定性。相對(duì)于網(wǎng)格搜索算法(Grid Search, GS)[14-15]、遺傳算法(Genetic Algorithm, GA)[16-17]，粒子群(Particle Swarm Optimization, PSO)尋優(yōu)算法[18-20]結(jié)構(gòu)較為簡單，且尋優(yōu)性能更勝一籌。因此，建立結(jié)合PSO算法確定SVR模型參數(shù)的最優(yōu)組合模型。以三峽庫區(qū)白水河滑坡為分析對(duì)象，在KPCA特征分析的基礎(chǔ)上，建立PSO-SVR聯(lián)合DES的組合預(yù)測模型實(shí)現(xiàn)滑坡位移預(yù)測分析。

1 模型原理

1.1 時(shí)間序列數(shù)據(jù)預(yù)處理

通常時(shí)間序列是由3部分組成(趨勢項(xiàng)、周期項(xiàng)和隨機(jī)項(xiàng))[21]。在滑坡位移的監(jiān)測數(shù)據(jù)中，趨勢項(xiàng)是受滑坡體內(nèi)部因素(地質(zhì)結(jié)構(gòu)、地形地貌)影響的線性部分，周期項(xiàng)是受滑坡區(qū)域環(huán)境的外部因素(降雨、庫水位波動(dòng))影響的非線性部分，而隨機(jī)項(xiàng)是隨機(jī)噪聲。因此，時(shí)間序列組成公式如下：

D=Dp+Dt+Dr

(1)

式中，D為總位移，Dp為周期項(xiàng)位移，Dt為趨勢項(xiàng)位移，Dr為隨機(jī)項(xiàng)位移。

本文采用小波去噪方法去除監(jiān)測數(shù)據(jù)中的隨機(jī)噪聲項(xiàng)。Hodrick-Prescott(HP)濾波器是Whittaker于1923年首次提出的[22]，這種方法被廣泛應(yīng)用于宏觀經(jīng)濟(jì)的趨勢分析研究。該方法采用對(duì)稱的數(shù)據(jù)移動(dòng)平均方法的原理，將時(shí)間序列中具有一定趨勢變化的平滑序列(趨勢項(xiàng))分離出來，然后原序列減去趨勢項(xiàng)得到周期項(xiàng)。得到趨勢項(xiàng)的公式如下：

(γt-1-γt-2)]2}

(2)

式中yt為時(shí)間序列，λ為平滑參數(shù)，γ為趨勢分量。在本研究中，滑坡位移每年急劇增加一次，因此，平滑參數(shù)λ的值被確定為100。

1.2 核主成分分析

為了更好地處理非線性數(shù)據(jù)，引入了非線性映射函數(shù)φ，將原空間的數(shù)據(jù)映射到高維空間中[12]，本文將滑坡影響因子序列xt映射到φ(xt)中，核主成分分析(kernel principal components analysis, KPCA)求解特征值如下：

(3)

(4)

(5)

1.3 二次指數(shù)平滑

指數(shù)平滑法是一種特殊的加權(quán)移動(dòng)平均法，其特點(diǎn)是最新數(shù)據(jù)的權(quán)重高于早期數(shù)據(jù)，此權(quán)重的因子隨著數(shù)據(jù)的老化依指數(shù)下降。二次指數(shù)平滑是指數(shù)平滑的衍生，它更適用于具有一定趨勢的時(shí)間序列的預(yù)測[13]。所以，本文采用DES擬合趨勢項(xiàng)位移。

對(duì)于一組時(shí)間序列{xi}，{si}代表第i時(shí)刻的指數(shù)平滑值，{bi}代表第i時(shí)刻的趨勢最佳估計(jì)。DES的預(yù)測值為Fi+m，代表第i次后的m次預(yù)測，通常m=1。DES算法的公式如下：

Fi+m=si+mbi

(6)

si=ζxi+(1-ζ)(si-1+bi-1)

(7)

bi=ξ(si-si-1)+(1-ξ)bi-1

(8)

式中，s1=x1，b1=x1-x0，ζ為數(shù)據(jù)平滑因子(0<ζ<1)，ξ為趨勢平滑因子(0<ξ<1)。

1.4 粒子群優(yōu)化

此算法在尋優(yōu)過程中，核心步驟是更新種群中每個(gè)粒子的位置和速度，而速度的更新最為關(guān)鍵。文獻(xiàn)[18-20]詳細(xì)說明了粒子群尋優(yōu)算法的尋優(yōu)原理。速度更新公式為：

vi=ω×vi+c1×rand()×(pbesti-xi)+

c2×rand()×(gbesti-xi)

(9)

式中，vi是粒子的當(dāng)前速度，w為慣性因子，rand()是隨機(jī)生成函數(shù)，生成值在(0-1) 之間，pbesti和gbesti分別為歷史最優(yōu)位和當(dāng)前最優(yōu)位，c1和c2為學(xué)習(xí)因子，通常c1=c2=2。位置更新的公式如下：

xi=xi+vi

(10)

1.5 支持向量回歸機(jī)

支持向量機(jī)(support vector machine，SVM)是一類按監(jiān)督學(xué)習(xí)方式對(duì)數(shù)據(jù)進(jìn)行二元分類的廣義線性分類器，最大邊距超平面是其決策邊界對(duì)學(xué)習(xí)樣本的求解[8-11]。它成功地解決了BP神經(jīng)網(wǎng)絡(luò)的過學(xué)習(xí)現(xiàn)象及局部最小問題，并在小樣本訓(xùn)練中體現(xiàn)出優(yōu)異的性能。因此，本文采用SVR模型對(duì)滑坡周期項(xiàng)位移進(jìn)行預(yù)測。

在非線性的變換之后，支持向量機(jī)使輸入的向量映射到高維特征空間中。通過構(gòu)造最佳的決策函數(shù)，將原有空間的核函數(shù)來替換高維特征空間中的點(diǎn)積運(yùn)算。應(yīng)用有限的樣本學(xué)習(xí)訓(xùn)練以獲得全局最優(yōu)值。SVR估計(jì)函數(shù)為:

f(x)=WTφ(x)+b

(11)

將估計(jì)函數(shù)通過ε不敏感損失函數(shù)轉(zhuǎn)化為優(yōu)化問題，SVR回歸預(yù)測模型可以通過二次規(guī)劃算法獲得：

(12)

2 滑坡變形特征分析

2.1 白水河滑坡GPS監(jiān)測點(diǎn)布設(shè)

根據(jù)白水河滑坡的地形、地質(zhì)條件、變形特征、觀測通視情況，確定該滑坡的監(jiān)測內(nèi)容主要為地表位移監(jiān)測、鉆孔測斜監(jiān)測、地下水位監(jiān)測等。監(jiān)測初期布設(shè)的7個(gè)GPS監(jiān)測點(diǎn)分布在3條縱向剖面上(圖1)，中間剖面有3個(gè)GPS監(jiān)測點(diǎn)(ZG118、ZG119、ZG120)，兩側(cè)剖面分別布設(shè)2個(gè)GPS監(jiān)測點(diǎn)(ZG91、ZG92、ZG93、ZG94)。2004年7月，由于白水河滑坡變形的突變產(chǎn)生較大的裂縫，經(jīng)專家評(píng)估建立了滑坡預(yù)警區(qū)，而后在預(yù)警區(qū)內(nèi)先后增建了4個(gè)GPS監(jiān)測點(diǎn)(XD-01、XD-02、XD-03、XD-04)，在滑體外圍東西兩側(cè)基巖脊上各建立1個(gè)GPS基準(zhǔn)點(diǎn)?；聟^(qū)內(nèi)共11個(gè)GPS變形監(jiān)測點(diǎn)，均采用美國產(chǎn)天寶GPS接收機(jī)(平面精度5+1ppm)進(jìn)行滑坡地表位移變形監(jiān)測。以上監(jiān)測自2003年6月上中旬三峽大壩壩前水位蓄水至135m開始獲得起始數(shù)據(jù)，以后每月監(jiān)測一次。其中，位于滑坡位移活躍區(qū)內(nèi)的監(jiān)測點(diǎn)有ZG93、ZG118和XD-01、XD-02、XD-03、XD-04。

圖1 白水河滑坡GPS監(jiān)測點(diǎn)Fig.1 Layout of GPS monitoring points of Baishuihe Landslide

2.2 滑坡體地質(zhì)工程特征

白水河滑坡位于三峽大壩以西約56km處，是三峽水庫重要的滑坡監(jiān)測區(qū)域。白水河滑坡為深層大型土質(zhì)滑坡，處于長江南岸單斜順層斜坡，該結(jié)構(gòu)單元位于白福坪背斜北翼，屬于侵蝕堆積谷的低山丘陵地貌?；碌钠矫嫘螒B(tài)為不規(guī)則矩形，左右寬度約450m，長度約352m，面積約16×104m2，體積約550萬m3，變形體的平均厚度約為35m，主滑方向?yàn)镹NE20°。如圖2所示，滑坡面為基巖與殘坡積層的接觸帶，厚度約為0.9～3.1m，呈折線形；滑坡體的基巖性為中厚巖砂夾薄層泥巖，產(chǎn)狀15°∠36°，巖層中節(jié)理裂隙發(fā)育[23]；滑坡體主要由第四系殘坡積碎石土組成，碎石含量20%～40%。

2.3 滑坡體變形特征

白水河滑坡屬于古滑坡堆積體，自2003年6月三峽水庫蓄水以來，滑坡開始復(fù)活，變形逐漸加劇。專業(yè)監(jiān)測數(shù)據(jù)顯示，白水河滑坡從2005年6月開始，發(fā)生了幾次嚴(yán)重的宏觀變形。其中，2007年出現(xiàn)了監(jiān)測以來最劇烈的一次宏觀變形；2008—2014年，變形逐漸趨于平緩，在強(qiáng)降雨的作用下，滑坡局部出現(xiàn)了淺層坍塌，但其宏觀變形跡象未出現(xiàn)明顯的變化。

圖2 滑坡工程地質(zhì)剖面Fig.2 Engineering geological section of landslide

圖3 白水河滑坡累積位移-庫水位-月降雨量-時(shí)間關(guān)系Fig.3 Cumulative displacement-reservoir water level-monthly rainfall-time diagram of Baishuihe Landslide

自2003年7月布設(shè)專業(yè)監(jiān)測以來，白水河滑坡監(jiān)測點(diǎn)的變形呈現(xiàn)階躍性特征(圖3)。結(jié)合圖1分析，預(yù)警區(qū)內(nèi)監(jiān)測點(diǎn)在每年5～9月，滑坡變形整體增大，10月到次年4月趨于平緩，說明庫水位漲落和降雨疊加作用對(duì)滑坡的穩(wěn)定性有較顯著的影響。每年5～9月份為三峽庫區(qū)的雨季，在此期間，該區(qū)降雨持續(xù)時(shí)間久、強(qiáng)度大，滑坡受季節(jié)性降雨的影響，表現(xiàn)出階躍性的變化特性，其中6～7月滑坡位移變化最為明顯，特別在2007、2008和2009年受強(qiáng)降雨的影響，出現(xiàn)非常明顯的位移劇增現(xiàn)象。每年5月開始，降雨急速增加，三峽水庫為防止洪災(zāi)開始放水，水庫水位明顯下降，因此水庫水位與降雨量呈負(fù)相關(guān)的關(guān)系。當(dāng)降雨量增加且水庫水位顯著下降時(shí)，往往伴隨著滑坡位移量的急劇變化，如2007年6月滑坡前緣受庫水作用強(qiáng)烈，監(jiān)測點(diǎn)XD-01、XD-03累積位移增長明顯高于中西部監(jiān)測點(diǎn)XD-02、XD-04、ZG93和ZG118，進(jìn)一步說明了滑坡體東北部結(jié)構(gòu)的不穩(wěn)定性。相比之下，非預(yù)警區(qū)監(jiān)測點(diǎn)多年變化較小且呈直線型。

通過上述分析可知，如果水庫水位下降期間伴隨強(qiáng)降雨作用，會(huì)對(duì)滑坡的穩(wěn)定性產(chǎn)生更不利的影響。因此，本文以活躍區(qū)內(nèi)典型位移監(jiān)測點(diǎn)ZG118為代表，結(jié)果同期的降雨量和庫水位高程對(duì)其滑坡位移時(shí)序進(jìn)行分析預(yù)測。

3 滑坡位移預(yù)測

3.1 位移預(yù)測技術(shù)路線

通過分析，本文以2003年7月至2013年2月的GPS滑坡位移監(jiān)測數(shù)據(jù)為分析對(duì)象，綜合考慮庫水位、降雨量對(duì)滑坡位移的影響，構(gòu)建滑坡位移組合優(yōu)化動(dòng)態(tài)預(yù)測模型，處理流程如圖4所示。

圖4 預(yù)測流程圖Fig.4 Forecast flow chart

(1)對(duì)位移監(jiān)測時(shí)間序列和初始影響因子進(jìn)行預(yù)處理，采用小波去噪方法去隨機(jī)噪聲，并進(jìn)行歸一化處理。

(2)利用HP濾波器將位移監(jiān)測時(shí)間序列分解為周期項(xiàng)位移和趨勢項(xiàng)位移。

(3)現(xiàn)有三峽庫區(qū)內(nèi)滑坡失穩(wěn)的持續(xù)形變都與庫水位漲落及降雨導(dǎo)致地下水滲流場改變有關(guān)，庫水位和降雨與位移變化的相關(guān)性強(qiáng)。楊背背[23]、Chao[24]等分析了變形影響因素的滯后性，擇取單期或多期庫水位和降雨數(shù)據(jù)作分析及模型輸入。因此，綜合考慮庫水位、降雨量的滯后影響，選擇每月庫水位、每月庫水位變化幅值、每月庫水位變化速率、每月降雨量和前兩月累積降雨量5個(gè)變量，構(gòu)建滑坡周期項(xiàng)位移的動(dòng)態(tài)響應(yīng)模型。

(4)對(duì)原始位移監(jiān)測數(shù)據(jù)與初始影響因素變量進(jìn)行核主成分分析，結(jié)果見表1。第一二主成分對(duì)原始數(shù)據(jù)的累積方差貢獻(xiàn)率達(dá)到82.37%，表明2個(gè)主成分包含了原始數(shù)據(jù)的大部分信息。因此選取這2個(gè)主成分作為新影響因子建立預(yù)測模型。

表 1 解釋的總方差

(5)趨勢項(xiàng)位移預(yù)測：選取2003年7月到2010年8月的86組白水河滑坡專業(yè)監(jiān)測點(diǎn)ZG118的觀測數(shù)據(jù)為訓(xùn)練樣本，并以2010年9月到2013年2月的30組數(shù)為測試樣本，利用MATLAB編程實(shí)現(xiàn)DES預(yù)測模型，數(shù)據(jù)平滑因子ζ為0.99，趨勢平滑因子ξ為0.98，模型預(yù)測結(jié)果如圖5所示。

(6)周期項(xiàng)位移預(yù)測：預(yù)測模型的樣本選取同上，使用MATLAB編程實(shí)現(xiàn)PSO參數(shù)尋優(yōu)，得到SVR徑向基核函數(shù)參數(shù)γ=0.65938，C=2.1656。根據(jù)最優(yōu)模型參數(shù)，學(xué)習(xí)訓(xùn)練樣本數(shù)據(jù)，構(gòu)建滑坡周期項(xiàng)位移與影響因子間的動(dòng)態(tài)響應(yīng)模型，同時(shí)對(duì)滑坡的周期項(xiàng)位移進(jìn)行預(yù)測，并通過測試樣本評(píng)估模型的預(yù)測精度和泛化性。

(7)將DES趨勢項(xiàng)預(yù)測結(jié)果和PSO-SVR周期項(xiàng)預(yù)測值通過加法模型，得到滑坡總累積位移結(jié)果，如圖7所示。

圖5 DES趨勢項(xiàng)預(yù)測結(jié)果Fig.5 Forecast results of trend items of DES

圖6 PSO-SVR周期項(xiàng)預(yù)測結(jié)果Fig.6 Forecast results of periodic items of PSO-SVR

圖7 滑坡累積位移預(yù)測結(jié)果Fig.7 Prediction results of the cumulative displacement of landslide

3.2 預(yù)測結(jié)果分析

為進(jìn)一步對(duì)比分析本文算法的優(yōu)越性，在同樣完成主成分分析處理后，使用3種不同的優(yōu)化算法來搜索SVR的模型參數(shù)，并建立相應(yīng)的響應(yīng)模型預(yù)測滑坡周期項(xiàng)位移。計(jì)算預(yù)測模型結(jié)果的均方根誤差(RMSE)、平均絕對(duì)百分比誤差(MAPE)、相關(guān)系數(shù)(R2)。

(13)

(14)

(15)

由表2可見，在相加了趨勢項(xiàng)位移的累積預(yù)測結(jié)果后，相較于基于網(wǎng)格搜索尋優(yōu)支持向量回歸機(jī)(GS-SVR)和遺傳算法尋優(yōu)支持向量回歸機(jī)(GA-SVR)的預(yù)測，粒子群尋優(yōu)支持向量回歸機(jī)(PSO-SVR)結(jié)果的RMSE和MAPE均為最小，而且R2最大，R2為0.89801，表明預(yù)測值與觀測值具有極強(qiáng)的相關(guān)性，即PSO-SVR模型在經(jīng)過KPCA分析后是最佳預(yù)測模型，并且該方法在分析影響因子對(duì)滑坡位移的變形影響作用中，具有極好的預(yù)測及泛化能力。從而證實(shí)了在參數(shù)優(yōu)化算法中，PSO算法對(duì)白水河滑坡預(yù)測的參數(shù)優(yōu)化明顯優(yōu)于傳統(tǒng)的網(wǎng)格搜索法和遺傳算法。同時(shí)，由于DES兼容了全期平均和移動(dòng)平均所長，通過定權(quán)的方式?jīng)Q定觀測數(shù)據(jù)的貢獻(xiàn)量，不舍棄過去的數(shù)據(jù)，使得其預(yù)測結(jié)果的收斂性明顯優(yōu)于多項(xiàng)式擬合。綜合以上分析可知，本文所提出的滑坡組合優(yōu)化預(yù)測模型具有很好的預(yù)測能力，滑坡位移預(yù)測效果可靠性高。

表 2 不同模型的累積位移預(yù)測誤差

4 結(jié)論

本文在綜合分析白水河滑坡GPS實(shí)時(shí)監(jiān)測位移及同時(shí)期降雨量、庫水位的監(jiān)測數(shù)據(jù)的基礎(chǔ)上，得出在降雨量和庫水位變化2個(gè)外部因子共同作用下影響滑坡的變形特征?？紤]到二者對(duì)滑坡作用具有滯后性，選取每月庫水位、每月庫水位變化幅值、每月庫水位變化速率、每月降雨量和前兩月累積降雨量5個(gè)變量作為分析的初始影響因子。據(jù)此，提出PSO-SVR聯(lián)合DES的滑坡位移預(yù)測優(yōu)化模型實(shí)現(xiàn)白水河滑坡的位移預(yù)測。首先，在對(duì)位移監(jiān)測時(shí)間序列進(jìn)行小波去噪處理的基礎(chǔ)上，利用HP濾波器將其分解為受地質(zhì)條件長期控制的趨勢項(xiàng)和受季節(jié)性因素影響的周期項(xiàng)，并通過核主成分分析將初始影響因子重構(gòu)為兩個(gè)主成分分量，用于SVR模型的構(gòu)建。模型構(gòu)建過程中，采用PSO算法尋優(yōu)SVR模型參數(shù)，學(xué)習(xí)訓(xùn)練樣本數(shù)據(jù)，對(duì)周期項(xiàng)位移進(jìn)行預(yù)測；并采用DES模型對(duì)趨勢項(xiàng)位移進(jìn)行預(yù)測。最后，將兩模型的預(yù)測值相加得到最終預(yù)測值，并通過與測試樣本序列對(duì)比分析進(jìn)行預(yù)測結(jié)果的精度評(píng)定。評(píng)定結(jié)果顯示，本文提出的組合優(yōu)化預(yù)測模型預(yù)測結(jié)果的RMSE為18.2031mm，R2為0.89801，明顯優(yōu)于基于網(wǎng)格搜索尋優(yōu)支持向量回歸機(jī)(GS-SVR)和遺傳算法尋優(yōu)支持向量回歸機(jī)(GA-SVR)模型預(yù)測結(jié)果，證實(shí)了該模型的預(yù)測能力和可靠性。