考慮降雨滯后效應(yīng)的八字門滑坡位移預測研究

尚 敏，張惠強，廖 芬，熊德兵，趙國飛

（1.湖北省地質(zhì)災害防治工程技術(shù)研究中心，湖北宜昌 443002；2.三峽庫區(qū)地質(zhì)災害教育部重點實驗室，湖北宜昌 443002；3.湖北長江三峽滑坡國家野外科學觀測研究站，湖北宜昌 443002；4.三峽大學土木與建筑學院，湖北宜昌 443002）

引言

我國地質(zhì)災害發(fā)生十分頻繁，滑坡是我國數(shù)量最多危害最嚴重的地質(zhì)災害類型［1－2］。根據(jù)2020年全國地質(zhì)災害通報數(shù)據(jù)：2020年我國共發(fā)生7 840起地質(zhì)災害，成功預報地質(zhì)災害534起，避免可能傷亡1.8萬人，避免直接經(jīng)濟損失10.2億元，其中，滑坡共發(fā)生4 810起，占比61.35%［3］。對于三峽庫區(qū)滑坡來講，變形趨勢明顯的滑坡處于長期變形預警監(jiān)測中，深入分析滑坡位移監(jiān)測數(shù)據(jù)與其影響因素之間的關(guān)系，開展滑坡預測預報研究，是判斷滑坡穩(wěn)定性的重要方法，也是實現(xiàn)防災減災的重要途徑［4－6］。

八字門滑坡是三峽庫區(qū)重點專業(yè)監(jiān)測滑坡之一，目前主要受到庫水位升降和降雨的周期性變化影響，該滑坡的累積位移－時間曲線呈臺階狀，是一個典型的階躍型滑坡［7－9］，至今仍處于間歇性蠕滑變形狀態(tài)。眾多學者對八字門滑坡開展了位移預測研究，周超等提出了一種基于誘發(fā)因素響應(yīng)分析的進化支持向量機位移預測模型，以當月降雨量、前兩月累積降雨量、當月庫水位平均值、當月庫水位最大變化量為影響因素開展位移預測，認為該方法適合于的階躍型滑坡位移預測［10］。尚敏等［11］認為八字門滑坡累積位移曲線的階躍性與降雨量、庫水位下降速率呈正相關(guān)，建立了一元線性回歸模型對滑坡階躍段的累積位移開展預測，取得了不錯的效果。梁陽等［12］認為八字門滑坡的當月降雨量會影響下一月的位移，降雨響應(yīng)的滯后期約為一個月。以上對研究數(shù)據(jù)的選擇通常以月為單位，時間跨度較大，在數(shù)據(jù)的選擇上存在較大的精度誤差，難以給出更為準確的滯后期。因此，提升數(shù)據(jù)的研究精度以及考慮降雨滯后效應(yīng)的滑坡位移預測是一個重要課題。

降雨是誘發(fā)滑坡變形破壞的重要因素，從時間的角度分析，滑坡的變形破壞多發(fā)生在降雨的中后期或略滯后［13－14］。田正國等［15］對三峽庫區(qū)滑坡誘發(fā)因素分析認為：降雨會導致滑坡體物質(zhì)軟化、抗剪強度降低、重力增大、動水壓力變化等效應(yīng)，由于滑坡體物質(zhì)結(jié)構(gòu)和滲透性不同，滑坡變形對不同強度和時長的降雨響應(yīng)存在不同的滯后時間，但沒有對降雨滯后效應(yīng)開展定量分析。張珍等［16］對重慶市的降雨和滑坡特征資料分析，發(fā)現(xiàn)滑坡變形量與降雨量大小正相關(guān)，4天內(nèi)的累積降雨量對滑坡發(fā)生與否和滑坡發(fā)生數(shù)量的相關(guān)程度很高，夏夢想等［18］對張家界地區(qū)降雨和滑坡關(guān)系分析認為：滑坡的發(fā)生與前10日的降雨均顯著相關(guān)，與前3日累積降雨量的相關(guān)程度最高，但對每個單體滑坡來講，滑坡形態(tài)、物質(zhì)組成和所處環(huán)境不盡相同，缺乏從降雨誘發(fā)滑坡機理上的解釋。唐棟等［19］基于三峽庫區(qū)降雨數(shù)據(jù)，運用Geo-Studio軟件建立邊坡模型，分析了不同前期降雨工況下邊坡的安全系數(shù)，認為前期累積降雨量可以作為判斷邊坡最小安全系數(shù)出現(xiàn)時刻的依據(jù)，砂土邊坡為前期10 d的累積降雨量，黏土邊坡為前期15 d的累積降雨量。但該研究僅基于數(shù)值模擬軟件，未從物理模型實驗的角度開展驗證分析。

以上學者在八字門滑坡位移預測、滑坡與降雨滯后效應(yīng)兩方面開展了研究。但以天為時間單位考慮降雨滯后效應(yīng)的位移預測研究較少，本文將基于八字門滑坡日降雨量和自動位移監(jiān)測數(shù)據(jù)分析降雨滯后期，采用灰色神經(jīng)網(wǎng)絡(luò)預測模型，對八字門滑坡開展位移預測研究，以期進一步提高八字門滑坡位移預測的精確度。

1 八字門滑坡概況

八字門滑坡地處三峽庫區(qū)香溪河河口右岸，地理坐標：110°45′30″E，30°58′16″N?；缕矫娉蚀榛?，主滑方向112°，剖面呈階梯狀，上陡下緩，地表坡度0～30°，由二級平臺組成，前緣平臺高程139～165 m，后緣平臺高程220～230 m［20］，前緣55～156 m被庫水淹沒，水上部分滑坡體長度80 m，寬度100～500 m，厚度10～35 m，體積約200萬立方米，是一發(fā)育于斜逆向巖層的堆積層滑坡，滑床基巖為侏羅紀早期紫紅色砂頁巖、棕色砂巖及燧石礫石，產(chǎn)狀：292∠29°。詳見圖1、圖2所示。八字門滑坡于1981年葛洲壩蓄水后復活，進入初期復活階段；受1998年長時間暴雨影響，滑坡變形加劇趨于整體復活，以牽引式為主；自2003年三峽大壩蓄水以來，滑坡前緣已被庫水淹沒，滑坡變形模式由牽引式轉(zhuǎn)變?yōu)橥埔剖剑?2］。經(jīng)現(xiàn)場調(diào)查可知：新的宏觀變形往往發(fā)生在每年的6～8月份，主要分布于滑坡中部、滑坡后緣等位置，坡體中后緣公路路面已受到嚴重破壞，影響通行。

圖1 八字門滑坡地質(zhì)圖［11］Fig.1 Geological map of Bazimen landslide

圖2 八字門滑坡剖面圖［21］Fig.2 Cross-section of Bazimen landslide

2 監(jiān)測數(shù)據(jù)分析

八字門滑坡于2016年5月安裝完成全自動地表位移監(jiān)測點GSCX3、GSCX5和、ZGX111，實現(xiàn)了監(jiān)測數(shù)據(jù)的實時采集，相對于以往每月1～2次的人工監(jiān)測，數(shù)據(jù)更加豐富全面，為提高原始數(shù)據(jù)的研究精度創(chuàng)造了條件。圖3為八字門滑坡自動位移監(jiān)測點的時間－累積位移－庫水位－降雨量曲線，分析可知：在每年降雨量較小的時間段（2016年10月～2017年3月、2017年11月～2018年5月，2018年8月～2019年3月，2019年7月～2020年4月），滑坡累積位移增長緩慢；而在每年降雨量較大且集中的時間段（2016年5～8月、2017年5月7月10月、2018年6～7月、2019年6～8月、2020年6～8月），累積位移增長迅速呈現(xiàn)階梯狀，且每年“階躍期”和集中降雨期存在時間上的滯后性，即滑坡位移“階躍”發(fā)生在連續(xù)降雨或強降雨之后，這種規(guī)律具有統(tǒng)一性。庫水位與滑坡累積位移曲線的“階躍段”出現(xiàn)規(guī)律關(guān)系并不明顯，在2016年6月～7月庫水為上升期、2017年3月～6月庫水位下降期、2018年6月庫水位低位運行期、2020年6～8月庫水位波動期，累積位移均迅速增加。在2017年8月～10月庫水位上升期內(nèi)，持續(xù)性強降雨出現(xiàn)在10月，但位移“階躍”僅出現(xiàn)在10月，由此可知引起八字門滑坡位移變化的主要因素是強降雨或者持續(xù)性降雨，次要因素是庫水位變化。

圖3 八字門滑坡自動監(jiān)測點時間-累積位移-庫水位-降雨量曲線圖Fig.3 Time-cumulative displacement-reservoir water level-rainfall curve of automatic monitoring point of Bazimen landslide

3 降雨響應(yīng)位移的分布滯后模型

3.1 數(shù)據(jù)選擇與模型建立

八字門滑坡的階躍變形對降雨的響應(yīng)往往不是瞬間發(fā)生，存在時間上的滯后性，即滑坡位移的變化不僅受到同期降雨作用的影響，還受到了前期降雨作用的影響。為了深入研究滑坡位移和降雨量的關(guān)系，選取2016年5月16日～2019年8月31日期間全自動位移數(shù)據(jù)和同期降雨數(shù)據(jù)建立分析，為保證研究精度，選擇天為時間單位，詳見圖3。

分布滯后模型是指被解釋變量的現(xiàn)期值不僅依賴于解釋變量現(xiàn)期值，而且依賴于該變量的過去值或滯后值影響的動態(tài)模型［23］?；挛灰坪徒涤甑捻憫?yīng)關(guān)系與之具有一致性，故將該模型應(yīng)用到降雨響應(yīng)位移的關(guān)系研究中。首先，定義日降雨量為解釋變量；滑坡日位移為被解釋變量，建立如下分布滯后模型：

式中：t=1,2,…,n；s代表滯后期數(shù)；α為模型系數(shù)；Y t表示位移現(xiàn)期值；X t表示降雨量現(xiàn)期值，X t-s表示降雨量的前期值；βs（s=0，1，2..s）為滯后系數(shù)，即滯后s期的降雨量對滑坡位移現(xiàn)期值的影響程度；u t代表誤差項，u t～(0,σ2)。

分布滯后模型從形式可看作是一個多元回歸模型，可采用最小二乘法估計解釋變量的系數(shù)β0～βs［24］。

3.2 模型的估計與檢驗

在分布滯后模型參數(shù)估計時，首先需要確定滯后期數(shù)，現(xiàn)式估計法適用于求解分布滯后模型中滯后期未知的情況，即通過依次增加滯后期數(shù)，直到滯后變量的回歸系數(shù)在統(tǒng)計上不再顯著時，便可確定滯后期數(shù)。為了簡化該過程采用皮爾遜相關(guān)性分析初步確定滯后期數(shù)［25］，結(jié)果詳見表1～表2。

表1 皮爾遜相關(guān)系數(shù)及相關(guān)強度表Table 1 Pearson correlation coefficient and correlation intensity table

表2 日降雨量與日位移的相關(guān)性系數(shù)表Table 2 Correlation coefficient between daily rainfall and daily displacement

表2中s=0代表當天降雨量與當天位移的相關(guān)系數(shù)；s=1表示前第1 d的降雨量與當天位移的相關(guān)系數(shù)；s=2表示前第2 d的降雨量與當天位移的相關(guān)系數(shù)，依次類推。分析可知：s=0～12時，相關(guān)系數(shù)先增大后減小，呈現(xiàn)中等強度至弱相關(guān)；s=13時，相關(guān)性系數(shù)小于0.2，呈極弱相關(guān)或無相關(guān)，即滑坡當天位移與前12 d降雨均具相關(guān)性，根據(jù)相關(guān)程度初步確定位移滯后降雨的滯后期為12 d。

以GSCX3監(jiān)測點為例，利用Eviews軟件的Estimate Equation命令，對滯后期為12的解釋變量參數(shù)估計并建立模型，最后對其進行方程顯著性和變量顯著性檢驗，結(jié)果見表3。

表3 日位移與日降雨量分布滯后模型結(jié)果表Table 3 Results of daily displacement and daily precipitation distribution lag model

由表3可知，日降雨量序列數(shù)據(jù)X t-s及常數(shù)α顯著性概率均小于0.05，置信度超過95%，變量顯著性檢驗通過；該模型的F統(tǒng)計值為90.379，F(xiàn)統(tǒng)計值的顯著性概率Prob（F-statistic）=0.000＜0.05，方程顯著性檢驗亦通過。增加滯后期到13期進行參數(shù)估計并檢驗，X t-13的顯著性概率Prob（t-statistic）=0.4067＞0.05，不滿足變量顯著性檢驗要求，由此確定滯后期為12 d。最后通過最小二乘估計得出滯后期為12 d的分布滯后模型，見式（2），在統(tǒng)計意義上是顯著的，其結(jié)果能代表實際滑坡位移的變化規(guī)律。

3.3 模型的驗證

將八字門滑坡2019年08月18日至2020年08月31日的實測日降雨量數(shù)據(jù)代入（式2），對2019年09月01日至2020年08月31日的日變形位移進行估算，并與位移實際值對比驗證該分布滯后模型的準確性，結(jié)果如圖4所示。

圖4 八字門滑坡分布滯后模型估算值與實際監(jiān)測值對比分析圖Fig.4 Comparison and analysis diagram of distributed lag model estimation and actual monitoring values of Bazimen landslide

通過圖4可發(fā)現(xiàn)，在2020年6～8月份雨量集中時，分布滯后模型位移估算值與監(jiān)測值的變化趨勢和突變點的位置一致，估算值曲線前移（如圖中黑色方框標示），說明在一定程度上抵消了降雨與位移之間的滯后效應(yīng)。該模型反映的降雨響應(yīng)滑坡位移變化的滯后關(guān)系是較為準確的，但估算值與監(jiān)測值在峰值處仍有較大偏差，這是因為降雨分布滯后模型只考慮了降雨因素，實際上滑坡變形還受到其他因素的影響，故需要進一步提高位移預測的精度。

4 八字門滑坡位移精確預測

4.1 預測方法的確定

八字門滑坡日變形量受到前期12天降雨的影響，將降雨滯后效應(yīng)考慮到滑坡日位移預測之中，可以提高滑坡預測精度，在滑坡位移預測中考慮滑坡的誘發(fā)因素也可以提高預測精度［26－27］，周超等基于滑坡誘發(fā)因素對八字門滑坡位移預測取得了很好的預測效果［10］，文章選取滯后期內(nèi)12 d累積降雨量、庫水位漲落速率、滑坡自身位移作為滑坡變形的影響因素，對八字門滑坡的日變形位移進行預測研究。

八字門滑坡的自動位移只有最近4年的數(shù)據(jù)，數(shù)據(jù)樣本量較少，滑坡變形受到多種因素影響呈非線性特征，灰色模型具備小樣本數(shù)據(jù)建模優(yōu)點［28］，神經(jīng)網(wǎng)絡(luò)適用于非線性和不確定性系統(tǒng)預測，灰色神經(jīng)網(wǎng)絡(luò)組合模型結(jié)合了灰色系統(tǒng)和神經(jīng)網(wǎng)絡(luò)的優(yōu)點，能夠?qū)Χ喑梢虻膹碗s未知系統(tǒng)進行建模，因此選擇灰色神經(jīng)網(wǎng)絡(luò)模型對八字門滑坡日變形位移的預測研究。

4.2 灰色神經(jīng)網(wǎng)絡(luò)基本原理

灰色模型是指對原始數(shù)據(jù)序列x(0)一次累加后得到新序列x(1)，新序列規(guī)律性會明顯增強，然后對該新序列進行擬合的模型［29］。

設(shè)時間數(shù)據(jù)序列x(0)：

累加之后的新序列x(1)為：

根據(jù)新的序列x()1，建立白化微分方程，即：

該方程的解為：

是x t(1)的估計值，對其做一次累減得到x(0)的預測值x t*(0)：

灰色神經(jīng)網(wǎng)絡(luò)模型是在灰色理論基礎(chǔ)上構(gòu)建的，為了方便表述，重新定義符號，用x(t)表示x(0)，y(t)代表一次累加后生成的數(shù)列x t(1)，z(t)表示預測結(jié)果x t*(1)，含有n個參數(shù)的灰色神經(jīng)網(wǎng)絡(luò)模型的微分方程：

式中，y2,...y n為系統(tǒng)輸入?yún)?shù)；y1為系統(tǒng)輸出參數(shù)；a,b1,b2,...b n-1為微分方程系數(shù)；微分方程式的解為：

式（10）映射到神經(jīng)網(wǎng)絡(luò)中，可得到含有n個輸入?yún)?shù)、1個輸出參數(shù)的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，如圖5所示。

圖5 灰色神經(jīng)網(wǎng)絡(luò)模型拓撲結(jié)構(gòu)Fig.5 Grey neural network modeltopology

其中t為輸入?yún)?shù)序號；y2(t)，...y n(t)為輸入?yún)?shù)；ω21,ω22,...,ω2n,ω31,ω32,...,ω3n為神經(jīng)網(wǎng)絡(luò)權(quán)值，LA、LB、LC、LD為灰色神經(jīng)網(wǎng)絡(luò)的層級，y1為網(wǎng)絡(luò)預測值。

4.3 基于灰色神經(jīng)網(wǎng)絡(luò)的位移預測

選擇八字門滑坡GSCX3全自動監(jiān)測點2016年5月24日－2020年6月12日數(shù)據(jù)（日變形位移、庫水位漲落速率、12 d累積降雨量）作為訓練樣本，對2020年6月24日～2020年8月28日的位移進行預測，借助Matlab編寫程序處理數(shù)據(jù)，訓練樣本原始數(shù)據(jù)見表4。

表4 訓練樣本原始數(shù)據(jù)表Table 4 Raw data table of training samples

首先，對輸入原始數(shù)據(jù)（訓練樣本）進行歸一化處理，消除數(shù)據(jù)之間數(shù)量級的差異，歸一化函數(shù)形式如下：

式中：xmin是x k序列的最小值；xmax是序列的最大值。

其次，構(gòu)建灰色神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)，位移預測輸入數(shù)據(jù)為3維，輸出數(shù)據(jù)為1維，所以灰色神經(jīng)網(wǎng)絡(luò)模型網(wǎng)絡(luò)結(jié)構(gòu)為1－1－4－1，即LA層有1個節(jié)點，時間序列t，LB層有1個節(jié)點，LC層有4個節(jié)點，在第2～4個節(jié)點分別輸入12天累積降雨量、庫水位漲落速率、滑坡自身位移的歸一化數(shù)值，輸出數(shù)據(jù)為預測位移。

然后，設(shè)置網(wǎng)絡(luò)訓練參數(shù)，對灰色神經(jīng)網(wǎng)絡(luò)進行訓練，使網(wǎng)絡(luò)具有滑坡日變形位移預測的能力，在訓練過程中網(wǎng)絡(luò)的連接權(quán)值也將不斷地被調(diào)整修正。

LD層中神經(jīng)元閾值為：

最后，用訓練好的灰色神經(jīng)網(wǎng)絡(luò)預測結(jié)果輸出，將輸出層進行一次累減并反歸一化得到預測值，預測結(jié)果分析詳見圖6和表5。

表5 八字門滑坡部分日變形位移預測分析表Table 5 Prediction table of partial daily deformation and displacement of Bazimen landslide

圖6 八字門滑坡日位移實際值與預測值對比分析圖Fig.6 Comparison and analysis diagram of actual and predicted daily displacement of Bazimen landslide

由圖6可知，考慮降雨滯后效應(yīng)、庫水位漲落速率、滑坡自身位移作為滑坡位移影響因素的灰色神經(jīng)網(wǎng)絡(luò)預測模型，能夠準確地預測滑坡日位移變化趨勢和突變點，計算可知平均絕對百分比誤差（MAPE）為11.98%，精確度88.02%，預測值與實際監(jiān)測值的最大偏差小于2 mm，符合模型預測的精度要求。可為八字門滑坡的監(jiān)測和防災減災工作提供理論支撐和參考建議，同時可為三峽庫區(qū)同類型滑坡的預測研究提供借鑒。

5 結(jié)論

（1）八字門滑坡的累積位移“階躍”發(fā)生在每年的降雨集中時期，八字門滑坡的變形主要受強降雨或持續(xù)性降雨影響。

（2）采用日降雨量和日變形位移數(shù)據(jù)建立了八字門滑坡位移分布滯后模型，模型估計值與滑坡實際監(jiān)測值變化趨勢一致，且估計值曲線前移，在一定程度上滯后效應(yīng)已被消除，八字門滑坡日變形位移響應(yīng)降雨的滯后期為12 d的分布滯后模型是較為準確的。

（3）在灰色神經(jīng)網(wǎng)絡(luò)預測模型考慮降雨滯后效應(yīng)，利用滯后期內(nèi)12 d的累積降雨量、庫水位漲落速率、滑坡自身位移作為影響因子，對八字門滑坡日變形位移預測，能夠準確地預測位移的變化趨勢及突變點。平均絕對百分比誤差11.98%，精確度88.02%，最大偏差2 mm?？蔀榘俗珠T滑坡的預警預報和防治提供可靠建議。