馮 諭,曾懷恩,2,3,涂鵬飛,3

(1.三峽大學 土木與建筑學院,湖北 宜昌 443002;2.三峽大學 湖北長江三峽滑坡國家野外科學觀測研究站,湖北 宜昌 443002; 3.三峽大學 湖北省水電工程施工與管理重點實驗室,湖北 宜昌 443002)

0 引 言

三峽庫區(qū)一直是我國地質(zhì)災害多發(fā)區(qū),滑坡災害的高發(fā)性與高危性使得其成為三峽庫區(qū)最具破壞力的自然災害,嚴重的滑坡災害對人民生活造成了重大影響?；碌乇硭轿灰祁A測一直是眾多專家學者研究防治滑坡災害的有效手段。對于“階躍式”滑坡位移的預測,目前大多數(shù)學者采用的是經(jīng)時序分解后的機器學習預測方法。因此,合理地對滑坡位移進行時序分解是決定模型預測精確性、解釋性、運用合理性的重要手段。

目前已有的時序分解方法主要有:2022年陳曦[1]提出的一種滑坡每期位移變化量與位移變化均值作比較,剔除階躍項,并應用多項式擬合的方法;二次移動平均提取趨勢項并應用多項式擬合方法(2022年曹博等[2]應用于八字門滑坡位移預測);三次移動平均提取趨勢項并應用多項式擬合方法(2015年周超等[3]應用于八字門滑坡位移預測);二次指數(shù)平滑法提取趨勢項并應用多項式擬合方法(2018年黃發(fā)明等[4]應用于白水河滑坡位移預測);變分模態(tài)分解(Variational Mode Decomposition,VMD)趨勢分解法(2019年趙力學等[5]用于河流流量預測);Cox-Stuart檢驗趨勢分解法(2020年胡成雨等[6]用于產(chǎn)品質(zhì)量檢驗);多項式最小二乘分解法等。

三峽庫區(qū)白水河、八字門等典型滑坡都屬于“階躍式”滑坡。上述方法下的滑坡位移預測雖然在短期內(nèi)相對精確,但其分解方法未考慮到滑坡最重要的階躍項位移。這致使在之后的滑坡位移預測中,可能因波動項位移包含有階躍項位移而造成預測不精確,同時,多項式擬合與線性擬合無法對滑坡位移的趨勢作出合理的力學解釋。因此,應對現(xiàn)有滑坡位移數(shù)據(jù)進行階躍項(突變點)檢測,并選擇滑坡蠕滑力學本構(gòu)模型作為時序分解趨勢項的合理擬合模型?，F(xiàn)有的突變點檢測方法主要有:Mann-Kendall突變檢測方法(2022年邵駿等[7]應用于金沙江上游河段水溫變化突變檢測研究);滑動t檢驗突變檢測方法(2021年歐陽衛(wèi)等[8]應用于灃河徑流量突變分析);啟發(fā)式分割(Bernaola Galvan, BG)算法突變檢驗方法(2005年封國林等[9]應用于氣候突變檢測研究);pettitt突變檢測方法(2014年傅麗昕[10]應用于豐縣氣溫與降水量趨勢性與突變特征檢測);Yamamoto突變檢測方法(2007年張建軍等[11]應用于合肥年平均氣溫與年平均降雨突變檢測分析);Bayes突變檢測方法(2007年隋學深等[12]應用于滬深兩市股指時序突變點研究)等。這些突變檢測方法對于具有單調(diào)變化趨勢的滑坡位移來說,適用性較差。

鑒于此,本文根據(jù)西原蠕變本構(gòu)模型與自適應改進遺傳算法模型,提出一種基于滑動t檢驗的滑動Rnl(曲線擬合優(yōu)度)階躍檢測方法及改進移動平均修正與提取階躍項位移方法,并將其應用于白水河滑坡水平位移時序分解。

1 理論和方法

1.1 滑坡位移時序分解

建立新型滑坡位移時序加法模型[13],即

F(t)=c(t)+w(t)+m(t) 。

(1)

式中:t為時間變量;F(t)為滑坡體變形水平位移;c(t)為以滑坡自重應力為主的地質(zhì)應力作用引起的蠕滑趨勢項位移;m(t)為強降雨、庫水位劇烈變動等地質(zhì)因素引起的階躍項位移;w(t)為其他因素引起的波動項位移。

c(t)為滑坡蠕滑位移本構(gòu)模型關(guān)系式;w(t)為周期項誤差波動關(guān)系式;m(t)表達式為

m(t)=

(2)

式中:x為未產(chǎn)生階躍的累積天數(shù);y為產(chǎn)生階躍的累積天數(shù);r(t)為降雨量;v(t)為庫水位變動速率;rq與vq分別為降雨量與庫水位的階躍閾值;f為以r(t)及v(t)為變量的階躍函數(shù)。

本文認為r(t)或v(t)達到階躍閾值(rq,vq)時,產(chǎn)生階躍項位移fi(r(t),v(t));階躍項位移具有累積效應,且不可逆,故有

x+y=n。

(3)

式中n為分析總期數(shù)。

1.2 滑坡蠕滑位移本構(gòu)模型

文獻[14]—文獻[16]驗證了西原模型對于滑坡位移-時間模型趨勢項擬合的適用性。因此,選取傳統(tǒng)西原蠕變模型為滑坡趨勢項位移c(t)的擬合模型。

據(jù)文獻[17]、文獻[18],傳統(tǒng)的西原蠕變模型方程為

(4)

式中:E0、E1分別為彈性模量與黏彈性模量;σ為總應力;σs為滑坡體屈服應力;e(t)為滑坡蠕變;η1、η2為黏滯系數(shù)。

據(jù)式(4),得到2條不同應力條件下的應變-時間(ε-t)曲線,如圖1所示。

圖1 不同應力條件下的蠕變曲線

考慮到滑坡位移處于不斷累增狀態(tài),取滑坡位移-時間模型趨勢項表達為變換參數(shù)后的西原模型關(guān)系式②,即

c(t)=a1(1-e-b1t)+c1t+D。

(5)

式中a1、b1、c1、D為趨勢項位移待求解力學參數(shù)。

1.3 滑動Rnl檢測

滑動Rnl檢測是基于滑動t檢驗與滑坡蠕滑位移本構(gòu)模型下的一種針對滑坡位移的新型突變檢測方法。具體原理如下。

Rnl是基于直線回歸擬合優(yōu)度決定系數(shù)R2的曲線回歸改進系數(shù),作用在于判斷曲線擬合的擬合優(yōu)度。Rnl定義[19]為

根據(jù)最小二乘擬合規(guī)則,對滑坡水平位移5組以上數(shù)據(jù)進行逐點擬合。按照圖1曲線②模式,控制c(t)參數(shù)的取值范圍,保留大范圍階躍后線性的存在可能。

根據(jù)Rnl定義式,階躍式位移產(chǎn)生之前,由于波動項的存在,Rnl表現(xiàn)為小范圍波動;階躍項位移產(chǎn)生時,突變點產(chǎn)生,最小二乘下的擬合函數(shù)在控制范圍內(nèi)無法收斂至突變點,使得Rnl有明顯的下降趨勢;階躍項位移消失時,由于突變消失,滑坡位移繼續(xù)向原趨勢發(fā)展,Rnl有一定的上升趨勢,并最終變?yōu)椴▌訝顟B(tài)。

1.4 自適應改進遺傳算法模型

曲線回歸最小二乘擬合大多數(shù)采用的是Newton迭代法,即選取參數(shù)的初值,運用Newton迭代式將其迭代至理想精度范圍內(nèi),完成求解。因此,對大范圍參數(shù)初值進行搜尋,進而得到最優(yōu)擬合曲線是保證滑動Rnl檢測順利進行的必要操作。本文選用自適應改進遺傳算法模型進行Rnl最優(yōu)值的求解。

遺傳算法(Genetic Algorithm,GA)是一種智能算法,該算法把需要求解的優(yōu)化問題模擬為一個生物自然進化過程,經(jīng)過不斷地迭代更新,最終找到一個最優(yōu)解[20-21]。遺傳算法可以直接進行搜尋求解計算,沒有求導或連續(xù)性的界定,其運算的隨機搜索性使其具有良好的全局尋優(yōu)能力[22]。

本文運用遺傳算法尋優(yōu)的次數(shù)較多,尋優(yōu)難度也有區(qū)別。因此,為了提高遺傳算法效率,建立自適應改進遺傳算法模型,如下:

設置自適應改進遺傳算法選擇模型,即

(7)

式中:p為被選擇的個體;pfitness_max為種群中的最大適應度;pfitness_ j為當前個體適應度;j為當前累計個體數(shù);i為當前累計數(shù)據(jù)個數(shù);q為Rnl開始有大幅度下降趨勢時的當前累計數(shù)據(jù)個數(shù)。

設置自適應改進遺傳算法交叉與變異算子,即

(8)

式中:NPC、NMU分別為交叉算子與變異算子;NPC_cloud、NMU_cloud分別為云模式(cloud)下的自適應交叉算子與變異算子。

(1)產(chǎn)生大階躍變形前,滑坡位移按蠕滑趨勢進行發(fā)展,非線性擬合較為明顯,遺傳算法尋優(yōu)困難。因此,設定選擇模型為輪盤賭選擇模型[23]并依據(jù)文獻[24]改進輪盤賭適應度模型;同時,設定云模式[25]下的自適應交叉算子與變異算子,設定方式依據(jù)文獻[26]。

(2)產(chǎn)生大階躍變形后,滑坡位移趨勢狀態(tài)被打破,擬合曲線可能會呈線性或者類線性形態(tài),遺傳算法尋優(yōu)較為簡單。因此,設定選擇模型為精英選擇模型;同時,設置交叉、變異模型為傳統(tǒng)單點交叉、變異模型,并設置交叉算子與變異算子為常用固定值0.9與0.1。

1.5 改進移動平均法剔除階躍項

滑坡位移呈遞增趨勢發(fā)展,因此,階躍項位移直接進行剔除會導致后續(xù)數(shù)據(jù)無法處理。提出一種基于移動平均原理的階躍項位移剔除方法,具體如下:

考慮到加權(quán)移動平均法有觀察期的近期觀察值對預測值有較大影響的特征,應用加權(quán)移動平均法改進式進行階躍項位移的剔除。

傳統(tǒng)加權(quán)移動平均法適用于平穩(wěn)時序,對于非線性遞增時序無法應用,故做適當改進,即

(9)

wL-i=L-i。

(10)

依據(jù)加權(quán)移動平均原理,保留接近階躍項的前5個權(quán)重,令

wL-6=wL-7=…=wL-n=0 。

(11)

綜上,式(9)可簡化為

(12)

階躍點剔除并重新設置后,階躍項位移為

m(L)=F(L)-mL。

(13)

式中:m(L)為L處剔除的階躍項位移;F(L)為L處的滑坡總位移。

重置階躍項后續(xù)位移數(shù)據(jù),即

F(L+e)1=F(L+e)-m(L) 。

(14)

式中:F(L+e)1為L+e處的重置位移;F(L+e)為L+e處的原滑坡總位移;e=1,2,…,u,u為下個識別階躍點的前一項。

2 白水河滑坡實例分析

2.1 白水河滑坡地質(zhì)概況

白水河滑坡位于湖北省宜昌市三峽庫區(qū)秭歸縣,該滑坡體前緣高程70 m,后緣高程390 m,南北長780 m,東西寬700 m,總體積1 260萬m3,坡體地層巖性為侏羅紀中厚層狀砂巖夾薄層狀泥巖,產(chǎn)狀為15°∠36°[27]。本文以ZG118監(jiān)測點2003年6月—2008年8月的滑坡累計水平位移監(jiān)測數(shù)據(jù)作為實例進行時序分解。ZG118監(jiān)測點累計水平位移曲線如圖2所示。

圖2 ZG118監(jiān)測點水平位移

據(jù)圖2,ZG118累計水平位移曲線表現(xiàn)為明顯的階躍遞增形態(tài)。

2.2 滑坡位移階躍點識別

滑坡位移-時間趨勢項表達見式(5),根據(jù)西原蠕變表達式,對式(5)中待求解參數(shù)進行初值搜索范圍的約束。各參數(shù)搜索取值范圍為

(15)

根據(jù)自適應改進遺傳算法模型,得到滑動Rnl檢測方法下的Rnl變化曲線,如圖3。

圖3 滑動Rnl檢測

據(jù)圖3,將Rnl降低的點作為階躍點,得到滑坡位移突變點分布,如圖4所示。

圖4 滑坡位移突變點分布

根據(jù)突變點分布,由于突變幅度較大,產(chǎn)生了部分識別錯誤點,對識別錯誤點51、52、53進行剔除。剔除的突變點Rnl降低,為保證擬合性能,仍對剔除突變點進行改進移動平均修正。

2.3 滑坡位移時序分解

根據(jù)改進移動平均法剔除識別的階躍點,修正階躍數(shù)據(jù)與后續(xù)位移數(shù)據(jù)。對剔除階躍點并修正后的位移數(shù)據(jù)做西原蠕變趨勢項擬合,將擬合曲線作為趨勢項位移曲線,擬合誤差曲線作為波動位移曲線,剔除的階躍曲線作為階躍項位移曲線。得到時序分解后的滑坡位移表達式如下:

據(jù)總位移表達式(式(1)),趨勢項位移、階躍項位移、波動項位移表達式分別為:

c(t)=127.4(1-e-1.095t)+4.685t-4.025; (16)

(17)

w(t)=F(t)-m(t)-c(t) 。

(18)

滑坡位移時序分解后的趨勢項位移曲線、波動項位移曲線、階躍項位移曲線如圖5所示。

圖5 時序分解后各位移曲線

滑坡總位移時序分解后的分解效果如圖6所示。

圖6 滑坡總位移分解效果

2.4 不同階躍點檢測方法對比分析

經(jīng)典突變檢測方法主要有:MK(Mann-Kendall)突變檢測方法、滑動t檢驗突變檢測方法、BG算法突變檢驗方法、pettitt突變檢測方法以及Bayes突變檢測方法等。其中,MK突變檢測方法、滑動t檢驗突變檢測方法、pettitt突變檢測方法以及Bayes突變檢測方法最為常用。

由于MK突變檢測方法與pettitt突變檢測方法檢測原理類似,因此,本文選用MK突變檢測方法、滑動t檢驗突變檢測方法、Bayes突變檢測方法與滑動Rnl檢測方法做對比分析。突變檢測結(jié)果如圖7所示。

圖7 突變檢測結(jié)果

根據(jù)圖7,MK突變檢測結(jié)果表明,滑坡位移有很強的單調(diào)遞增趨勢,只檢驗出其具有趨勢性,無法檢測到突變點;滑動t檢驗檢測出的突變點存在有明顯且較多平穩(wěn)序列,檢測局部正確;Bayes突變檢測結(jié)果表明,滑坡位移突變點產(chǎn)生概率均相等,無法識別特殊突變點。相較于上述方法,對于滑坡位移而言,滑動Rnl突變檢測方法更加適用。

2.5 不同時序分解方法對比分析

現(xiàn)階段較為常用的時序分解方法有移動平均法,指數(shù)平滑法、VMD方法、多項式最小二乘分解法等。由于多項式最小二乘分解較為簡單,且數(shù)據(jù)量有限,本文選用二次移動平均法、三次指數(shù)平滑法、VMD分解方法與滑動Rnl檢測時序分解方法做對比。

選取前5個數(shù)據(jù)作為一次移動平均計算數(shù)據(jù),一次移動平均前5個數(shù)據(jù)作為二次移動平均計算數(shù)據(jù);為保證指數(shù)平滑分解后曲線平滑性,指數(shù)平滑系數(shù)α定為0.3。VMD分解時不對分解后的模態(tài)個數(shù)k進行討論,采用傳統(tǒng)分解方法,分解出的波動曲線集合作為波動項位移曲線。白水河滑坡二次移動平均分解、VMD時序分解、三次指數(shù)平滑分解以及新型時序分解結(jié)果如圖8所示。

圖8 時序分解結(jié)果

根據(jù)圖8可知:

(1)二次移動平均法時序分解后丟失了一部分原始數(shù)據(jù),對后續(xù)的時序預測會造成一定的影響;分解出的趨勢項位移無規(guī)律性,也沒有力學理論支撐,可能會造成由強降雨、庫水位急劇變動或者偶然因素引起的非趨勢項位移也計算在趨勢項位移內(nèi);分解后的波動項位移變化幅度較大,預測的誤差可能會超出允許誤差范圍;由于處理后趨勢項數(shù)據(jù)分布仍不規(guī)則,三次多項式擬合的效果也較差,均方根誤差RMSE為59.07 mm,且擬合曲線發(fā)展趨勢與原始數(shù)據(jù)發(fā)展趨勢也不相同。二次移動平均趨勢項擬合曲線如圖9所示。

圖9 二次移動平均趨勢項擬合曲線

(2)三次指數(shù)平滑法時序分解后同樣有無規(guī)律性,將非趨勢項位移計算在趨勢項位移內(nèi)的可能;分解后的波動項位移變化幅度仍然較大;三次多項式擬合的效果相較二次移動平均法好,但均方根誤差RMSE仍有34.14 mm,且擬合曲線發(fā)展趨勢與原始數(shù)據(jù)發(fā)展趨勢不同的問題仍然存在。三次指數(shù)平滑趨勢項擬合曲線如圖10所示。

圖10 三次指數(shù)平滑趨勢項擬合曲線

(3)相較于前2種時序分解方法,VMD時序分解也存在趨勢項無規(guī)律性等問題。同時,趨勢項位移前期出現(xiàn)了大量負值,與實際現(xiàn)象不符。此方法并不適用。

綜上并結(jié)合圖8(d)可知,滑動Rnl檢測時序分解方法有如下優(yōu)勢:①解決了趨勢項位移無規(guī)律、無力學依據(jù)的問題; ②分解后的趨勢項位移無負值;③波動項位移變化幅度較小,預測精度有所提高;④西原蠕變模型擬合效果較好,均方根誤差RMSE為18.56 mm,發(fā)展趨勢具有遞增性且符合力學規(guī)律;⑤引入“階躍式”滑坡位移預測中最重要的階躍項位移,分析預測更具有針對性。

3 結(jié) 論

本文根據(jù)西原蠕變本構(gòu)模型與自適應改進遺傳算法模型,提出了滑動Rnl階躍項位移檢測方法與改進加權(quán)移動平均修正階躍項位移方法。并將其應用于白水河滑坡位移時序分解,結(jié)果表明:

新型時序分解方法優(yōu)勢:①解決了以往時序數(shù)學分解方法上的無規(guī)律、無力學理論依據(jù)的問題;②對于“階躍式”滑坡,引入代表其特征的階躍項位移,分析預測更具針對性;③波動項位移變化幅度減小,預測精度有所提高。

新型時序分解方法存在的問題:雖然分解后的趨勢項與波動項位移的預測精度均會提高,且較其他分解方法更具有合理性,但是分解出的階躍項位移如何進行分析預測是較大的難點,也是作者未來的研究方向。

文中提出的新型時序分解方法的應用案例僅有白水河滑坡一例,新方法是否可以廣泛應用于“階躍式”滑坡位移的分解預測還有待于更多實例驗證。