基于主控因子分析與GM-IAGA-WNN聯合模型的平推式滑坡位移預測研究*
——以垮梁子滑坡為例

黃 健 李 橋 巨能攀 許 強 王昌明

(①地質災害防治與地質環境保護國家重點實驗室(成都理工大學) 成都 610059)(②成都理工大學 環境與土木工程學院 成都 610059)

0 引 言

平推式滑坡是廣泛分布在我國西南部三峽庫區、川東和川北紅層地區的一種特殊結構形式的滑坡(張倬元等， 2016)。此類滑坡多發育在近水平砂、泥巖互層的巖體中，巖層傾角一般僅3～10°、斜坡坡面一般也小于20°，如四川德陽中江縣垮梁子滑坡(1949，1981)、重慶巴南麻柳嘴滑坡(1998)、四川達州宣漢縣天臺鄉滑坡(2004)、青寧鄉巖門村滑坡(2007)及四川巴中南江縣黑山坡滑坡(2012)等。可見，該類滑坡普遍存在，但是針對其形成機制與主控因子的研究尚存在不少爭議，對其如何進行預警預報，也成為眾多國內學者重點關注的焦點之一(Xu et al.，2016)。

滑坡變形曲線是一條易受內外界因素影響且具有波動性的非平穩、不光滑曲線，籠統性對滑坡歷史變形數據進行單因素或單模型分析，無法達到理想預測效果。因此，現階段的主要研究方向是基于監測數據，提取主控因子建立主控因子-滑坡位移的非線性多因子綜合預測模型(李驊錦等， 2017)，相關學者已在該方面獲得諸多成果。吳益平等(2007)引入層次分析的思想，將滑坡累計位移劃分為趨勢項位移和具有不確定性的隨機項位移，提出了灰色-神經網絡聯合模型，取得了較好的預測結果。張俊等(2015)運用滑坡演化規律，分析滑坡主控因子，并基于時間序列將累計位移分解，分別利用不同方法構建預測模型，提高了預測精度(徐峰等， 2011； Du et al.，2013； 李驊錦等， 2016)。但是上述成果大多是基于從外部因子去分析滑坡變形規律，未考慮對滑坡變形更為敏感的坡體內部因子，“平推式滑坡”基于此進行預測結果可靠性有待提高。因此，如何將平推式滑坡各項因子合理整合并納入滑坡位移預測模型是研究的重難點。

鑒于此，在綜合分析各模型優劣的基礎上，本文以垮梁子滑坡為例，基于現場多年的監測數據，將滑坡累計位移分成趨勢項位移和變化項位移兩部分。首先利用GM(1， 1)灰色模型對滑坡趨勢項位移進行提取； 再利用灰色關聯度分析和相關性分析方法對坡體內外各項因子進行篩選，將優選的主控因子導入基于改進型自適應遺傳算法優化的小波神經網絡模型對滑坡變化項位移進行預測； 最后，疊加趨勢項位移和變化項位移得到滑坡累計位移預測值； 通過與傳統BP神經網絡模型、小波神經網絡模型和未優化遺傳算法-小波神經網絡聯合模型預測值的對比分析，驗證了本文所構建模型的合理性和優越性。

1 GM-IAGA-WNN聯合模型

滑坡在外部誘發因素(降雨、溫度和人類活動等)、坡體內部變化因素(滲壓、地下水位和土體含水率等)和自身地質條件(地形地貌、地質構造和地層巖性)的共同影響下，變形位移呈現出不確定性和隨機性，且不易預測。影響因素對滑坡變形的影響模式各不相同，自身地質條件決定了滑坡變形的總體趨勢； 坡體內部和外部的影響因素導致了滑坡位移的不規則變化(唐天國等， 2005)。因此，可以將滑坡累計位移分解成趨勢項和變化項兩部分：趨勢項位移是隨時間變化的近似單調增長的曲線，變化項位移為一個較為復雜的非線性時間序列。

1.1 GM(1， 1)灰色模型

本文運用GM(1， 1)灰色模型對滑坡位移的趨勢項進行提取。設原始滑坡累計位移時間序列為X(0)={x(0)(i)，i=1， 2，…，n}，n為數據個數。將x(0)進行一次累加(Aucumulated Generating Operator， 1-AGO)以便弱化隨機序列的波動性和隨機性，得到新數據序列：X(1)={x(1)(i)，i=1， 2，…，n}，其中

(1)

對x(1)(t)建立下述一階線性微分方程：

(2)

1.2 小波神經網絡模型

小波神經網絡(Wavelet Neural Network， WNN)是一種基于小波變換理論和人工神經網絡所構建的新型神經網絡模型。其基本思想是利用非線性小波基函數取代BP神經網絡中通常采用的非線性函數作為隱含層節點的傳遞函數(Szu et al.，1992； Zhang et al.，1992)，因此小波神經網絡較BP神經網絡具有更好的學習泛化能力，在滑坡變化項位移的預測中預測精度更高。本文選擇3層函數型小波神經網絡，模型結構如圖 1所示(Zhang et al.，2011)。

圖 1 小波神經網絡模型結構Fig. 1 Architecture of the wavelet neural network

設系統輸入數據集為X={xi，i=1， 2，…，n}，n為輸入數據個數，輸出數據集為Y={yi，i=1， 2，…，n}，n為輸出數據個數，則第i線路的輸出計算公式為：

(3)

式中，K為隱含層的神經元個數;N為輸出層的神經元個數;w1k，n、w2i，k為各層網絡的權值；b1k、b2i為各層網絡的閾值；ψa， b為隱含層激活函數為Morlet母小波基函數)，即：

(4)

(5)

式中，ak、bk為各層網絡的伸縮平移系數。

σ(t)為輸出層激活函數(Sigmoid函數)，即：

(6)

1.3 IAGA-WNN模型

對于WNN的優化問題可以表述為：基于樣本數據集P={xi，yi，i=1， 2，…，t}，t為數據個數，通過某種手段尋找最佳網絡參數，使得誤差函數Error最小且滿足精度要求(張松華等，2011)：

(7)

遺傳算法是通過模擬自然界進化過程的一種搜索最優解的方法，具有不受空間限制和能同時搜索多個優良點，并在附近繼續優化繁殖的特點，在復雜的工程優化問題中得到了廣泛引用，但也存在算法收斂速度慢、計算穩定性差和效率底下的問題。為進一步提高模型預測準確度和穩定性，避免參數的盲目試算，本文應用一種改進型自適應遺傳算法(Improved Adaptive Genetic Algorithm， IAGA)，考慮利用IAGA優化WNN的參數，參數主要為：伸縮系數、平移系數和各層間的權值和閾值。IAGA-WNN神經網絡運算的步驟(Khoa et al.，2004)：

(1)隨機初始化。對染色體進行編碼，本文采用實數編碼法。3層IAGA-WNN神經網絡各個體編碼構成如下：

akbkw1k,nw2i,kb1kb2i

(2)確定適應度值f。計算公式(Van et al.，1992)為：

(8)

(3)算法運算。由父代群體適應度利用遺傳算法的3種算法(選擇、交叉和變異)得到子代群體。選擇算法基于適應度比例的方法進行操作，各個體被選擇概率與其適應度大小成正比。交叉算法中依據一定的交叉概率Pc隨機地從被復制群體中選擇兩個個體進行交叉，Pc自適應調整公式(李延梅， 2012)為：

(9)

最后以變異概率Pm對每一個染色體內的各個基因(個體)進行變異算法計算，Pm自適應調整公式(李延梅， 2012)為：

(10)

式中，Pm1>Pm2>Pm3，取(0，1)之間的值。

(4)循環操作。將新一代群體插入原始種群中，并計算綜合體的適應度。

(5)誤差判定。計算模型誤差值，若達到預設值，則終止循環，得到最佳編碼個體，否則回轉步驟(3)。

(6)參數轉化與導入。將最佳編碼個體轉化成相應的網絡參數，并將其賦值給作為初始值進行仿真預測。

2 垮梁子滑坡變形預測應用

2.1 滑坡工程地質概況及監測設備布置情況

根據坡表變形強烈程度差異，將滑坡分為如圖 2的3個區域，坡南半側區域為滑坡Ⅰ區，滑坡北半側區域為滑坡Ⅱ區，孤立于北部山梁上發育的小滑坡為滑坡Ⅲ區。受地形地貌控制，滑坡區域地下水滲流場為獨立水文地質單元，主要補給來源為大氣降水，并且坡表的洼地積水和坡體沖積水也會對地下水進行入滲補給，沿風化裂隙、構造裂隙、層面裂隙從高處向低處運移，以泉水和滲流的形式向坡體前緣溝谷內排泄。

垮梁子滑坡于2011年初開始在滑坡Ⅰ區布設滑坡監測裝置，并于2013年6月先后完成了雨量監測裝置2個，地表位移監測裝置2個，滲透壓監測裝置8個，土壤含水率監測裝置9個，地下水位監測裝置7個。考慮數據可靠性和完整性，選取2013年6月至2017年3月作為研究時間區間，相關監測點的監測平面和剖面布置圖如圖 2和圖 3所示。

圖 2 垮梁子滑坡監測平面布置圖(單位：m)Fig. 2 Floor plan of monitoring arrangement in Kualiangzi landslide(unit： m)

圖 3 垮梁子滑坡監測剖面(1-1′)布置圖(單位：m)Fig. 3 Section plan of monitoring arrangement in Kualiangzi landslide(unit： m)

2.2 基于GM(1， 1)灰色模型的趨勢項位移預測

通過GPS01和GPS02累計位移曲線(圖 4)，可見滑坡位移隨時間增大而呈上升趨勢，并在局部由于內外界因素的影響呈現出不規律波動和階躍現象。為研究其規律性，利用GM(1， 1)灰色模型從滑坡累計位移中提取出不含隨機性變化的趨勢項位移，以便挖掘其潛在規律。

考慮到GM(1， 1)灰色模型短期預測精度較好，而長期預測由于不可期的隨機擾動效果一般，為此本文采用等維灰數遞補動態預測模型，在處理數據時進行“新陳代謝”，其過程為：讀入數據時設立變化周期t，到達周期數時在原始數據集中加入最新預測數據，并去掉周期內最舊的數據，形成等維的新數據集，進行下一步計算，如此循環往復，直至完成所有預測。

為避免原始離散序列的光滑度影響預測精度，通過觀察兩個位移監測點的曲線變化，將其分成2個部分： 2012年9月至2013年9月， 2013年10月至2017年3月(考慮監測數據總體完整程度，本文重點研究2013年6月以后部分)，并設置t=5，利用Matlab R 2016年軟件編程分段計算，得到GPS01和GPS02基于GM(1， 1)灰色模型的趨勢項位移與實際觀測值的對比圖(圖 4)：

圖 4 趨勢項位移提取值Fig. 4 Extracted values of periodic displacementsa. 監測點GPS01；b. 監測點GPS02

2.3 基于IAGA-WNN模型的變化項位移預測

滑坡變化項位移，即為滑坡累計位移實際監測值與上節提取的趨勢項位移的差值。

由于平推式滑坡的演化機理尚未得到統一的解釋(范宣梅等， 2008； 郭曉光等， 2013)，本文以多種監測數據為基礎，展開主控因子篩選分析：初步劃分為外部主控因子、內部主控因子和坡體自身穩定性因子，盡可能全面地挖掘出各因子與滑坡位移的關系，實現平推式滑坡變形趨勢預測，也在一定程度上探索其演化規律。同時為了避免所建模型參數繁雜、信息冗余，每項只選擇出最優因子參與計算。

如圖 5所示，在每年的雨季階段，都有向上躍遷的趨勢，因此，降雨是影響坡體變形的關鍵外部因子。同降雨相關的滲壓、地下水位及含水率作為影響坡體變形的內部因子進行分析。滑坡變形的大小作為坡體自身穩定性因子參與綜合預測模型的構建。各因子的具體量化關系主要是通過利用灰色關聯度分析和相關性分析(包含Person相關系數和t檢驗值)進行計算篩選，獲得最終影響坡體變形的關鍵因子序列。

圖 5 變化項位移提取值Fig. 5 Extracted values of Variable displacements

灰色關聯度和相關性分析應用十分廣泛，其在滑坡敏感性分析中也有實現(王貝莉， 2013； 鄧冬梅等， 2017)。本文給定灰色關聯度計算中分辨系數取ρ=0.5，灰色關聯度系數CG取值為0～1之間，且當CG越接近1則表明該因素對滑坡位移的影響越敏感，反之亦然。在計算Person相關性系數時，CP絕對值小于0.4時，可認定該因素與滑坡位移弱相關；CP絕對值介于0.4～0.6之間，可認定為中等相關；CP絕對值介于0.6～0.8時，可認定為強相關； 當CP絕對值大于0.8時，可認定為極強相關。在進行t檢驗時設置顯著性水平為?=0.05，即當t檢驗概率值Ct小于?時，則拒絕原假設，認為該因素與滑坡位移總體是相關的，且Ct越小說明相關性越顯著，并約定：當Ct小于0.01時，認定該因素與滑坡位移總體具有顯著相關性； 當Ct介于0.01～0.05之間，認定該因素與滑坡位移具有一般相關性。

在進行主控因子分析時，先基于灰色關聯度系數大小進行排列，再判斷Person相關系數，最后基于t檢驗值確定相關系數是否滿足顯著性檢驗，以杜絕虛假相關現象。

2.3.1 降雨主控因子

降雨事件對滑坡的影響方式主要分為持續性降雨和突發性降雨兩種形式。為了確定何種降雨形式對滑坡的影響程度更大，并避免相似影響因子的冗余，本文提出4種降雨類因子作為備選： ①單月降雨量：當月降雨總量； ②連續2月降雨量： 2個相鄰月降雨總量； ③截至年內累計降雨量：從年初至當月降雨總量； ④當月最大日降雨量：當月內記錄的最大日降雨量。

為分析降雨影響因子與滑坡位移之間的關系，以實際監測數據處理得到滑坡月位移、單月降雨量、連續2月降雨量、截止年內累計降雨量和當月最大日降雨量隨時間變化曲線，如圖 6所示。

圖 6 滑坡單月位移及多型降雨量監測曲線Fig. 6 Monitored monthly displacement of landslide and multiple types of rainfall

由圖 6可見，降雨主要集中在每年的雨季(6～10月)，在雨季期間內滑坡位移速率明顯增大，各降雨備選因子曲線的波動與滑坡位移速率波動曲線有較好的一致性。

基于灰色關聯度系數CG、Person相關性系數CP和t檢驗值Ct進行綜合篩選，最終計算結果如表 1所示：

表 1 降雨影響因子相關性分析結果統計Table 1 Correlation analysis of rainfall influencing factors

結合灰色關聯度分析和相關性分析，確定GPS01和GPS02的降雨影響因子為連續2月降雨量。

2.3.2 坡體內部主控因子

通過在坡體內部多種類型的監測裝置而得到的監測曲線可以發現，其變化與滑坡變化項位移有一定關聯性，但關聯程度各不相同，結合數據的完整性提出備選因子： ①滲透壓：坡體后緣滲透壓(SY01)、坡體中部滲透壓(SY03、SY04)、坡體前緣滲透壓(SY05、SY08)； ②地下水位：坡體中部地下水位(DX05)、坡體前緣地下水位(DX06、DX07)； ③土體含水率：坡體中部土體含水率(HH02)、坡體前緣土體含水率(HH03、HH05)。

基于原始監測數據得到滲透壓-時間變化曲線(圖 7)，土體含水率-時間變化曲線(圖 8)和地下水位-時間變化曲線(圖 9)：

圖 7 滑坡單月位移及多個滲透壓監測曲線Fig. 7 Monitored monthly displacement of landslide and multiple types of osmotic pressure

圖 8 滑坡單月位移及多個地下水位監測曲線Fig. 8 Monitored monthly displacement of landslide and multiple types of groundwater table

圖 9 滑坡單月位移及多個土體含水率監測曲線Fig. 9 Monitored monthly displacement of landslide and multiple types of soil moisture contenta. HH02； b. HH03 & HH05

從圖 7可以得出SY01監測曲線與滑坡月位移量曲線的變化規律較為一致，但具有一定的“滯后效應”，如2013年7月滑坡月位移達到最高值274.29imm和289.36imm，而滲透壓值則在1個月之后2013年8月達到最高值37.69ikPa和86.23ikPa。SY03、SY04和SY05位于坡體中部，該部分滲透壓一直處于平穩狀態與時間變化關系不大。SY08呈現出明顯的季節性變化，每年的6月份隨著滑坡區域雨季的來臨呈現上揚趨勢，而后在9月份到達峰值，再緩慢回落。除SY01外其他監測點與月位移量曲線協同性不強。

GPS01和GPS02的變化項位移與滲透壓備選因子的灰色關聯度分析和相關性分析，最終計算結果如表 2所示：

表 2 滲透壓備選因子相關性分析結果統計Table 2 Correlation analysis of osmotic pressure alternative factors

GPS01和GPS02的變化項位移與地下水位備選因子的灰色關聯度分析和相關性分析，最終計算結果如表 3所示：

表 3 地下水位備選因子相關性分析結果統計Table 3 Correlation analysis of underwater alternative factors

由圖 9a可見，HH02對外部降雨有明顯的響應變化，呈現“下降-上升-下降”的總體趨勢，且隨時間推移，土體含水率變化范圍不斷變小，趨于平穩。HH02在2013年、2014年、2015年和2016年的極值差(最大含水率-最小含水率)為1.56， 1.24， 1.18和0.45。HH03和HH05的兩條監測曲線，如圖 9b所示，自監測裝置安裝后即開始走高至2014年3月保持平穩，后小幅回落，于2015年6月再次抬升，結合降雨監測數據未發現明顯規律，與滑坡位移變化曲線無明確響應關系。

GPS01和GPS02的變化項位移與土體含水率備選因子的相關性分析，最終計算結果如表 4所示：

表 4 土體含水率備選因子相關性分析結果統計Table 4 Correlation analysis of alternative factors of soil moisture content

綜合灰色關聯度分析和相關性分析，確定GPS01和GPS02的坡體內部影響因子為SY01、DX06和HH02。

2.3.3 坡體自身穩定性因子

滑坡的變化項位移的大小不僅與坡體內外影響因子的作用有關，也要考慮滑坡自身的穩定性。滑坡越穩定，對坡體內外影響因子的響應就越弱，反之亦然。結合前述對垮梁子滑坡變形趨勢分析，可以明顯地發現滑坡在雨季位移速率增大，而在非雨季期間位移速率較小且平穩，周期內的總體變化趨勢大致一致，但位移量值存在差異。因此，本文將滑坡的變形趨勢過程以年(當年6月至次年5月)為周期進行劃分，周期內滑坡位移量越小，滑坡越穩定。統計周期內單月滑坡位移量，并逐月疊加，得到監測點周期內位移累計增量，以此作為衡量坡體自身穩定性因子。經計算該因子與變化項位移的CG為0.570，CP為0.601，Ct為0.001，因而可以作為滑坡變形的主控因子之一。

2.3.4 IAGA-WNN模型預測結果

基于滑坡累計位移與主控因子的確定，進一步采用IAGA-WNN模型對滑坡變化項位移進行預測，模型預測過程如下：

(1)框架搭建。以2013年6月至2016年5月共計36組監測數據作為訓練學習樣本，以2016年6月至2017年3月共計10組為測試樣本。模型采用3層網絡設置，輸入降雨影響因子、滲透壓影響因子、地下水位影響因子、土體含水率影響因子、坡體自身穩定性因子共計5個節點，并將所有影響因子及其變化項位移進行歸一化處理。

(2)網絡參數的確定。依據國際上較為認可的確定隱含層個數的公式(楊濤， 2012)：

(14)

式中，m、n分別為輸入層和輸出層的個數。

取m=5，n=1，計算可得S取在5～11之間，最終通過試算確定網絡結構為5-7-1，即輸入層節點為5個，隱含層節點為7個，輸出層節點為1個，學習率Ir1=0.01，Ir2=0.001，遺傳優化算法中種群規模P=20，最大進化代數為100代，交叉率Pc=0.3，變異率Pm=0.05，誤差精度為0.001。

(3)模型預測與對比。用已確定的網絡參數對訓練樣本進行學習，得到預測結果，并將結果與傳統BP神經網絡模型和小波神經網絡模型進行對比，結果如圖 10所示。

圖 10 變化項位移預測及對比Fig. 10 Perdiction and comparison of variable displacementsa. 監測點GPS01；b. 監測點GPS02

表 5 各模型預測精度及誤差Table 5 Contrast of accuracy and error for prediction of each model

各模型預測精度及誤差對比如表 5所示，結果表明：IAGA-WNN模型在GPS01和GPS02兩個監測點的確定系數R2分別達到了0.976和0.943，均方根誤差RMSE分別為4.297和6.916，最大相對誤差分別MaxE分別為5.931imm和5.975imm，總體精度較高，與傳統BP神經網絡模型和小波神經網絡相比有明顯優勢。且與GA-WNN模型進行效率對比，耗時相差無幾，R2值分別提高了2.5%和2.1%，RMSE和MaxE也相應降低。

2.4 滑坡累計位移預測結果與誤差分析

通過疊加GM(1， 1)灰色模型得到的趨勢項位移與IAGA-WNN模型得到的變化項位移，最終獲得聯合模型的滑坡累計位移預測值，如圖 11所示。模型預測值變化趨勢與實際監測值較為一致，總體預測精度較好。GPS01模型預測值的確定系數R2值和均方根誤差RMSE分別為0.968， 27.689，GPS02模型預測值的確定系數R2值和均方根誤差RMSE分別為0.971， 26.383。

圖 11 滑坡實際監測值與模型預測值的比較Fig. 11 Comparison of the observing displacements and the predicted valuesa. 監測點GPS01；b. 監測點GPS02

由圖 11可見，誤差主要發生在2013年9月和2013年10月，分析其原因可能是由于坡體內影響因子與滑坡位移變化協同具有一定的滯后性，由外部強降雨引起的坡體位移速率較大幅度增快，經過一段時間后才會反應在坡體內部監測點的數據變化上，因此模型的最終預測結果上產生了誤差較大的點。

3 結 論

本文以四川德陽中江縣垮梁子滑坡為例，從滑坡監測數據著手，優化篩選出滑坡位移關鍵控制因子，應用GM-IAGA-WNN聯合模型對滑坡進行位移預測研究。取得以下主要結論：

(1)通過對監測數據的分析，發現了坡體各部位變化大致同步，但坡體中部變形較坡體前緣更劇烈。坡體變形主要受雙月累計降雨量、SY01、DX06和HH02的控制。因此對于此類“平推式滑坡”的監測重點應放在外部降雨量、坡體后緣滲透壓、中部土體含水率和前緣地下水位上。

(2)本文提出的聯合模型，與傳統BP神經網絡模型、未添加遺傳算法的小波神經網絡模型和未優化遺傳算法-小波神經網絡聯合模型預測結果對比，具有穩定性強、精度高的優勢。盡管由于外部強降雨引起坡體內部影響因子變化具有一定的滯后現象，個別時間點的位移預測值未達到預期預測，但是模型有效降低了預測誤差，總體效果較好，在實際運用中具有較高的推廣應用價值。

(3)通過對監測數據的挖掘，揭示了滑坡變形與坡體內外影響因子的相互關系。為進一步開展平推式滑坡成因機理深入分析和滑坡預警預報提供了一定的技術支撐。