尚 敏,張惠強,廖 芬,熊德兵,趙國飛
(1.湖北省地質(zhì)災(zāi)害防治工程技術(shù)研究中心,湖北宜昌 443002;2.三峽庫區(qū)地質(zhì)災(zāi)害教育部重點實驗室,湖北宜昌 443002;3.湖北長江三峽滑坡國家野外科學(xué)觀測研究站,湖北宜昌 443002;4.三峽大學(xué)土木與建筑學(xué)院,湖北宜昌 443002)
我國地質(zhì)災(zāi)害發(fā)生十分頻繁,滑坡是我國數(shù)量最多危害最嚴(yán)重的地質(zhì)災(zāi)害類型[1-2]。根據(jù)2020年全國地質(zhì)災(zāi)害通報數(shù)據(jù):2020年我國共發(fā)生7 840起地質(zhì)災(zāi)害,成功預(yù)報地質(zhì)災(zāi)害534起,避免可能傷亡1.8萬人,避免直接經(jīng)濟損失10.2億元,其中,滑坡共發(fā)生4 810起,占比61.35%[3]。對于三峽庫區(qū)滑坡來講,變形趨勢明顯的滑坡處于長期變形預(yù)警監(jiān)測中,深入分析滑坡位移監(jiān)測數(shù)據(jù)與其影響因素之間的關(guān)系,開展滑坡預(yù)測預(yù)報研究,是判斷滑坡穩(wěn)定性的重要方法,也是實現(xiàn)防災(zāi)減災(zāi)的重要途徑[4-6]。
八字門滑坡是三峽庫區(qū)重點專業(yè)監(jiān)測滑坡之一,目前主要受到庫水位升降和降雨的周期性變化影響,該滑坡的累積位移-時間曲線呈臺階狀,是一個典型的階躍型滑坡[7-9],至今仍處于間歇性蠕滑變形狀態(tài)。眾多學(xué)者對八字門滑坡開展了位移預(yù)測研究,周超等提出了一種基于誘發(fā)因素響應(yīng)分析的進化支持向量機位移預(yù)測模型,以當(dāng)月降雨量、前兩月累積降雨量、當(dāng)月庫水位平均值、當(dāng)月庫水位最大變化量為影響因素開展位移預(yù)測,認為該方法適合于的階躍型滑坡位移預(yù)測[10]。尚敏等[11]認為八字門滑坡累積位移曲線的階躍性與降雨量、庫水位下降速率呈正相關(guān),建立了一元線性回歸模型對滑坡階躍段的累積位移開展預(yù)測,取得了不錯的效果。梁陽等[12]認為八字門滑坡的當(dāng)月降雨量會影響下一月的位移,降雨響應(yīng)的滯后期約為一個月。以上對研究數(shù)據(jù)的選擇通常以月為單位,時間跨度較大,在數(shù)據(jù)的選擇上存在較大的精度誤差,難以給出更為準(zhǔn)確的滯后期。因此,提升數(shù)據(jù)的研究精度以及考慮降雨滯后效應(yīng)的滑坡位移預(yù)測是一個重要課題。
降雨是誘發(fā)滑坡變形破壞的重要因素,從時間的角度分析,滑坡的變形破壞多發(fā)生在降雨的中后期或略滯后[13-14]。田正國等[15]對三峽庫區(qū)滑坡誘發(fā)因素分析認為:降雨會導(dǎo)致滑坡體物質(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ù)值模擬軟件,未從物理模型實驗的角度開展驗證分析。
以上學(xué)者在八字門滑坡位移預(yù)測、滑坡與降雨滯后效應(yīng)兩方面開展了研究。但以天為時間單位考慮降雨滯后效應(yīng)的位移預(yù)測研究較少,本文將基于八字門滑坡日降雨量和自動位移監(jiān)測數(shù)據(jù)分析降雨滯后期,采用灰色神經(jīng)網(wǎng)絡(luò)預(yù)測模型,對八字門滑坡開展位移預(yù)測研究,以期進一步提高八字門滑坡位移預(yù)測的精確度。
八字門滑坡地處三峽庫區(qū)香溪河河口右岸,地理坐標(biāo):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ā)育于斜逆向巖層的堆積層滑坡,滑床基巖為侏羅紀(jì)早期紫紅色砂頁巖、棕色砂巖及燧石礫石,產(chǎn)狀:292∠29°。詳見圖1、圖2所示。八字門滑坡于1981年葛洲壩蓄水后復(fù)活,進入初期復(fù)活階段;受1998年長時間暴雨影響,滑坡變形加劇趨于整體復(fù)活,以牽引式為主;自2003年三峽大壩蓄水以來,滑坡前緣已被庫水淹沒,滑坡變形模式由牽引式轉(zhuǎn)變?yōu)橥埔剖剑?2]。經(jīng)現(xiàn)場調(diào)查可知:新的宏觀變形往往發(fā)生在每年的6~8月份,主要分布于滑坡中部、滑坡后緣等位置,坡體中后緣公路路面已受到嚴(yán)重破壞,影響通行。
圖1 八字門滑坡地質(zhì)圖[11]Fig.1 Geological map of Bazimen landslide
圖2 八字門滑坡剖面圖[21]Fig.2 Cross-section of Bazimen landslide
八字門滑坡于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
八字門滑坡的階躍變形對降雨的響應(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]。
在分布滯后模型參數(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代表當(dāng)天降雨量與當(dāng)天位移的相關(guān)系數(shù);s=1表示前第1 d的降雨量與當(dāng)天位移的相關(guān)系數(shù);s=2表示前第2 d的降雨量與當(dāng)天位移的相關(guān)系數(shù),依次類推。分析可知:s=0~12時,相關(guān)系數(shù)先增大后減小,呈現(xiàn)中等強度至弱相關(guān);s=13時,相關(guān)性系數(shù)小于0.2,呈極弱相關(guān)或無相關(guān),即滑坡當(dāng)天位移與前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ī)律。
將八字門滑坡2019年08月18日至2020年08月31日的實測日降雨量數(shù)據(jù)代入(式2),對2019年09月01日至2020年08月31日的日變形位移進行估算,并與位移實際值對比驗證該分布滯后模型的準(zhǔn)確性,結(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)測值的變化趨勢和突變點的位置一致,估算值曲線前移(如圖中黑色方框標(biāo)示),說明在一定程度上抵消了降雨與位移之間的滯后效應(yīng)。該模型反映的降雨響應(yīng)滑坡位移變化的滯后關(guān)系是較為準(zhǔn)確的,但估算值與監(jiān)測值在峰值處仍有較大偏差,這是因為降雨分布滯后模型只考慮了降雨因素,實際上滑坡變形還受到其他因素的影響,故需要進一步提高位移預(yù)測的精度。
八字門滑坡日變形量受到前期12天降雨的影響,將降雨滯后效應(yīng)考慮到滑坡日位移預(yù)測之中,可以提高滑坡預(yù)測精度,在滑坡位移預(yù)測中考慮滑坡的誘發(fā)因素也可以提高預(yù)測精度[26-27],周超等基于滑坡誘發(fā)因素對八字門滑坡位移預(yù)測取得了很好的預(yù)測效果[10],文章選取滯后期內(nèi)12 d累積降雨量、庫水位漲落速率、滑坡自身位移作為滑坡變形的影響因素,對八字門滑坡的日變形位移進行預(yù)測研究。
八字門滑坡的自動位移只有最近4年的數(shù)據(jù),數(shù)據(jù)樣本量較少,滑坡變形受到多種因素影響呈非線性特征,灰色模型具備小樣本數(shù)據(jù)建模優(yōu)點[28],神經(jīng)網(wǎng)絡(luò)適用于非線性和不確定性系統(tǒng)預(yù)測,灰色神經(jīng)網(wǎng)絡(luò)組合模型結(jié)合了灰色系統(tǒng)和神經(jīng)網(wǎng)絡(luò)的優(yōu)點,能夠?qū)Χ喑梢虻膹?fù)雜未知系統(tǒng)進行建模,因此選擇灰色神經(jīng)網(wǎng)絡(luò)模型對八字門滑坡日變形位移的預(yù)測研究。
灰色模型是指對原始數(shù)據(jù)序列x(0)一次累加后得到新序列x(1),新序列規(guī)律性會明顯增強,然后對該新序列進行擬合的模型[29]。
設(shè)時間數(shù)據(jù)序列x(0):
累加之后的新序列x(1)為:
根據(jù)新的序列x()1,建立白化微分方程,即:
該方程的解為:
是x t(1)的估計值,對其做一次累減得到x(0)的預(yù)測值x t*(0):
灰色神經(jīng)網(wǎng)絡(luò)模型是在灰色理論基礎(chǔ)上構(gòu)建的,為了方便表述,重新定義符號,用x(t)表示x(0),y(t)代表一次累加后生成的數(shù)列x t(1),z(t)表示預(yù)測結(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ò)預(yù)測值。
選擇八字門滑坡GSCX3全自動監(jiān)測點2016年5月24日-2020年6月12日數(shù)據(jù)(日變形位移、庫水位漲落速率、12 d累積降雨量)作為訓(xùn)練樣本,對2020年6月24日~2020年8月28日的位移進行預(yù)測,借助Matlab編寫程序處理數(shù)據(jù),訓(xùn)練樣本原始數(shù)據(jù)見表4。
表4 訓(xùn)練樣本原始數(shù)據(jù)表Table 4 Raw data table of training samples
首先,對輸入原始數(shù)據(jù)(訓(xùn)練樣本)進行歸一化處理,消除數(shù)據(jù)之間數(shù)量級的差異,歸一化函數(shù)形式如下:
式中:xmin是x k序列的最小值;xmax是序列的最大值。
其次,構(gòu)建灰色神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu),位移預(yù)測輸入數(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ù)為預(yù)測位移。
然后,設(shè)置網(wǎng)絡(luò)訓(xùn)練參數(shù),對灰色神經(jīng)網(wǎng)絡(luò)進行訓(xùn)練,使網(wǎng)絡(luò)具有滑坡日變形位移預(yù)測的能力,在訓(xùn)練過程中網(wǎng)絡(luò)的連接權(quán)值也將不斷地被調(diào)整修正。
LD層中神經(jīng)元閾值為:
最后,用訓(xùn)練好的灰色神經(jīng)網(wǎng)絡(luò)預(yù)測結(jié)果輸出,將輸出層進行一次累減并反歸一化得到預(yù)測值,預(yù)測結(jié)果分析詳見圖6和表5。
表5 八字門滑坡部分日變形位移預(yù)測分析表Table 5 Prediction table of partial daily deformation and displacement of Bazimen landslide
圖6 八字門滑坡日位移實際值與預(yù)測值對比分析圖Fig.6 Comparison and analysis diagram of actual and predicted daily displacement of Bazimen landslide
由圖6可知,考慮降雨滯后效應(yīng)、庫水位漲落速率、滑坡自身位移作為滑坡位移影響因素的灰色神經(jīng)網(wǎng)絡(luò)預(yù)測模型,能夠準(zhǔn)確地預(yù)測滑坡日位移變化趨勢和突變點,計算可知平均絕對百分比誤差(MAPE)為11.98%,精確度88.02%,預(yù)測值與實際監(jiān)測值的最大偏差小于2 mm,符合模型預(yù)測的精度要求??蔀榘俗珠T滑坡的監(jiān)測和防災(zāi)減災(zāi)工作提供理論支撐和參考建議,同時可為三峽庫區(qū)同類型滑坡的預(yù)測研究提供借鑒。
(1)八字門滑坡的累積位移“階躍”發(fā)生在每年的降雨集中時期,八字門滑坡的變形主要受強降雨或持續(xù)性降雨影響。
(2)采用日降雨量和日變形位移數(shù)據(jù)建立了八字門滑坡位移分布滯后模型,模型估計值與滑坡實際監(jiān)測值變化趨勢一致,且估計值曲線前移,在一定程度上滯后效應(yīng)已被消除,八字門滑坡日變形位移響應(yīng)降雨的滯后期為12 d的分布滯后模型是較為準(zhǔn)確的。
(3)在灰色神經(jīng)網(wǎng)絡(luò)預(yù)測模型考慮降雨滯后效應(yīng),利用滯后期內(nèi)12 d的累積降雨量、庫水位漲落速率、滑坡自身位移作為影響因子,對八字門滑坡日變形位移預(yù)測,能夠準(zhǔn)確地預(yù)測位移的變化趨勢及突變點。平均絕對百分比誤差11.98%,精確度88.02%,最大偏差2 mm。可為八字門滑坡的預(yù)警預(yù)報和防治提供可靠建議。