梁 陽(yáng) 肖 婷 胡 程 任世聰 曾 亮

（1.重慶三峽學(xué)院電子與信息工程學(xué)院，重慶 404000；2.北京理工大學(xué)重慶創(chuàng)新中心，重慶 401135；3.北京理工大學(xué)信息與電子學(xué)院雷達(dá)技術(shù)研究所，北京 100081；4.重慶市地質(zhì)礦產(chǎn)研究院，重慶 400042）

1 引言

我國(guó)滑坡災(zāi)害具有分布廣、頻率高和危害大的特點(diǎn)，每年都會(huì)造成不同程度的人員傷亡和嚴(yán)重的經(jīng)濟(jì)損失，因此需要加強(qiáng)對(duì)滑坡的預(yù)測(cè)預(yù)警。滑坡地質(zhì)災(zāi)害的形成原因較為復(fù)雜，滑坡影響因素多、隨機(jī)性強(qiáng)、預(yù)測(cè)難度大，難以及時(shí)定點(diǎn)地準(zhǔn)確預(yù)報(bào)。因此，最好的方式是在監(jiān)測(cè)地質(zhì)環(huán)境變化的同時(shí)結(jié)合現(xiàn)有的長(zhǎng)期數(shù)據(jù)進(jìn)行有效地預(yù)測(cè)，對(duì)滑坡環(huán)境的變化情況和發(fā)展趨勢(shì)進(jìn)行分析，同時(shí)計(jì)算氣象環(huán)境變化與滑坡發(fā)生的內(nèi)在聯(lián)系，進(jìn)行有效地提前預(yù)防。

滑坡預(yù)測(cè)預(yù)報(bào)研究發(fā)展階段大致可分為20 世紀(jì)60 年代的自然現(xiàn)象，經(jīng)驗(yàn)方程預(yù)測(cè)時(shí)期，20 世紀(jì)80 年代的統(tǒng)計(jì)分析預(yù)測(cè)時(shí)期和20 世紀(jì)90 年代的非線性及綜合預(yù)測(cè)時(shí)期［1］。進(jìn)入21 世紀(jì)，隨著多項(xiàng)式擬合、神經(jīng)網(wǎng)絡(luò)、極限學(xué)習(xí)機(jī)、機(jī)器學(xué)習(xí)等智能算法的發(fā)展和逐步成熟，其在非線性映射能力和高精度函數(shù)逼近等方面表現(xiàn)出優(yōu)越性能，于是地質(zhì)災(zāi)害預(yù)警領(lǐng)域也開始大量應(yīng)用［2］。其中，靜態(tài)預(yù)測(cè)模型多用于以往的實(shí)驗(yàn)中，雖然其表現(xiàn)出良好的非線性映射能力，但將滑坡位移預(yù)測(cè)視為靜態(tài)回歸問題忽略了滑坡演化的動(dòng)態(tài)系統(tǒng)本質(zhì)［3］。為了避免穩(wěn)態(tài)模型對(duì)滑坡預(yù)測(cè)造成的這一問題，可以考慮采用動(dòng)態(tài)模型來(lái)對(duì)滑坡進(jìn)行預(yù)測(cè)。在目前的動(dòng)態(tài)系統(tǒng)模型中，RNN（Recurrent Neural Network）是一種經(jīng)典的算法，主要用于處理時(shí)間序列的數(shù)據(jù)，已大量應(yīng)用于搜索查詢、股票市場(chǎng)分析和用戶視頻推送等方面［4］。但是由于RNN 內(nèi)部結(jié)構(gòu)的問題，在傳播過程中會(huì)出現(xiàn)梯度消失或梯度爆炸問題。由此，學(xué)術(shù)界提出了長(zhǎng)短期記憶（Long Short-Term Memory，LSTM）神經(jīng)網(wǎng)絡(luò)來(lái)解決這個(gè)問題。在傳播過程中會(huì)出現(xiàn)梯度消失或梯度爆炸問題，而長(zhǎng)短期記憶（Long Short-Term Memory，LSTM）神經(jīng)網(wǎng)絡(luò)的出現(xiàn)有效地解決了這個(gè)問題，相比于RNN，LSTM 結(jié)構(gòu)中新增了遺忘門、輸入門和輸出門等單元［5-6］，解決了神經(jīng)網(wǎng)絡(luò)固有的易陷入局部最小值、梯度缺失問題［7］。因此，本文選取LSTM 和RNN 兩個(gè)神經(jīng)網(wǎng)絡(luò)模型進(jìn)行位移預(yù)測(cè)及模型性能對(duì)比。

以八字門滑坡為例，先將滑坡點(diǎn)總位移分解為受自身地質(zhì)條件演化控制的趨勢(shì)項(xiàng)和受降雨、庫(kù)水位變化等外界誘發(fā)因素影響的的周期項(xiàng)；然后，采用多項(xiàng)式最小二乘法擬合預(yù)測(cè)趨勢(shì)項(xiàng)位移，利用LSTM 和RNN神經(jīng)網(wǎng)絡(luò)模型預(yù)測(cè)周期項(xiàng)位移。對(duì)總預(yù)測(cè)值進(jìn)行誤差分析，結(jié)果表明LSTM 模型通過控制歷史信息的流向，在面對(duì)數(shù)十年的長(zhǎng)時(shí)間序列時(shí)表現(xiàn)出明顯更好的預(yù)測(cè)能力。

2 滑坡位移時(shí)間序列分析

2.1 滑坡累計(jì)位移分解

研究表明，滑坡體和周圍不動(dòng)巖體自身的層位，地質(zhì)結(jié)構(gòu)結(jié)合其周圍的氣象環(huán)境條件共同導(dǎo)致了滑坡體的變形位移，其位移表現(xiàn)為非線性函數(shù)，滑坡累計(jì)位移時(shí)間曲線可表達(dá)為：

式中：X(t)為位移時(shí)間序列；s(t)為趨勢(shì)項(xiàng)位移，是滑坡位移序列變動(dòng)總方向，它受滑坡體和周圍不動(dòng)巖體的自身地質(zhì)結(jié)構(gòu)影響；c(t)為周期項(xiàng)位移，它是由降雨、庫(kù)水位變化等因素決定的；ε(t)為隨機(jī)項(xiàng)位移，一般為監(jiān)測(cè)點(diǎn)的局部擾動(dòng)等，視為曲線噪音，一般不考慮［8］。

2.2 模型介紹

（1）循環(huán)神經(jīng)網(wǎng)絡(luò)模型（RNN）

RNN 模型通過對(duì)序列型的數(shù)據(jù)進(jìn)行高度擬合從而建立深度模型。在傳統(tǒng)的全連接神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)中，神經(jīng)元之間互不影響，并沒有直接聯(lián)系，神經(jīng)元與神經(jīng)元之間相互獨(dú)立。而在RNN 結(jié)構(gòu)中，隱藏層的神經(jīng)元開始通過一個(gè)隱藏狀態(tài)所相連，通常會(huì)被表示為h(t)。其結(jié)構(gòu)圖如圖1所示。

上圖為經(jīng)典的RNN 結(jié)構(gòu)，其運(yùn)算過程可以表示為：

式中y(t)表示神經(jīng)網(wǎng)絡(luò)的輸出，h(t-1)表示前一個(gè)時(shí)間點(diǎn)的狀態(tài)；x(1)，x(2)...x(t)為序列數(shù)據(jù)，即神經(jīng)網(wǎng)絡(luò)的輸入；U、W、V是參數(shù)矩陣，b、c是偏置項(xiàng)，f是激活函數(shù)，通常采用”熱擼”、tanh 函數(shù)作為激活函數(shù)，用Softmax將輸出轉(zhuǎn)換成各個(gè)類別的概率。

RNN 在反向傳導(dǎo)時(shí)，需要計(jì)算各個(gè)參數(shù)的梯度，在梯度一步步傳播過程中，隨著參數(shù)U、W、V的主特征值在1 附近的變化與積累，最終會(huì)導(dǎo)致梯度的消失或爆炸，無(wú)法學(xué)到遠(yuǎn)距離的數(shù)據(jù)。所以RNN單元在面對(duì)長(zhǎng)序列數(shù)據(jù)時(shí)，很容易便遭遇梯度彌散，使得RNN 只具備短期記憶，即RNN 面對(duì)長(zhǎng)序列數(shù)據(jù)，僅可獲取較近的序列的信息，而對(duì)較早期的序列不具備記憶功能，從而丟失信息。因此，Hochreiter 和Schmidhuber 提出了一種改進(jìn)方法，即長(zhǎng)短時(shí)記憶網(wǎng)絡(luò)（LSTM）9］。

（2）長(zhǎng)短時(shí)間記憶網(wǎng)絡(luò)模型（LSTM）

LSTM 是一種特殊的RNN，相比于RNN，LSTM模型新增加了記憶單元、細(xì)胞狀態(tài)、輸入門、遺忘門和輸出門等機(jī)制，從而有效地解決了RNN 的短期依賴瓶頸［10］，其結(jié)構(gòu)圖如圖2所示。

如圖2 所示，LSTM 的關(guān)鍵是整個(gè)單元的狀態(tài)（綠色框代表一個(gè)單元或者一個(gè)細(xì)胞），其中橫穿單元的水平線稱為細(xì)胞鏈，它起著計(jì)算的作用。細(xì)胞中的信息狀態(tài)直接穿過細(xì)胞鏈，在細(xì)胞鏈中存在一些線性相互作用，以保持信息的容易傳遞。使用細(xì)胞鏈的優(yōu)點(diǎn)是可以直接添加或刪除信息。因此，可以選擇性地允許信息通過并與sigmoid 神經(jīng)元層或逐點(diǎn)乘法操作一起使用。LSTM 有三個(gè)值用于保護(hù)和控制單元狀態(tài)的門：

遺忘門：首先，LSTM 是計(jì)算從細(xì)胞狀態(tài)中丟棄哪些信息。這個(gè)決定是由遺忘門做出的。在單元狀態(tài)為Ct-1時(shí)，該門可以獲得值為1 或0 的ht-1（前一層的隱藏狀態(tài)）和xt的輸出。1 表示“完全保留”，0表示“完全放棄”。

輸入門：這個(gè)步驟是計(jì)算有多少新信息應(yīng)該添加到這個(gè)細(xì)胞狀態(tài)中。實(shí)現(xiàn)這一目標(biāo)有兩個(gè)步驟：首先，將前一層隱藏狀態(tài)的信息和當(dāng)前輸入的信息傳遞到sigmoid函數(shù)中去；然后使用前一層隱藏狀態(tài)的信息和當(dāng)前輸入的信息來(lái)計(jì)算需要更新的信息并生成一個(gè)向量。總之，通過組合這兩個(gè)步驟來(lái)更新單元的狀態(tài)。我們將舊狀態(tài)Ct-1乘以ft（ft是需要丟棄的信息），然后Ct-1更新為Ct。ft的方程，Ct-1，Ct為分別如下。

式中，ht-1是前一層的輸出，xt是輸入向量，σ是sigmoid 函數(shù)，Wf是ft和it的系數(shù)矩陣，bf是ft的偏置，bi是it的偏置，Wc是的系數(shù)矩陣，bc是的偏置，Ct-1是來(lái)自前一塊的存儲(chǔ)器，Ct是來(lái)自當(dāng)前塊的存儲(chǔ)器。

輸出門：最終，應(yīng)該確定輸出值，此輸出取決于細(xì)胞狀態(tài)。首先，計(jì)算一個(gè)sigmoid 函數(shù)來(lái)決定輸出單元的哪個(gè)部分，接下來(lái)，通過前一層的隱藏狀態(tài)選擇細(xì)胞狀態(tài)，并將其與sigmoid 函數(shù)的輸出相乘。

式中，Wo是ot的系數(shù)矩陣，bo是ot的偏置，ht是當(dāng)前模塊的輸出。

輸入層有三個(gè)輸入，隱藏層有10 個(gè)神經(jīng)元，輸出層預(yù)測(cè)一個(gè)值。激活函數(shù)用sigmoid 迭代100 次。批量損失函數(shù)為全批量學(xué)習(xí)損失函數(shù)和交叉熵。

式中T是周期時(shí)間，yi是樣本數(shù)據(jù)，ht是LSTM 的輸出。

2.3 模型預(yù)測(cè)精度指標(biāo)

本文采用均方根誤差（Root mean square error，RMSE）和相關(guān)性系數(shù)（Correlation coefficient，R）評(píng)價(jià)模型預(yù)測(cè)精度［11］，絕對(duì)誤差和相對(duì)誤差用于評(píng)價(jià)最終預(yù)測(cè)結(jié)果的誤差。

式中：xi和分別為真實(shí)值和預(yù)測(cè)值；和分別為真實(shí)值和預(yù)測(cè)值的平均值；N為樣本數(shù)。

3 八字門滑坡

3.1 滑坡地質(zhì)概況

八字門滑坡地處湖北省秭歸縣歸州鎮(zhèn)，三峽庫(kù)區(qū)香溪河右岸河口處。香溪河在此呈南北走向，與長(zhǎng)江近直交，三峽水庫(kù)已淹沒滑坡體前緣55～156 m段［12］?；麦w呈撮箕狀展布于岸坡坡腳，分布高程139～280 m 西高東低，向東傾斜［13］。滑坡于1981 年出現(xiàn)復(fù)活跡象，此后幾年周邊相繼出現(xiàn)了4 條大裂縫。1987 年到1989 期間暴雨導(dǎo)致變形加劇，坡面房屋開裂，公路出現(xiàn)裂縫。20 世紀(jì)末長(zhǎng)江爆發(fā)洪水，該區(qū)域186 m 高程處產(chǎn)生多條裂縫，路面下沉，房屋遭到破壞。21 世紀(jì)初，隨著三峽水庫(kù)不斷地蓄水，滑坡各部相繼出現(xiàn)多條裂縫，原有裂縫進(jìn)一步加劇變形。此后滑坡自身的地質(zhì)結(jié)構(gòu)發(fā)生改變，穩(wěn)定性變差，在降雨和庫(kù)水位下降時(shí)期，滑坡多次出現(xiàn)裂縫，公路路面發(fā)生小型坍塌。

3.2 監(jiān)測(cè)數(shù)據(jù)分析

八字門滑坡從20 世紀(jì)80 年代開始，陸續(xù)發(fā)生位移變形現(xiàn)象。在三峽工程逐步推進(jìn)的同時(shí)，庫(kù)水位也逐漸升高，其對(duì)八字門滑坡體及其周圍不動(dòng)巖體的地質(zhì)結(jié)構(gòu)造成了一定程度的影響，導(dǎo)致滑坡的不穩(wěn)定性增加。為了了解掌握滑坡的具體情況和監(jiān)控其位移變形，滑體上設(shè)置有ZG109，ZG110，ZG111，ZG112四個(gè)監(jiān)測(cè)點(diǎn)，如圖3所示。

如圖4 所示，2003～2013 年間4 個(gè)監(jiān)測(cè)點(diǎn)的位移時(shí)間序列。2006 年庫(kù)水位從145 m 升至155 m，導(dǎo)致海拔較低的ZG109 監(jiān)測(cè)點(diǎn)失效。2008 年庫(kù)水位上升至175 m，導(dǎo)致ZG112 監(jiān)測(cè)點(diǎn)失效。ZG110和ZG111 都具有明顯的階躍式位移特性，變形更大的ZG111 監(jiān)測(cè)點(diǎn)能更好地代表滑坡的危險(xiǎn)狀態(tài)，因此，選取該點(diǎn)的位移時(shí)間曲線進(jìn)行位移預(yù)測(cè)分析。監(jiān)測(cè)點(diǎn)數(shù)據(jù)的總時(shí)間跨度為2003 年7 月～2013 年12 月，將2011 年8 月以前的數(shù)據(jù)作為建模數(shù)據(jù)，預(yù)測(cè)2011 年9 月～2012 年12 月的位移值，最后將預(yù)測(cè)值與真實(shí)值進(jìn)行比對(duì)。

圖5為ZG111監(jiān)測(cè)點(diǎn)累計(jì)位移量與庫(kù)水位和降雨量的關(guān)系圖，ZG111 的監(jiān)測(cè)點(diǎn)年變形量在122～268 mm 之間，且變形量逐年增加，其中2008 年到2009 年變化量最大為268 mm。庫(kù)水位自2006 年9 月開始呈周期性波動(dòng)，年變化量在20.14～28.43 m 之間，其中2011 年變化量最大為28.43 m。每年7 月份左右是庫(kù)水位的谷值，12 月份左右是庫(kù)水位的峰值。降雨量呈季節(jié)性周期性變化，每年4到9月份為強(qiáng)降雨季節(jié)，7、8月降雨量最大，1、2月降雨量最小。

4 滑坡位移預(yù)測(cè)

4.1 趨勢(shì)項(xiàng)預(yù)測(cè)分析

將滑坡總位移分解為趨勢(shì)項(xiàng)和周期性，趨勢(shì)項(xiàng)采用一次移動(dòng)平均法計(jì)算。一次移動(dòng)平均法是把時(shí)間序列數(shù)據(jù)按一定的長(zhǎng)度計(jì)算一次移動(dòng)平均，能較好地反映時(shí)間序列的趨勢(shì)及其變化，計(jì)算公式如下：

式中Mi為時(shí)間序列，為原序列中Mt所對(duì)應(yīng)的趨勢(shì)項(xiàng)提取值，t=N，N+1，...，n，分段點(diǎn)數(shù)N為跨越期數(shù)，即由遠(yuǎn)及近用于計(jì)算的數(shù)據(jù)個(gè)數(shù)，用于平滑隨機(jī)性帶來(lái)的偏差，一般N=6～200。本文分析的周期項(xiàng)主要由庫(kù)水位，降雨量等以年為周期進(jìn)行波動(dòng)的因子組成，同時(shí)文章中N以月為單位，所以為了更好地平滑趨勢(shì)項(xiàng)和保留周期項(xiàng)的特征，本文N取12。提取出趨勢(shì)項(xiàng)位移后，利用三次多項(xiàng)式函數(shù)對(duì)其進(jìn)行擬合及預(yù)測(cè)［14］。由于庫(kù)水位從2008 年開始正式周期性波動(dòng)，數(shù)據(jù)以2008 年9 月為界線分兩段進(jìn)行擬合（圖6、圖7），采用第二段的擬合公式進(jìn)行預(yù)測(cè)（公式（17））。趨勢(shì)項(xiàng)預(yù)測(cè)值如圖8 所示，預(yù)測(cè)值與真實(shí)值高度重合。

擬合函數(shù)為：

4.2 周期項(xiàng)預(yù)測(cè)分析

周期項(xiàng)影響因素主要為前期總位移、降雨量和庫(kù)水位，將影響因素初步用10 個(gè)指標(biāo)因子表示，采用神經(jīng)網(wǎng)絡(luò)建立指標(biāo)因子與周期項(xiàng)位移值的關(guān)系模型。為避免不相干因子對(duì)建模產(chǎn)生的干擾，首先需要對(duì)特征因子進(jìn)行篩選。

4.2.1 特征因子篩選

本文通過Pearson相關(guān)性分析來(lái)篩選特征因子，Pearson 相關(guān)系數(shù)是衡量?jī)蓚€(gè)連續(xù)變量間關(guān)系的大小和方向的，其取值范圍為［-1，1］，-1代表負(fù)相關(guān)，1代表正相關(guān)，0則代表不存在相關(guān)關(guān)系。相關(guān)系數(shù)越接近0，相關(guān)關(guān)系越弱；越接近-1 或1，相關(guān)關(guān)系越強(qiáng)，公式如下。

其中，X，Y為兩個(gè)等長(zhǎng)向量，N為向量元素個(gè)數(shù)。

對(duì)2003年7月～2011年8月的103組數(shù)據(jù)（每組都是一個(gè)周期項(xiàng)對(duì)應(yīng)10個(gè)特征因子）作各特征因子與周期項(xiàng)的Pearson相關(guān)性分析，相關(guān)性分析見表1。表中當(dāng)月降雨量、前三月庫(kù)水位變化量?jī)蓚€(gè)特征因子相關(guān)性較低，說(shuō)明降雨入滲、軟化巖土體參數(shù)等存在一定滯后性，近兩月的庫(kù)水位波動(dòng)對(duì)滑坡應(yīng)力狀態(tài)改變明顯。因此，剔除當(dāng)月降雨量和前三月庫(kù)水位變化量，剩下的8個(gè)作為有效特征因子。

表1 影響因子篩選結(jié)果Tab.1 Screening results of influencing factors

4.2.2 建立模型

本文實(shí)驗(yàn)樣本集屬于小規(guī)模樣本集，為了保證預(yù)測(cè)集的樣本數(shù)量滿足實(shí)驗(yàn)要求，把樣本數(shù)據(jù)按6∶2∶2 的比例分為訓(xùn)練集，驗(yàn)證集和預(yù)測(cè)集。其中訓(xùn)練集通過不斷地迭代訓(xùn)練來(lái)得到合理的神經(jīng)網(wǎng)絡(luò)的參數(shù)，驗(yàn)證集用于檢驗(yàn)訓(xùn)練后的模型是否符合要求，測(cè)試集用于評(píng)價(jià)神經(jīng)網(wǎng)絡(luò)的性能。LSTM 和RNN模型預(yù)測(cè)結(jié)果，分別如圖9、圖10所示。對(duì)預(yù)測(cè)段的預(yù)測(cè)值與真實(shí)值進(jìn)行誤差分析可知（表2），LSTM 模型的預(yù)測(cè)結(jié)果誤差約10 mm，相關(guān)系數(shù)R高達(dá)0.985，預(yù)測(cè)值及趨勢(shì)線與真實(shí)情況幾乎吻合（圖9）。RNN模型的預(yù)測(cè)結(jié)果誤差約24 mm，相關(guān)系數(shù)R為0.741，預(yù)測(cè)值及趨勢(shì)線與真實(shí)情況存在一定偏差（圖10）。由預(yù)測(cè)結(jié)果誤差圖（圖11）可以看出，本案例的預(yù)測(cè)結(jié)果LSTM 模型明顯優(yōu)于RNN 模型，因此周期項(xiàng)預(yù)測(cè)結(jié)果采用LSTM 模型所計(jì)算的預(yù)測(cè)值。

表2 誤差分析Tab.2 Error analysis

4.3 累計(jì)位移預(yù)測(cè)

將預(yù)測(cè)得到的趨勢(shì)項(xiàng)與周期項(xiàng)進(jìn)行數(shù)值相加，即得到預(yù)測(cè)的總位移值（圖12）。由圖可知，預(yù)測(cè)曲線與真實(shí)曲線高度吻合，在庫(kù)水位與強(qiáng)降雨期間，準(zhǔn)確預(yù)測(cè)到了位移發(fā)生的時(shí)間和階躍位移量。

5 結(jié)論

文章以滑坡位移這一非線性系統(tǒng)為對(duì)象，采用神經(jīng)網(wǎng)絡(luò)模型，結(jié)合十年監(jiān)測(cè)數(shù)據(jù)，通過選取合適的指標(biāo)因子，對(duì)八字門滑坡ZG111 監(jiān)測(cè)點(diǎn)進(jìn)行位移預(yù)測(cè)，并主要得出如下結(jié)論：

（1）八字門滑坡位移呈臺(tái)階式階躍上升，位移階躍期與強(qiáng)降雨期、庫(kù)水位下降期高度重合，說(shuō)明降雨、庫(kù)水位波動(dòng)是該滑坡位移主要影響因素。

（2）由指標(biāo)因子與周期項(xiàng)的相關(guān)性分析可知，“當(dāng)月降雨量”、“前三月庫(kù)水位變化量”兩個(gè)因子與周期項(xiàng)的相關(guān)性較低，說(shuō)明降雨入滲、軟化巖土體參數(shù)等存在一定滯后性，當(dāng)月降雨量主要影響下月位移，庫(kù)水位波動(dòng)對(duì)滑坡應(yīng)力狀態(tài)的影響時(shí)間為兩個(gè)月。

（3）滑坡的位移不僅受到外界因素的影響（如：降雨、庫(kù)水位），還受到前期位移量的影響，位移時(shí)間序列具有時(shí)間指向性。與傳統(tǒng)的循環(huán)神經(jīng)網(wǎng)絡(luò)（RNN）相比，長(zhǎng)短時(shí)間記憶網(wǎng)絡(luò)（LSTM）模型通過控制歷史信息的流向，更適用于有時(shí)間順序的長(zhǎng)時(shí)間序列。