基于多重分形特征和分項組合預測聯合響應的滑坡預警預測研究

雷 恒 周曉嵐 王永強

1 小流域水利河南省高校工程技術研究中心，河南省開封市東京大道1號，475004 2 黃河水利職業技術學院水利工程學院，河南省開封市東京大道1號， 475004 3 長江科學院，武漢市黃浦大街23號，430010

滑坡是常見的地質災害之一，具有頻發性高、危險性大等特點[1-2]。為切實保證滑坡災害的安全防治，開展滑坡災害預警預測研究具有重要意義，已成為目前研究熱點。

在滑坡預警研究方面，董遠峰等[3]在滑坡變形特征分析基礎上開展監測預警研究，周辀等[4]結合數值模擬構建滑坡預警判據。上述研究雖涉及滑坡預警，但缺乏滑坡多重分形特征研究，且鄧小鵬[5]研究認為，基于滑坡變形監測成果的滑坡預警分級具有可行性。因此，仍可進一步拓展基于滑坡變形數據的多重分形特征研究及預警分級評價。

在滑坡預測研究方面，黃曉虎等[6]構建具有階躍型變形特征的滑坡預測模型；李秋全等[7]則在滑坡變形趨勢分析基礎上，利用優化支持向量機進行滑坡變形預測。上述研究在滑坡變形預測方面取得一定成果，但忽略了滑坡變形數據中誤差信息對預測精度的影響，且趙淑敏[8]研究認為，通過對誤差信息進行分離處理，能有效提高預測精度。因此，進一步開展滑坡變形的分項組合預測具有必要性。

基于上述研究成果認為，滑坡變形預警預測研究十分必要，且仍有進一步拓展的空間。因此，本文基于滑坡變形數據，首先利用多重分形消除趨勢波動分析(multifractal detrended fluctuation analysis, MF-DFA)和M-K分析進行滑坡變形的多重分形特征研究及預警分級評價；其次，在利用集成經驗模態分解法對滑坡變形數據信息進行分離處理基礎上，通過優化遞歸神經網絡和混沌理論實現滑坡變形的分項組合預測；最后，將多重分形特征研究結果和變形預測分析結果進行聯合響應，綜合評價滑坡變形規律，以便為其防治提供一定的理論指導。

1 基本原理

在前述分析基礎上，將論文分析思路進一步細分為：1)先利用MF-DFA模型開展滑坡變形數據的多重分形特征分析，在此基礎上構建預警判據及指標，以實現滑坡預警分級，旨在為其防治奠定理論基礎。2)先通過集成經驗模態分解法(ensemble empirical mode decomposition，EEMD)實現滑坡變形數據的分離處理，并通過優化遞歸神經網絡(recurrent neural networks，RNN)和混沌理論(chaos theory，CT)構建其分項組合預測模型，以實現滑坡變形預測研究。

1.1 多重分形特征分析模型構建

MF-DFA模型屬于多重非均勻分形方法，其不僅能揭示變形序列的多重分形特征，還能有效評價變形的發展趨勢，優越性較為明顯[9]。MF-DFA模型已被廣泛應用于巖土領域，因此利用其構建滑坡多重分形特征分析模型具有可行性。

結合MF-DFA模型基本原理，以滑坡變形數據為基礎，首先對其進行累積離差求解，并對累積離差序列進行子序列劃分，劃分依據為：子序列長度為s，個數為Ns=N/s(N為累積離差序列總數)，由于Ns存在非整數可能，難以發揮剩余累積離差節點的作用，因此提出對累積離差序列進行逆序重新劃分，則所得子序列個數為2Ns。其次，在q階波動函數階次條件下，求解其對應的波動函數F(q,s)為：

(1)

式中，F2(s,v)為第v個子序列的方差。

據前述可知，子序列長度s值與波動函數F(q,s)值具有一一對應關系，因此，通過改變子序列長度s可得到若干(s，F(q,s))，且兩者具有如下線性關系：

lnF(q,s)=C+h(q)·lns

(2)

式中，h(q)為q階波動函數條件下的Hurst指數，C為擬合常數。

將波動函數q值取值范圍設定為-8～8之間的偶數值，若h(q)值隨波動函數q值變化而變化，則說明滑坡變形序列具有多重分形特征，且通過h(2)值可判斷滑坡變形趨勢。具體判據為：1)當h(2)值屬于[0，0.5)區間時，滑坡變形具有反向持續性，變形會趨于減小，且h(2)值越小，趨勢性越強。2)當h(2)=0.5時，滑坡變形具有游離性，無法判斷其發展趨勢。3)當h(2)值屬于(0.5，1.0)區間時，滑坡變形具有正向持續性，變形會趨于增加，且h(2)值越大，趨勢性越強。

最后，由于波動函數q值與h(q)值具有對應關系，通過兩者求解奇異指數a(q)：

a(q)=h(q)+qh′(q)

(3)

式中，h′(q)為h(q)導函數。

以奇異指數a(q)為基礎，進一步計算多重分形譜寬度參數Δa、波形中大小波動所占比例參數Δf(a)：

Δa=amax-amin

(4)

Δf(a)=Δf(amax)-Δf(amin)

(5)

Δa參數主要用于評價滑坡變形序列的多重分形譜寬度，若Δa值越大，則多重分形強度也隨之越強，即波動越劇烈；Δf(a)參數主要用于評價滑坡變形序列波形中大小波動所占比例，Δf(a)值越小，大波動波形所占比例越大。依據文獻[10]研究成果，通過Δa和Δf(a)參數構建滑坡預警判據，以實現滑坡預警等級劃分，具體標準設定如表1所示。

表1 滑坡預警等級劃分標準

為實現表1中Δa和Δf(a)參數的趨勢判斷，再引入M-K分析方法。根據M-K分析方法的基本原理，將其評價參數Z值表示為：

(6)

式中，S為初步統計量，var(S)為初步統計量的方差值。

通過Z值大小即可判斷對應評價對象的趨勢特征：Z≥Za(Za值為對應顯著水平a條件下的臨界值，本文將顯著水平a設置為0.05，則Z0.05=1.960)說明評價對象具有增大趨勢，其越大表明趨勢性越強；-Za

綜上所述，通過M-K分析可判斷Δa和Δf(a)參數的趨勢，即可實現滑坡預警分級。

1.2 分項組合預測模型構建

滑坡變形影響因素較多，致使滑坡變形數據具有較大的波動性，且受不確定因素影響，滑坡變形數據會含有一定的誤差信息：

y(t)=q(t)+w(t)

(7)

式中，y(t)為滑坡變形值，q(t)為趨勢序列，w(t)為誤差序列。

由于存在誤差信息，在一定程度上會影響預測模型的訓練過程，即降低預測精度，這也體現出分形組合預測的優勢。經總結分析，將滑坡分項組合預測模型的預測流程分為3個步驟：1)先利用集成經驗模態分解將滑坡變形序列分離為趨勢序列和誤差序列。2)利用GOA-RNN模型實現趨勢序列的變形預測，以得到趨勢序列的預測值。3)將趨勢序列預測結果的預測誤差疊加至誤差序列中，組成新的殘差序列，并通過混沌理論實現其預測處理。將趨勢序列的預測值和殘差序列的預測值進行疊加，即得到滑坡變形的最終預測值。結合分項組合預測流程，將各階段的方法原理詳述如下。

首先，需對滑坡變形數據進行趨勢序列與誤差序列的分離處理。考慮到經驗模態分解可通過自適應正交基進行時頻信號處理，對非線性不平穩信號的處理能力較強，因此，以其為基礎構建滑坡變形數據的信息分離模型。需要指出的是，傳統經驗模態分解存在模態混疊現象，對分離效果具有一定影響，為克服該問題，集成經驗模態分解隨之產生，其通過增加白噪聲來解決頻率混疊問題。具體分析步驟為：1)將具有正態分布特征的白噪聲添加至滑坡原始序列中，使之重組成具有白噪聲信號的新序列。2)利用傳統經驗模態分解對新序列進行分離，以得到若干本征分量和一個剩余項。3)對本征分量和剩余項進行均值求解，得到EEMD最終分離結果：

(8)

式中，j為本征分量個數，IMF(t)為對應的本征分量，R(t)為剩余分量。

在前述信息分離基礎上，需進一步構建信息分離效果評價指標。以往研究多是利用單一指標進行評價，如均方根誤差、信噪比及平滑度指標等，由于不同指標的評價原理具有一定差異，單一評價指標已難以滿足應用需求，本文提出以3個基礎指標為基礎，構建綜合評價指標k：

k=g1+g2+g3

(9)

式中，g1、g2、g3為3個基礎指標對應的歸一化值。k值范圍為0～3，其值越大，說明信息分離效果越優；反之，分離效果越差。

其次，在信息分離基礎上，再構建趨勢序列的預測模型。考慮到RNN模型可一次性激活不同神經元，具有較強的即時處理能力，同時RNN模型在擬合過程中具有較快的收斂速度和更優的訓練精度，因此通過RNN模型構建趨勢序列的預測模型具有可行性[11]。

值得注意的是，RNN模型雖然具有顯著的優越性，但也存在一定不足，如連接權值及閾值具有較強的隨機性，會對預測精度造成一定影響。因此，為保證趨勢序列的預測精度，采用蝗蟲優化算法(grasshopper optimization algorithm，GOA)對其進行優化處理。結合GOA算法的基本原理，將其優化過程詳述為：1)對粒子位置進行初始化設置，如設定搜索空間上、下界，并設置最大迭代次數等。2)對每個粒子進行適應度計算，并保存最優粒子。3)改變粒子位置，并重新計算粒子適應度值，再對比前期保存的最優粒子，若前者更優，則對最優粒子進行替換；反之，繼續迭代尋優。4)當滿足最大迭代次數后，終止迭代，并輸出最優粒子對應的尋優參數，即完成RNN模型的連接權值及閾值優化處理。

GOA-RNN模型雖具有較優的預測精度，但由于滑坡變形的非線性特征，使得其仍會存在一定預測誤差，將其預測誤差疊加至誤差序列中，組成新的殘差序列。

最后，由于殘差序列具有較強的波動性和隨機性，混沌特性明顯，為實現高精度預測，提出利用混沌理論構建其預測模型。首先計算殘差序列的混沌指數λmax，λmax>0時表明所得殘差序列具有混沌特性，可通過混沌理論進行殘差序列的預測處理；然后通過延遲時間參數τ和嵌入維數M進行殘差序列的空間重構處理，并以相空間中的第i個相點ψi為預測中心，計算其與最近鄰點ψl之間的距離d：

(10)

當d值最小時，通過反推可得到殘差序列的預測值。

將趨勢序列GOA-RNN模型的預測結果和殘差序列CT模型的預測結果疊加，即為滑坡變形的最終預測結果。

2 實例分析

2.1 工程概況

白家包滑坡位于宜昌市秭歸縣，屬長江一級支流香溪河流域，距三峽工程41.2 km。據現場調查成果，滑坡縱向長度約550 m，前寬后窄，均寬約350 m，面積約19.25萬m2，平均厚度約45 m(厚度變化差異較大，中部較厚)，總體積約900萬m3，屬深層大型滑坡。滑坡區地形起伏較大，總體呈東低西高，其東側相對較平緩，坡度介于0°～25°；西側相對陡峻，坡度介于35°～50°。

自2003年三峽水庫蓄水以來，白家包滑坡歷年均會出現不同程度變形，且隨著時間持續，局部變形具有持續發展特征。經以往資料統計[12]，近年來白家包滑坡變形特征較為顯著。為實時掌握其穩定狀態，對其進行變形監測，其中主滑方向上布設ZG324和ZG325監測點，兩者監測成果較為完備。因此，本文以其2013～2018年監測成果作為后續分析的數據來源，監測頻率為1次/月，共計得到72期變形數據，并對兩個監測點的變形-時間曲線進行統計(圖1)。由圖可知，隨著時間持續，兩個監測點的變形持續增加，且具有階梯狀特征，這是由于，一方面滑坡月變形具有集中性特點；另一方面，從側面說明滑坡變形具有顯著非線性特征，對其變形進行分項組合預測具有必要性。

圖1 滑坡變形-時間序列Fig.1 Time series of landslide deformation

2.2 滑坡預警分級研究

利用MF-DFA模型對滑坡變形數據進行多重分形特征分析及預警分級研究。 首先，計算得到兩個監測點在不同階次q條件下的h(q)值(表2)。由表可見，隨著階次q減小，h(q)值也隨之減小，說明兩個監測點變形均具有多重分形特征；同時，兩個監測點h(2)值分別為0.714和0.694，均大于0.5，說明滑坡變形具有正向持續性，變形趨于增加，且ZG324監測點h(2)值相對更大，說明其較ZG325監測點趨勢性更強。

表2 多重分形特征分析結果

其次，利用M-K分析對Δa和Δf(a)參數進行趨勢判斷，以實現滑坡預警分級。結果分析如下。

1)Δa指標判據結果分析。通過計算統計，得到Δa指標判據結果(表3)。由表可見，兩個監測點的Z值不僅均大于0，還均大于Z0.05值，表明兩者均呈增大趨勢，且ZG324較ZG325具有相對更大的Z值，說明前者趨勢性更強。結合表1中判據標準可知，兩個監測點在Δa指標判據條件下的預警等級為Ⅱ級。

表3 Δa指標判據結果

2)Δf(a)指標判據結果分析。類比前述，再對Δf(a)指標判據結果進行統計(表4)。由表可見，ZG324監測點Z值為-2.304，具有減小趨勢；ZG325監測點Z值為-1.862，具有平穩趨勢，且其趨勢性小于ZG324監測點。結合表1中判據標準可知，兩個監測點在Δf(a)指標判據條件下的預警等級為Ⅱ～Ⅲ級。

表4 Δf(a)指標判據結果

3)最終預警結果分析。在Δa指標判據和Δf(a)指標判據結果基礎上，對兩個監測點的最終預警結果進行分析(表5)。由表可見，在兩個判據的預警結果中，僅ZG325監測點在Δf(a)指標判據條件下的預警等級為Ⅲ級，其余均為Ⅱ級。因此按照不利原則綜合得出，兩個監測點的最終預警等級均為Ⅱ級，即滑坡變形趨向不利方向發展，破壞風險一般，建議加強監測頻率及巡視，并做好防災預案。

表5 滑坡最終預警結果

2.3 滑坡變形預測分析

在滑坡預警分級研究基礎上，再利用分項組合預測模型進行滑坡變形預測分析，即將分析過程分為如下兩步。

2.3.1 滑坡變形數據分離處理

對兩個監測點優化前后的分離效果進行統計(表6)。由表可見，經優化處理，EEMD模型在兩個監測點中的綜合評價指標k值均有較大提高，且EMD模型的平均綜合評價指標k值為2.416，而EEMD模型為2.711，表明后者具有更明顯的優勢。通過前述分析可知，EEMD模型具有更好的分離效果，說明其優化處理過程具備有效性。

表6 EMD模型優化前后的分離結果

2.3.2 變形分項組合預測分析

在滑坡變形數據分離處理基礎上，為充分驗證分項組合預測模型的穩定性及滾動預測能力，加之考慮到監測成果的周期較長，因此，將預測過程劃分為兩個階段，即前期預測和后期預測，其中，前期預測樣本范圍為1～36周期，后期預測樣本范圍包含所有樣本(1～72周期)。同時，在預測過程中，將相應預測階段中的后5個樣本作為驗證樣本，并以ZG324監測點為例，詳述不同預測階段的預測效果。

1)前期預測結果分析

先利用GOA-RNN模型對ZG324監測點趨勢項序列進行預測處理，且為驗證GOA算法對RNN模型的預測效果，對優化前后的預測結果進行統計(表7)。由表可見，經GOA算法優化處理，5個驗證樣本的相對誤差值均出現不同程度減小，充分說明GOA算法對RNN模型的參數優化具備有效性；在GOA-RNN模型的預測結果中，相對誤差變化范圍為2.54%～2.91%，預測精度一般，從側面說明進行后續殘差序列預測的必要性。

表7 ZG324監測點前期趨勢項預測結果

由表7可知，趨勢項序列的預測結果存在一定預測誤差，將其疊加至誤差項序列中，得到ZG324監測點在前期的殘差序列，并通過計算得到λmax=0.057>0，說明殘差序列具有混沌特性，可利用混沌理論實現其預測處理。

通過CT模型對殘差序列進行預測處理，得到ZG324監測點在前期的最終預測結果(表8)。由表可見，經CT模型處理，ZG324監測點在前期預測結果中的相對誤差范圍為2.13%～2.27%，平均相對誤差為2.21%，相較其趨勢項預測結果，預測精度得到一定提高，表明CT模型具備有效性，也初步說明GOA-RNN-CT模型在滑坡變形預測中的適用性。

表8 ZG324監測點前期最終預測結果

類比ZG324監測點前期預測過程，對ZG325監測點進行前期預測，表9為預測結果。由表可見，在ZG325監測點的前期預測結果中，相對誤差范圍為1.99%～2.18%，平均相對誤差為2.12%，其預測效果略優于ZG324監測點。上述分析充分說明，GOA-RNN-CT模型在滑坡變形預測中的穩定性較強。

表9 ZG325監測點前期最終預測結果

2)后期預測結果分析

在前期預測基礎上，再通過后期預測來驗證組合預測模型的滾動預測能力，并實現外推預測，以評價滑坡變形的發展趨勢。經計算統計，兩個監測點的后期預測結果見表10。由表可見，兩個監測點后期預測效果相當，其中，ZG324監測點后期預測結果的相對誤差范圍為1.96%～2.13%，平均相對誤差為2.05%；ZG325監測點后期預測結果的相對誤差范圍為1.85%～2.04%，平均相對誤差為1.99%，均具有較高的預測精度，表明GOA-RNN-CT模型在滑坡預測中不僅具有較優的預測效果，還具有較強的滾動預測能力。同時通過外推預測可知，滑坡變形仍會進一步增加，與前述趨勢判斷結果一致。

表10 滑坡后期預測結果

最后，類比前述后期預測流程，再以支持向量機和GM(1,1)模型進行預測，以進一步佐證對比本文預測模型的有效性，所得結果如表11(單位%)所示。由表可見，在相應監測點的預測結果中，本文預測模型均具有更小的平均相對誤差值，充分說明本文預測模型相較傳統預測模型具有更優的預測效果。

表11 不同預測模型的預測結果

綜上分析可知，分項組合預測在滑坡變形預測中具有較優的預測效果及穩定性，且預測結果顯示，滑坡變形會繼續加大，并向不利方向發展。

將滑坡預警分級結果和變形預測結果進行聯合響應可知，滑坡變形現處于不利等級，并向不穩定方向持續發展，應采取切實有效措施，避免滑坡災害發生。

3 結 語

通過多重分形特征及分項組合預測在滑坡預警預測中的聯合響應分析，主要得出以下結論：

1)通過MF-DFA模型對滑坡變形的多重分形特征進行分析可知，滑坡變形具有多重分形特征，并具有正向持續性，變形會趨于增加。

2)通過變形預測分析得出，分項組合預測模型在滑坡變形預測中具有較優的預測效果及穩定性，且通過外推預測可知，滑坡變形會繼續增加，無收斂趨勢，趨向不利方向發展。

3)將多重分形特征研究結果和變形預測分析結果進行聯合響應綜合得出，滑坡現有預警等級相對不利，且后續變形還會進一步增加，趨向不穩定方向發展，建議對滑坡采取必要防治措施。