孫曉云，段 綽，王明明，鄭海青，靳 強

(1.石家莊鐵道大學電氣與電子工程學院，河北 石家莊 050043；2.河北金隅鼎新水泥有限公司，河北 石家莊 050200)

滑坡作為一種分布廣泛、發生頻繁且影響嚴重的地質災害，直接威脅著人們的生命安全和財產安全。邊坡位移預測作為一種避免滑坡地質災害的有效手段，已經成為國內外眾多學者的熱點課題[1-2]，具有重要的經濟價值和社會意義。目前，常用的邊坡位移預測模型可以分為兩種，一種是基于力學模型和有限元方法，模擬邊坡受多因素影響下的行為特征，進而采用數值方法預測邊坡變形[3]；另一種是以邊坡變形監測資料為基礎，結合數學模型進行分析預測，數學模型又可以分為時間序列模型[4]、神經網絡模型[5]和各種耦合模型[6-7]。研究表明，對位移序列采用信號分析的方法進行數據預處理可以提高預測性能，基于信號分析的耦合模型受到研究者的廣泛關注，周超等[8]提出了小波分解-極限學習機預測模型，但小波分析的分解效果取決于基函數的選取，自適應性差；韓永亮等[9]提出了局域均值分解模型(LMD)、蝙蝠算法優化的極限學習機模型，取得了不錯的預測效果，然而LMD模型的迭代次數少，端點效應較輕，判斷純調頻信號的條件需要試湊，影響了算法的精度；變分模態分解模型(VMD)相比于EMD模型、LMD模型等遞歸篩選模式[10]，可以更好地分解原始序列。張妍等[11]對比了EEMD-LSSVM模型和VMD-LSSVM模型兩種模型對風速預測結果，表明采用VMD模型進行分解可以獲得更高的預測精度。

鑒于上述研究，本文以河北省某水泥廠邊坡數據為例，首先采用VMD模型將邊坡位移序列分解為一系列相對平穩的分量，再利用鳥群算法優化核極限學習機的模型參數，通過對每個分量構建BSA-KELM模型進行滾動預測，最后疊加各分量預測結果，實現邊坡累計位移預測。

1 VMD和BSA-KELM預測模型

1.1 變分模態分解

變分模態分解(VMD)是2014年由DRAGOMIRETSKIY和ZOSSO在EMD模型的基礎上提出的，不同于EMD模型的遞歸求解方式，VMD模型是一種非遞歸、自適應的信號分解方法[12]。其基本原理是把原始時間序列信號分解為K個調頻調幅的子信號(IMF)，每個IMF是中心頻率不同的有限帶寬。子信號表示公式見式(1)。

uk(t)=Ak(t)cos(φk(t))

(1)

式中：Ak(t)為瞬時幅值；φk(t)為相位。

VMD模型的具體步驟可以分為構造變分問題和求解變分問題。在設置好模態數K、上升步長τ和懲罰參數α等參數后，變分問題可以表述為滿足K個IMF之和等于原始輸入信號f的前提下，使得尋求的各IMF估計帶寬之和最小，計算公式見式(2)。

(2)

式中：uk為分解得到的K個子信號；wk為各個子信號的中心頻率。

引入二次懲罰項α和拉格朗日算子λ將式(2)變成無約束問題(式(3))。

L({uk},{wk},λ)=

(3)

(4)

(5)

1.2 核極限學習機

ELM模型是HUANG等[13]根據廣義逆矩陣理論提出的單隱層前饋型網絡。由于ELM模型的輸入層到隱含層的連接權值和隱含層的偏置無需人為設定，所以其具有結構簡單、運行速度快、泛化性能好等優點，因此在分類和回歸問題上得到了廣泛的應用。當給定訓練樣本S={(xn,yn),n=1,2,…,N},其模型表示見式(6)。

(6)

式中：ω為輸入層和隱含層的權值；b為隱含層偏置；g(·)為激活函數；β為輸出權重。

通過最小二乘法，求解線性方程組Y=Hβ，并且為了提高ELM模型的泛化性能，引入了正則化系數C，輸出權重β表達式見式(7)。

(7)

式中：I為對角矩陣;Y為期望輸出。

在此基礎上，HUANG等又提出了KELM模型，以核映射替代ELM模型中的隨機映射。 定義核矩陣Ω=HHT,矩陣元素ΩELMi,j=h(xi)h(xj)=K(xi,xj),其中,K(·)為核函數。此時KELM模型的輸出函數可以表示為式(8)。

(8)

本文選取局部性強、泛化性好的徑向基核函數，得到的表達式見式(9)。

K(xi,xj)=exp(-γ‖xi,xj‖2)

(9)

1.3 相空間重構

邊坡位移含有非周期運動，具有混沌特性，研究混沌時間序列的基礎就是采用相空間重構，在相空間中進行混沌模型的建立和預測。給定滑坡位移序列{x1,x2,…,xN},當設定好嵌入維數m和延遲時間τ，便可將一維時間序列重構成多維的狀態空間，表達式見式(10)和式(11)。

Zt={y1,y2,…,yt}

(10)

(11)

式中，t=N-mτ+1。

選擇合適的嵌入維數和延遲時間是相空間重構的關鍵，邊坡位移中存在噪聲干擾和估計誤差，延遲時間不宜過大，并且由于邊坡位移數據采樣周期長、數據樣本少，一般選取1作為延遲時間。而對于嵌入維數一般常采用試算法、飽和嵌入維法和虛假鄰近點法，但是嵌入維數的選取與核極限學習機的拓撲結構有關，因此，本文采用鳥群算法確定嵌入維數的值。

1.4 鳥群算法

鳥群算法(BSA)是2015年基于鳥群行為提出的一種群智能優化算法。與粒子群優化算法、遺傳算法相比，BSA算法具有優化精度高、收斂速度快、魯棒性好等優點。BSA算法主要有鳥群的飛行、覓食和警戒三個群體行為[14]，其依賴飛行行為跳出局部最優進行全局搜索，通過覓食行為記錄個體和群體最后的解，通過覓食和警戒行為的隨機切換搜索當前局部的最優解，通過飛行間隔FQ來平衡算法全局搜索和局部搜索能力。鳥群算法流程圖如圖1所示。

圖1 鳥群算法流程Fig.1 Flow of bird swarm algorithm

2 評價指標和算法流程

2.1 評價指標

為了定量、準確地表示預測模型的精度，本文采用平均絕對誤差(MAE)、均方根誤差(RMSE)和擬合優度(R2)這三個指標對預測結果進行誤差分析，具體表達式見式(12)～式(14)。

(12)

(13)

(14)

2.2 建模流程

本文提出的VMD-BSA-KELM預測模型建模流程如圖2所示，具體步驟如下所述。

圖2 VMD-BSA-KELM模型位移預測流程圖Fig.2 Displacement prediction flow chart of VMD-BSA-KELM model

1) 使用VMD將原始的邊坡位移序列分解為K個子序列。

2) 對各分量分別建立BSA-KELM預測模型，采用BSA算法對相空間重構的嵌入維數m和KELM模型中的懲罰系數C和核參數γ三個數值進行聯合優化。

3) 將優化后確定的嵌入維數m代入相空間中對原始序列進行重構，將懲罰系數C和核參數γ帶入KELM模型中建立最優模型進行預測。

4) 將各個子序列的預測結果疊加得到最終的邊坡位移預測值。

5) 運用式(12)～式(14)進行誤差分析。

3 應用實例

3.1 監測點實例分析

本文選取了河北省某水泥廠作為邊坡位移監測地點，實際監測中劃分了5個監測斷面、13個監測點，監測點的布置如圖3所示。每個監測點安裝一套GNSS監測設備，點位數據通過本地局域網絡傳輸到數據處理中心進行實時存儲及數據處理，所有數據均可長期保存、不易丟失，數據可恢復性強，GNSS觀測墩如圖4所示。

選取了2019年6月—2020年5月實際測量的邊坡位移數據進行案例分析，采用了變形量較大的監測點G123作為研究對象。采集數據樣本共314個，采樣周期為1 d，采用VMD-BSA-KELM預測模型對邊坡位移進行提前一步的滾動預測。因采樣周期為1 d，提前一步的預測即可提前一天預測出第二天的位移值，多步預測可以得出后幾天的值，但是隨著步數的增加，累計誤差也會激增，因此本文只針對單步預測進行討論。選取前250個數據為訓練樣本，后64個數據為測試樣本，得到的邊坡位移曲線如圖5所示。

圖5 邊坡位移曲線Fig.5 Slope displacement curve

采用VMD對邊坡位移序列進行分解，在分解過程中，K的取值直接影響著信號序列的分解效果，二次懲罰系數α的值也直接影響模態分量的帶寬。本文采用觀察中心頻率的方法確定K的大小，取K值為5，具體規則可參見文獻[15]，α取VMD默認值2 000，分解后各分量的頻譜如圖6所示。

原始數據經過VMD模型分解之后分為5個模態分量，如圖6所示。 從圖6中可以看出，IMF1的平均振幅是5個模態分量中最大的，達到了4.259 mm，且變化趨勢較為平緩，變化周期長，可以看作是邊坡受自身內在因素影響的趨勢項位移；IMF2和IMF3的平均振幅分別為0.522 mm和0.572 mm，變化趨勢有規律性，變化周期短，可以看作是受降雨量、地下水位等外在影響的周期性位移；IMF4和IMF5的平均振幅分別為0.387 mm和0.394 mm，振幅小周期短， 隨機波動性強， 且含有噪聲分量較多，可以看作是受人為因素等影響的波動性位移。最后，對這5個模態分量分別建立BSA-KLEM預測模型，采用鳥群算法優化嵌入維數m和懲罰系數C和核參數γ三個數值。本文BSA初始參數見表1。

表1 BSA基本參數設置Table 1 BSA basic parameter settings

圖6 VMD分解結果Fig.6 VMD decomposition results

預測模型取各分量前250個采樣點進行訓練，得到最優模型，將建立的最優BSA-KELM模型對5個模態分量分別進行預測，預測值如圖7所示。從圖7中可以看出，各個模態分量的預測精度都較高，曲線也很好地擬合了原始曲線，各分量評價指標見表2。將5個分量的預測值進行疊加，即可得到邊坡位移的預測值。

圖7 各分量預測結果Fig.7 Prediction results of each component

表2 各分量評價指標Table 2 Evaluation indexes of each component

為了驗證不同模型的預測性能，建立了以下四種模型，并對結果與VMD-BSA-KELM預測模型進行對比分析。①KELM模型：直接建立KELM模型，對邊坡位移序列進行1步的滾動預測；②BSA-KELM模型：采用BSA算法優化嵌入維數m、懲罰系數C和核參數γ三個數值的BSA-KELM模型；③EEMD-BSA-KELM模型：采用集合經驗模態分解方法(EEMD)對位移序列進行分解，再使用BSA-KELM模型進行預測,EEMD分解結果如圖8所示。

圖8 EEMD分解結果Fig.8 EEMD decomposition results

所有的訓練和測試均在Matlab R2016a軟件下運行，采用Intel Core i7-9750處理器，16 G內存的計算機平臺。各種預測模型結果如圖9所示。

圖9 四種模型預測結果曲線Fig.9 Curve of prediction results of four models

從圖9中可以看出，直接對原始位移數據進行預測，預測值和實際值偏差較大，在各個突變點上均不能達到預測精度要求。當采用鳥群算法對預測模型進行優化，雖然精度有所提升，但是改善效果并不明顯，其原因可能是在KELM模型參數選取時，參數本身偏離最優參數不大，導致鳥群算法提升的空間不大。但是，當采用了VMD模型或者EEMD模型對信號進行分解后，雖然預測時間增長，但是預測精度均優于直接預測的模型，說明合適的信號預處理技術可以有效降低邊坡位移序列的突變特性，提高模型的預測性能和精度。

在兩種采用了信號分解的模型中，VMD-BSA-KELM模型預測精度明顯高于EEMD-BSA-KELM模型。主要原因是EEMD模型分解的高頻部分仍表現出復雜的非線性，使得預測精度不佳，而高頻部分所占幅值較大，導致整體預測精度下降，其根本原因還是基于遞歸過程的EEMD模型對原始信號分解不徹底。并且由于EEMD模型分解后的分量較多，各分量都要建立模型進行預測再求和，所以其運行時間也長于VMD-BSA-KELM模型。由于VMD模型對原始序列分解更為徹底，有效地降低高頻分量的幅值，使得預測模型在31點、44點、61點等突變樣本點的預測值更加接近實際位移值，有效地跟蹤位移序列的變化。

從表3中可以看出，VMD-BSA-KELM模型平均絕對誤差為0.165 9 mm，均方根誤差為0.212 0，擬合優度達到了0.988，各種評價指標均優于其他模型，說明VMD-BSA-KELM模型在一定程度上可以提高邊坡位移預測的準確性。

表3 四種預測模型評價指標Table 3 Evaluation indexes of four prediction models

3.2 其他監測點實例分析

為了驗證模型的泛化性，選取了三個不同測量點的邊坡位移數據分別建立了VMD-BSA-KELM模型，仍然使用前250個數據樣本作為訓練集，后64個數據樣本作為測試集，進行提前一步的滾動預測。 預測曲線如圖10所示，評價指標見表4。 三個測量點的平均絕對誤差均在0.20 mm以下，均方根誤差均在0.23 mm以下，擬合優度均達到0.98以上。

表4 不同測量點的預測結果Table 4 Prediction results of different measurement points

圖10 不同測量點的預測曲線Fig.10 Prediction curves at different measurement points

4 結 論

1)直接對邊坡位移序列建立預測模型，預測效果并不理想，精度較差。

2) BSA算法對KELM模型參數優化可以提高預測精度，但是提升空間有限。

3) 與EEMD模型相比，VMD模型更適合邊坡位移序列的分解，分解更為徹底，VMD-BSA-KELM模型比EEMD-BSA-KELM模型預測精度更高，且運行時間更短。