基于混沌時間序列的黃土滑坡變形預測方法及應用

王 利，岳 聰，舒 寶*，張耀輝，許 豪，義 琛

(1. 長安大學 地質工程與測繪學院，陜西 西安 710054； 2. 地理信息工程國家重點實驗室，陜西 西安 710054； 3. 長安大學 西部礦產資源與地質工程教育部重點實驗室，陜西 西安 710054； 4. 自然資源部第一大地測量隊，陜西 西安 710054)

0 引 言

滑坡本質上是一個開放、耗散和復雜的非線性動力學系統。在采用基于變形數據的滑坡預測模型(灰色系統模型[1-3]、指數平滑模型[4-6]等)進行變形預測時，由于沒有考慮所研究系統確定性與隨機性的關系，所以在一定程度上會影響預測結果的準確性[7]?；煦缋碚撌且环N非線性系統預測方法，主要研究自然界復雜系統的發展規律，可以直接根據序列本身的客觀規律進行預測，減少了人為因素的干擾，使得預測結果可靠性更高。盡管混沌系統自身對初始條件的敏感性導致了對其行為進行長期預測的困難，但是并沒有影響到對混沌系統發展規律的研究[8]。目前，混沌理論已廣泛應用于經濟學、氣候、電力等領域[9-16]。

已有研究表明，滑(邊)坡在演化過程中會通過一定途徑進入混沌狀態[17-21]?；诖?，本文基于陜西涇陽地區廟店滑坡GNSS變形監測數據序列與經過抑噪處理后的數據序列，在判定兩種變形監測數據序列為混沌時間序列的基礎上，引入幾種常見的混沌時間序列預測方法，并對這些方法的預測結果進行比較，以便選取適用于黃土滑坡變形預測的方法；同時，對這兩種數據序列分別進行預測結果的對比分析，判斷對原始監測序列進行抑噪處理后是否可以提高滑坡變形監測結果的預測精度。

1 理論模型

1.1 混沌相空間重構理論

混沌理論作為一種非線性系統理論，主要通過研究如何由時間序列通過相空間重構，從另一個維度和視角來辨識系統，挖掘系統中蘊藏的規律，并預測系統的未來走勢。在這個過程中，相空間重構的質量直接影響到系統的未來走勢。因此，需要采取有效且合適的方法重構相空間。Takens認為相空間重構最有效的方法為延遲坐標法[22]，其基本原理如下。

設時間序列為{x(i),i=0,1,2,…,N}，N是時間序列長度，以時間延遲(Δt)和嵌入維數(m)進行相空間重構，相空間中相點(Xi)的表達式為

Xi=

(1)

也可以表示為

Xi=[x(i),x(i+Δt),x(i+2Δt),…,

x(i+(m-1)Δt)]i=1,2,…M

(2)

M=N-(m-1)Δt

式中：M為相空間中相點的個數。

時間延遲和嵌入維數是相空間重構的兩個重要參數，這兩個參數的選取是否合理將直接決定相空間重構的質量。

相空間重構的延遲坐標法原理表明，由系統原始狀態變量構成的相空間和一維觀測值重構相空間里的動力學行為等價(拓撲學認為兩個空間是相同的)，即一維觀測值中包含系統所有狀態變量演化的全部信息。由此演化規律可得到系統下一時刻的狀態，從而得到時間序列下一時刻的預測值，這就為采用混沌時間序列預測提供了依據。

1.2 S-變換時頻濾波

變形序列的噪聲會直接影響變形預測結果，因此，本文采用S-變換時頻濾波方法削弱GNSS坐標序列的噪聲，其原理如下。

S-變換濾波過程的表達式為

h(t)=S-1(S(h(t))×HTF(τ,f))

(3)

(4)

式中：t為時間;f為頻率;τ為控制窗口移動的參數;h(t)為給定的時間信號;HTF(τ,f)為時頻濾波器;S和S-1分別為正變換與逆變換；D為指定信號的時頻通域。

由于變形監測過程中隨機噪聲分布在各個時間點處，所以D為整個時頻通域。

時間信號h(t)經過S-變換后得到二維復時頻矩陣S(t,f),將S(t,f)中各元素求模后的矩陣記為S|(t,f)|。由S-變換結果可知，時頻濾波器的設計主要有兩種方式：①根據S(t，f)設計濾波器；②根據S|(t，f)|設計濾波器。第一種方式會改變系數的相位，從而使S-逆變換重構得到的時域信號產生附加噪聲，因此，本文根據最大最小值理論，按照S|(t，f)|設計濾波器。定義如下時頻濾波器，其表達式為

(5)

式中：Mi為ti時刻振幅譜中S|(ti,f)|的最小值。

具體抑噪過程如圖1所示。

圖1 基于S-變換的抑噪過程

2 實驗分析

2.1 數據來源

本文所采用的觀測數據來自陜西省涇陽縣太平鎮廟店村，該村位于涇河南岸，是黃土滑坡頻發地帶?；麦w平面形態呈“上陡下緩”的形態[圖2(a)]。

圖2 陜西涇陽地區廟店滑坡現場及監測網

滑坡監測網[圖2(b)]中使用北斗/GNSS接收機進行觀測，并結合長安大學北斗分析中心自主研發的北斗/GNSS實時監測系統獲得的監測點N(北)、E(東)、U(高程)方向點坐標變化序列的小時解，數據采集時間為2018年8月28日17:00到2019年4月17日2:00，共1 913個歷元，基準站設在距離滑坡體約400 m處的穩定地帶。本文選取MD09監測點坐標變化序列進行分析。由于GNSS信號易受干擾，所以對MD09監測點GNSS變形監測結果采用S-變換時頻濾波方法進行抑噪處理[23]，結果如圖3所示。

圖3 原始數據序列和S-變換抑噪后的數據序列

2.2 滑坡變形時間序列的相空間參數重構

時間延遲和嵌入維數是相空間重構的兩個重要參數。本文以MD09監測點E方向GNSS變形監測序列和經S-變換抑噪處理后的GNSS變形監測結果序列為例進行相空間參數確定。

2.2.1 時間延遲的確定

在實際應用中，時間延遲的選取主要有平均位移法、自相關函數法以及互信息量法[24]。其中，自相關函數法主要利用相空間中各相點的線性相關性求解時間延遲。盡管該方法簡單易實現，但該方法不能推廣到高維相空間重構過程?；バ畔⒘糠▌t是在考慮了以上方法的局限性之后提出的一種解決時間延遲比較有效的方法。本文采用互信息量法求解MD09監測點E方向滑坡變形時間序列的時間延遲(圖4)。

圖4 MD09監測點E方向變形時間序列時間延遲t

從圖4可以看出，MD09監測點的滑坡變形時間序列互信息函數隨時間延遲變化明顯。隨著時間延遲的增大，互信息函數有明顯波動，這主要是因為利用互信息量法能夠很好地反映MD09監測點各時間序列的相關性以及整個時間序列的非線性性質。因此，該方法適用于確定滑坡變形時間序列的時間延遲。當時間延遲為2時，原始時間序列的互信息函數取得第一個極小值[圖4(a)]，因此，原始時間序列的時間延遲為2；當時間延遲為5時，原始時間序列的互信息函數取得第一個極小值[圖4(b)]，因此，原始時間序列的時間延遲為5。

2.2.2 嵌入維數的確定

常用確定嵌入維數的方法有虛假鄰近點法、改進的虛假鄰近點法(Cao算法)以及飽和關聯維法(G-P算法)。Cao算法是在虛假鄰近點法的基礎上針對其存在的不足改進而來，G-P算法和Cao算法在計算最佳嵌入維數時都容易實現，而且兩種方法所得結果可以相互論證[25]。因此，本文采用Cao算法獲取MD09監測點E方向變形時間序列的嵌入維數(圖5)。

圖5 MD09監測點E方向變形時間序列的嵌入維數m

從圖5可以看出：采用Cao算法計算所得嵌入維數的均值(E(m))受嵌入維數影響較小，其值在1附近有小幅度的波動，整體趨近于1。令E1(m)=E(m+1)/E(m)，E1(m)隨嵌入維數的變化而變換；當嵌入維數增加到某一個值后，E1(m)也趨于穩定。這說明監測點變形時間序列在計算之前所假定嵌入維數的范圍內，沒有因取值太小出現數據冗余，也沒有因取值太大遺漏部分數據或丟失數據間的相關關系。因此，計算結果較好，進而說明在涇陽地區廟店滑坡變形時間序列的相空間重構過程中，可以采用Cao算法確定最佳嵌入維數。

當嵌入維數增加到10時，E1(m)不再隨著維數的增加而改變[圖5(a)]，因此，原始時間序列E方向上的嵌入維數為10；當嵌入維數增加到8時，E1(m)不再隨著維數的增加而改變[圖5(b)]，因此，原始時間序列E方向上的嵌入維數為8。

2.3 黃土滑坡混沌時間序列預測性能分析

2.3.1 黃土滑坡混沌時間序列的判別

運用混沌理論研究時間序列的前提是該序列具有混沌特性，因此，需要進行MD09監測點GNSS變形序列混沌特性的識別。時間序列的混沌特性判定方法主要分為定性和定量兩個方面。定性判定方法簡單、直觀、方便，但該方法受人為主觀因素影響較大，導致對時間序列的混沌特性判別不準確；而定量判定方法主要是對通過計算滿足混沌特性的相關量值進行判定，避免了人為因素的影響。因此，本文采取定量判定方法中的Lyapunov指數進行時間序列的混沌特性識別。

Lyapunov指數的求解主要有wolf法和小數據量法，其中，小數據量法計算簡單且易于實現，適合處理有限長度的時間序列[26]。因此，本文采用小數據量法求解時間序列的Lyapunov指數。圖6為采用小數據量法分別計算的MD09監測點E方向原始時間序列和經過S-變換抑噪處理的時間序列最大Lyapunov指數。計算結果表明：原始時間序列最大Lyapunov指數為0.001 2，經過S-變換抑噪處理的時間序列最大Lyapunov指數為0.001 0，兩種時間序列的最大Lyapunov指數均大于0，為正數。因此，采用GNSS技術獲取的滑坡變形時間序列具有混沌特性，可以采用混沌理論分析GNSS滑坡監測數據。

圖6 MD09監測點E方向變形時間序列最大Lyapunov指數

2.3.2 GNSS混沌時間序列變形預測

本文采用的變形監測數據來自陜西涇陽地區廟店滑坡，監測時間序列共計1 913歷元(即每小時采一次樣)，其中前1 883歷元為訓練數據，后30歷元為預測數據。對MD09監測點的變形監測數據分別采用加權一階局域預測方法、最大Lyapunov指數預測方法以及BP神經網絡預測方法進行預測，預測方式選擇多步預測。為了更好地分析監測數據經不同方法預測后的精度，采用平均絕對誤差(MAE,δMAE)和平均相對誤差(MRE,δMRE)兩個評價指標來評價預測方法的精度。其表達式為

(6)

(7)

對MD09監測點E方向混沌時間序列進行預測，根據已經求出的時間延遲和嵌入維數對監測點原始序列進行相空間重構，并通過不同預測方法對序列后30歷元數據進行預測，預測結果如圖7所示。由圖7可以看出：采用加權一階局域預測方法和最大Lyapunov指數預測方法的預測結果變化趨勢與原始時間序列的變化趨勢大致相同，但預測值與原始時間序列實際值相差較大，最大偏差分別為2.4 mm和1.5 mm。相比較而言，BP神經網絡預測方法的預測值較為接近原始時間序列實際值，最大偏差為1.1 mm。為了定量比較3種預測方法的預測精度，本文分別計算了不同預測方法的絕對誤差、相對誤差兩個評價指標序列(圖8、9)。從3種預測方法對應的評價指標序列可以看出，BP神經網絡預測方法優于加權一階局域預測方法與最大Lyapunov指數預測方法，尤其是對前10歷元的預測，其預測結果的相對誤差在15%以內，絕對誤差在0.5 mm以內。整體而言，原始時間序列的預測值與實際值偏差較大。

圖7 MD09監測點E方向變形時間序列預測結果

圖8 MD09監測點E方向變形時間序列各模型相對誤差

2.3.3 基于S-變換的GNSS混沌時間序列變形預測

為了分析經S-變換抑噪處理后的預測效果，根據已經求出的時間延遲和嵌入維數對監測點原始序列進行相空間重構，并建立不同預測方法對序列后30歷元數據進行預測，預測結果如圖10所示。從圖10可以看出：加權一階局域預測方法與最大Lyapunov指數預測方法預測值與實際值偏差較大，最大偏差分別為0.8 mm和1.3 mm。而BP神經網絡預測方法的預測值比較接近實際值，最大偏差為0.3 mm，且與實際值序列具有相同的波動趨勢。為了定量比較3種預測方法的預測精度，本文分別計算了不同預測方法的絕對誤差、相對誤差兩個評價指標序列(圖11、12)。從3種預測方法對應的評價指標序列可以看出，BP神經網絡預測方法優于加權一階局域預測方法與最大Lyapunov指數預測方法，尤其是對前10歷元的預測，其預測結果的相對誤差在5%以內，絕對誤差在0.3 mm以內，基本能反映原始時間序列的變化特征。在前10歷元預測精度高，主要是因為混沌的不穩定性以及初值敏感性，導致混沌時間序列只能是短期預測。相比較而言，在S-變換抑噪處理之前，GNSS混沌時間序列的預測值與實際值偏差較大，其主要是在對GNSS原始時間序列進行相空間重構時，由于受到觀測噪聲的影響，混沌時間序列的奇異吸引子不能完全恢復，所以預測結果偏差較大。

圖9 MD09監測點E方向變形時間序列各模型絕對誤差

圖10 S-變換抑噪后時間序列預測結果

圖11 S-變換抑噪后時間序列各模型相對誤差

通過對比圖11與圖12可知，相較于原始數據預測結果，經過S-變換抑噪處理后的預測結果更加接近實際值，說明抑噪處理能夠提高數據的預測精度。為了進一步說明抑噪處理后的數據有較高的預測精度，本文通過平均絕對誤差和平均相對誤差兩項指標對預測結果進行分析，結果如表1所示。

圖12 S-變換抑噪后時間序列各模型絕對誤差

通過表1可以看到，經過S-變換抑噪處理后的時間序列預測結果的平均絕對誤差與平均相對誤差均小于原始時間序列預測結果的評定指標。其中，加權一階局域預測方法、最大Lyapunov預測方法和BP神經網絡預測方法經S-變換抑噪處理后的時間序列的平均絕對誤差與原始時間序列相比分別提高了66.7%、57.1%和75.0%，平均相對誤差分別提高了58.5%、56.5%和65.5%。這說明原始數據存在的噪聲對預測結果有較大的影響，且經過抑噪處理后預測效果有較明顯的提高。不管是原始時間序列，還是經過S-變換抑噪處理后的時間序列，采用3種預測方法后，BP神經網絡預測方法的平均絕對誤差與平均相對誤差均小于其他兩種預測方法預測結果的評定指標。經過原始時間序列預測的平均絕對誤差和平均相對誤差分別為0.4 mm和11.9%；經過S-變換抑噪處理后的時間序列預測的平均絕對誤差和平均相對誤差分別為0.1 mm和4.1%。

表1 不同時間序列各模型預測精度統計結果

3 結 語

由于滑(邊)坡變形普遍具有非線性性質，且在演化過程中會通過一定途徑進入混沌狀態。如果不考慮滑(邊)坡系統確定性與隨機性之間的關系，在一定程度上會影響預測結果的準確性。為此，本文結合陜西涇陽地區廟店滑坡GNSS變形監測數據，采用基于混沌理論的滑坡變形時間序列預測方法對該滑坡的變形時間序列進行了分析和預測。

(1)采用相同方法對GNSS滑坡變形監測的原始時間序列和經過S-變換抑噪處理后的時間序列求解各自的相空間重構參數，并進行混沌識別。結果表明：GNSS滑坡變形監測的原始時間序列與經過S-變換抑噪處理后的時間序列均具有混沌特性。

(2)采用不同的混沌時間序列預測方法對廟店滑坡GNSS變形監測的原始時間序列與經過S-變換抑噪處理后的時間序列分別進行預測。結果表明：經過S-變換抑噪處理后的時間序列預測結果更接近實際值，平均絕對誤差、平均相對誤差均優于原始時間序列的指標評定值；多種預測方法中，BP神經網絡預測方法的預測精度較好；GNSS變形監測原始時間序列預測的平均絕對誤差和平均相對誤差分別為0.4 mm和11.9%，經過S-變換抑噪處理后的時間序列預測的平均絕對誤差和平均相對誤差分別為0.1 mm和4.1%。綜上所述，GNSS變形監測原始時間序列中存在的噪聲對變形預測結果有較大影響；經過S-變換抑噪處理后，預測效果得到明顯改善。