曹 博， 汪 帥， 宋丹青， 杜 涵， 劉光偉， 周志偉

(1.遼寧工程技術(shù)大學(xué) 礦業(yè)學(xué)院， 遼寧 阜新 123000； 2.清華大學(xué) 水沙科學(xué)與水利水電工程國家重點(diǎn)試驗(yàn)室，北京 100084； 3.神華寶日希勒能源有限公司露天煤礦， 內(nèi)蒙古 呼倫貝爾 021000)

1 研究背景

我國是世界上滑坡災(zāi)害發(fā)育最為嚴(yán)重的國家之一[1-2]，它的頻繁發(fā)生和快速運(yùn)動(dòng)，每年都會(huì)造成大量的人員傷亡和財(cái)產(chǎn)損失[3]。2018年的《全國地質(zhì)災(zāi)害通報(bào)》表明，2017年全年我國共發(fā)生地質(zhì)災(zāi)害2 966起，其中滑坡災(zāi)害1 631起，占比為55.0%，因此，如何有效防止滑坡災(zāi)害成為了該研究領(lǐng)域的熱點(diǎn)問題。

近年來，滑坡位移預(yù)測(cè)已逐漸成為滑坡早期預(yù)警的有效手段，它允許研究人員通過室內(nèi)計(jì)算來分析滑坡歷史變形的主要特征，并對(duì)未來一段時(shí)間內(nèi)滑坡的位移趨勢(shì)進(jìn)行預(yù)測(cè)[4]。由于該方法所需成本較低并且精度較高，因此得到了快速發(fā)展。目前常用的滑坡位移定量預(yù)測(cè)方法主要包括數(shù)理統(tǒng)計(jì)模型和非線性理論[5]，2009年Li等[6]提出在滑坡預(yù)測(cè)過程中使用尖點(diǎn)突變模型，王珣等[7]在2017年建立切線角模型進(jìn)行了臨滑預(yù)報(bào)判據(jù)研究，但是這些方法通常需要對(duì)滑坡的演化狀態(tài)進(jìn)行較為深入準(zhǔn)確的分析，不確定性較大。與這些模型相比，機(jī)器學(xué)習(xí)模型不需要分析復(fù)雜的滑坡機(jī)理且對(duì)于以滑坡為代表的非線性動(dòng)態(tài)系統(tǒng)具有較好的擬合效果，因此應(yīng)用更加廣泛，最具代表性的模型包括2019年柳青等[8]進(jìn)行滑坡位移預(yù)測(cè)預(yù)報(bào)所采用的神經(jīng)網(wǎng)絡(luò)模型、Zhou等[9]和李仕波等[10]先后于2016及2019年預(yù)測(cè)位移時(shí)采用的支持向量機(jī)模型、2015年周超等[11]預(yù)測(cè)滑坡位移時(shí)引入的混沌理論。極限學(xué)習(xí)機(jī)模型(extreme learning machine, ELM)作為近年來應(yīng)用廣泛的智能算法，其在滑坡位移預(yù)測(cè)方面的適用性得到了不少研究的驗(yàn)證，如Cao等[12]利用ELM模型對(duì)白家包滑坡的位移進(jìn)行了預(yù)測(cè)，同時(shí)與支持向量機(jī)模型進(jìn)行了預(yù)測(cè)效果比較；Huang等[13]對(duì)滑坡位移利用混沌理論和小波變換進(jìn)行分解后，再基于ELM模型進(jìn)行了非平穩(wěn)時(shí)間序列的預(yù)測(cè)。在ELM模型的計(jì)算過程中，由于有個(gè)別參數(shù)是隨機(jī)選取的初始值，因此為加快收斂速度并提高準(zhǔn)確性，需要加入一定的算法進(jìn)行選參優(yōu)化，如粒子群算法[10]、遺傳算法[5]、蟻群算法[14]等常見方法。然而，對(duì)于多個(gè)算法優(yōu)化效果的對(duì)比，前人的研究較為欠缺。

綜上所述，為驗(yàn)證ELM模型在我國三峽庫區(qū)滑坡位移預(yù)測(cè)中的適用性，本文選擇湖北秭歸八字門滑坡為研究對(duì)象，在深入分析其變形歷史和影響因素后，完成了基于蟻群算法優(yōu)化ELM模型的位移預(yù)測(cè)過程，并通過對(duì)比驗(yàn)證了不同優(yōu)化算法的優(yōu)化效果，以期為后續(xù)研究提供科學(xué)依據(jù)。

2 蟻群算法優(yōu)化極限學(xué)習(xí)機(jī)模型

2.1 基于移動(dòng)平均法的時(shí)間序列分解

基于時(shí)間序列分析原理，對(duì)于通過監(jiān)測(cè)或遙感等手段得到的滑坡位移數(shù)據(jù)，可以將其看成是一種非線性時(shí)間序列，并且有如下關(guān)系[15]：

yt=St+Ct+εt

(1)

式中：yt為滑坡累計(jì)位移，mm；St、Ct和εt分別為滑體的趨勢(shì)項(xiàng)、周期項(xiàng)和隨機(jī)項(xiàng)位移，mm。其中，趨勢(shì)項(xiàng)位移代表了滑體在長時(shí)間段內(nèi)的變形趨勢(shì)，通常為單調(diào)遞增；周期項(xiàng)位移表示滑坡位移在周期性因素影響下的變形；隨機(jī)項(xiàng)位移為滑坡在隨機(jī)因素影響下產(chǎn)生的位移，由于影響因素較難量化，因此以往研究通常將該項(xiàng)忽略。

對(duì)于滑坡位移預(yù)測(cè)的研究，通常要利用一定的方法將滑坡的累計(jì)位移分解成不同的分量，通過對(duì)子序列進(jìn)行分析預(yù)測(cè)，疊加所有子序列的預(yù)測(cè)值來得到總的滑坡位移預(yù)測(cè)值。分解方法主要包括兩種，即數(shù)理統(tǒng)計(jì)法和智能方法，后者又包括小波分解[11]、經(jīng)驗(yàn)?zāi)B(tài)分解[16]、集合經(jīng)驗(yàn)?zāi)B(tài)分解[17]等方法，但是這類方法總體來說較為復(fù)雜。針對(duì)滑坡位移時(shí)序的特點(diǎn)，本文選擇移動(dòng)平均法進(jìn)行總位移的分解，其計(jì)算公式如下：

Mi={M1,M2,…,Mt,…,Mn}

(2)

(3)

2.2 極限學(xué)習(xí)機(jī)

2006年Huang等[18]首先提出了極限學(xué)習(xí)機(jī)，其本質(zhì)是一種單層神經(jīng)網(wǎng)絡(luò)，但相比較而言，ELM模型的結(jié)構(gòu)更為簡單，且運(yùn)算速度更快[19]，并且其中的層間連接權(quán)值及隱含層偏差不需要事先設(shè)置。在ELM算法中，隱含層神經(jīng)元閾值及輸入層與隱含層的連接權(quán)值由函數(shù)生成，訓(xùn)練過程中無需人為干預(yù)，只需進(jìn)行隱含層神經(jīng)元個(gè)數(shù)設(shè)置即可獲得該條件下的最優(yōu)解。

若設(shè)模型中輸入層與隱含層之間連接權(quán)值為w，同時(shí)隱含層與輸出層間的連接權(quán)值為β，給定訓(xùn)練樣本總數(shù)為Q，ELM神經(jīng)網(wǎng)絡(luò)的輸出為T，激活函數(shù)為g(x)時(shí)，隱含層輸出矩陣H可由下式求解：

Hβ=T′

(4)

(5)

式中：H+為輸出矩陣H的廣義逆矩陣。

極限學(xué)習(xí)機(jī)模型在理論上可以獲得最小的訓(xùn)練誤差，利用ELM模型進(jìn)行滑坡位移預(yù)測(cè)模型構(gòu)建的過程中，所需考慮的參數(shù)僅為隱含層神經(jīng)元個(gè)數(shù)，利用ELM模型進(jìn)行位移預(yù)測(cè)主要包含以下步驟：

(1)對(duì)隱含層相關(guān)參數(shù)隨機(jī)進(jìn)行初始化，確定隱含層神經(jīng)元個(gè)數(shù)，并隨機(jī)生成輸入層與隱含層間的連接權(quán)值及隱含層神經(jīng)元偏差的初始值；

(2)合理選擇一個(gè)無限可微函數(shù)作為隱含層激活函數(shù)，以便求解隱含層輸出矩陣H；

2.3 蟻群算法優(yōu)化極限學(xué)習(xí)機(jī)

Parpinelli等[20]于2002年首次提出了蟻群算法，蟻群算法類似于自然界中螞蟻的覓食過程，螞蟻群體會(huì)通過不同路徑搜尋食物，對(duì)于正確的路徑，螞蟻的數(shù)量和密度會(huì)越來越大，而對(duì)于其他路徑則相反。最后螞蟻能夠搜尋到一條抵達(dá)食物源頭的最佳路徑，即問題的最優(yōu)解[20-21]。蟻群算法最初用于旅行商難題上，隨著科研人員的不斷探索，應(yīng)用范疇幾乎涉及各個(gè)領(lǐng)域。為更好更快地找到最優(yōu)的連接權(quán)值與神經(jīng)元閾值，本文擬采用蟻群優(yōu)化算法對(duì)極限學(xué)習(xí)機(jī)模型進(jìn)行改進(jìn)，以達(dá)到更優(yōu)的計(jì)算效果。

使用蟻群算法優(yōu)化極限學(xué)習(xí)機(jī)模型的流程如下所示：

(1)初始化參數(shù)，計(jì)算螞蟻從一個(gè)點(diǎn)轉(zhuǎn)移到另一個(gè)點(diǎn)的概率，概率計(jì)算公式如下：

(6)

(2)通過不停迭代，尋求一個(gè)最佳的覓食路徑，確定蟻群算法全局最優(yōu)路徑，即極限學(xué)習(xí)機(jī)的權(quán)值和閾值；

(3)基于蟻群算法最優(yōu)路徑的計(jì)算，合理選擇隱含層激活函數(shù)，求解隱含層輸出矩陣H；

(4)對(duì)輸出層權(quán)值矩陣進(jìn)行求解，以此為基礎(chǔ)進(jìn)行滑坡位移預(yù)測(cè)。

3 研究實(shí)例

3.1 八字門滑坡概況

八字門滑坡位于湖北省宜昌市秭歸縣歸州鎮(zhèn)香溪河右岸，距離三峽大壩38 km。滑坡高程分布在140～270 m之間，滑體平均長度和平均寬度分別為350和400 m，厚度約為30 m，體積達(dá)400×104m3，屬于大型堆積層滑坡。根據(jù)野外鉆孔揭露的情況顯示，滑坡包括兩層滑帶，上部為次級(jí)滑面，厚度約為3 m，下部的主滑帶連續(xù)分布，中部的厚度最大?；禄鶐r為侏羅系香溪組(J1x)的砂巖和泥質(zhì)粉砂巖，呈強(qiáng)分化特性，局部夾雜有炭質(zhì)泥巖。八字門滑坡平面及剖面圖如圖1所示。

3.2 滑坡變形特征及影響因素分析

自從八字門滑坡發(fā)生較為明顯的變形開始，地質(zhì)災(zāi)害管理部門就對(duì)其變形進(jìn)行了持續(xù)的監(jiān)測(cè)，主要包括深部位移和地表GPS位移。在滑坡的監(jiān)測(cè)剖面Ⅰ-Ⅰ′上(圖1)，安裝了3個(gè)地表監(jiān)測(cè)點(diǎn)(ZG110、ZG111、ZG112)和1個(gè)深部位移監(jiān)測(cè)點(diǎn)(QTZK3)，監(jiān)測(cè)點(diǎn)多年監(jiān)測(cè)結(jié)果以及相應(yīng)的庫水位和月降雨量變化見圖2，2005年3月與2004年9月相比滑坡深部監(jiān)測(cè)點(diǎn)的相對(duì)位移如3所示。

3個(gè)地表監(jiān)測(cè)點(diǎn)的高程分別約為165、191和215 m，可以代表滑坡前、中、后部的變形特點(diǎn)和趨勢(shì)。圖2中滑坡累計(jì)位移曲線顯示，更靠近滑坡后緣的ZG111監(jiān)測(cè)點(diǎn)的位移量最大，總位移超過了1 700 mm，因?yàn)樵擖c(diǎn)所在位置是滑坡的陡坎處，在每年的雨季和庫水位下降期(5-9月)，該點(diǎn)都會(huì)發(fā)生較為明顯的階躍式變形增長，如從2009年的5月開始，在5個(gè)月的時(shí)間內(nèi)，該點(diǎn)的累計(jì)位移增長超過了200 mm，增長率達(dá)到了24.9%。但是位于滑坡中部的ZG110監(jiān)測(cè)點(diǎn)和靠近滑坡前緣的ZG112監(jiān)測(cè)點(diǎn)的總位移卻較小，分別約為1 200 mm和360 mm，且該兩個(gè)監(jiān)測(cè)點(diǎn)的位移增長趨勢(shì)較為平穩(wěn)。上述結(jié)果說明八字門滑坡的運(yùn)動(dòng)模式為“前進(jìn)式”，即最開始由滑坡上部發(fā)生蠕變并向下擠壓，隨著位移的逐漸增加，變形開始向滑坡前緣發(fā)展，最終形成整體破壞。

圖1 香溪河八字門滑坡平面及剖面圖

圖2 2004-2012年八字門滑坡各監(jiān)測(cè)點(diǎn)地表位移監(jiān)測(cè)曲線以及相應(yīng)庫水位和月降雨量變化

圖3 八字門滑坡2005年3月相對(duì)于2004年9月深部監(jiān)測(cè)點(diǎn)的位移

選取變形量最大的ZG111監(jiān)測(cè)點(diǎn)對(duì)滑坡影響因素和變形間的響應(yīng)關(guān)系進(jìn)行分析，在三峽水庫蓄水前，滑坡前緣被庫水淹沒，加之地下水上升帶來的侵蝕作用，滑動(dòng)面前部的抗剪強(qiáng)度逐漸降低，引起了滑坡前部的緩慢蠕變，并在前緣逐漸形成了少量張拉裂縫，因此在這段時(shí)間內(nèi)滑坡存在漸進(jìn)式變形的特征。2003年三峽水庫開始蓄水，在八字門滑坡上出現(xiàn)了較多的變形跡象；同年8月，滑坡體上道路兩側(cè)出現(xiàn)了兩個(gè)較大的裂縫，長度分別達(dá)到了30和45 m；2004年的5月至8月，庫水位短暫下降，在道路上又發(fā)現(xiàn)了多處裂縫，如圖4所示。

圖4 八字門滑坡體上公路變形實(shí)拍圖

由于水庫的周期蓄水，滑坡的累計(jì)位移也隨之發(fā)生變化。在2007年2-4月，水庫水位由156 m下降至152 m，但累計(jì)位移幾乎不變，但是在該年的5-9月，由于總降水量達(dá)到近600 mm，該點(diǎn)的累計(jì)位移迅速增加，變化量超過了120 mm，同時(shí)伴隨有多處地表裂縫的出現(xiàn)。這是因?yàn)槎虝r(shí)間內(nèi)的強(qiáng)降水導(dǎo)致地下水位的抬升，而該滑坡主要由堆積物構(gòu)成，高孔隙率和良好的滲透性使得地表水的入滲更加容易，導(dǎo)致了異常的孔隙水壓力和水動(dòng)力壓力誘發(fā)的滑坡變形。綜上所述，滑坡的變形主要受到了庫水位波動(dòng)和降雨入滲的影響。

3.3 滑坡位移分解及預(yù)測(cè)

3.3.1 基于移動(dòng)平均法的累計(jì)位移分解 由于ZG111測(cè)點(diǎn)處的累計(jì)位移最大，而且可以反映滑坡的階躍式變形特點(diǎn)，因此選取該點(diǎn)進(jìn)行位移預(yù)測(cè)?？傆?jì)獲取102組數(shù)據(jù)，將前90組數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)集，后12組數(shù)據(jù)作為測(cè)試數(shù)據(jù)集，首先根據(jù)公式(1)、(2)將2005年7月至2013年1月的全部數(shù)據(jù)分為趨勢(shì)項(xiàng)位移和周期項(xiàng)位移，如圖5所示。由圖5可看出，提取出的趨勢(shì)項(xiàng)位移較為平滑，基本不存在突變現(xiàn)象，說明移動(dòng)平均法可以很好地消除滑坡階躍式曲線的影響，從而反映出真正的滑坡長期變形特征。另一方面，滑坡的周期項(xiàng)位移也具有較明顯的周期特征，而且初步估計(jì)該曲線的周期近似為12個(gè)月，也符合實(shí)際情況。

圖5 2005年7月至2013年1月八字門滑坡ZG111測(cè)點(diǎn)累計(jì)位移分解結(jié)果

3.3.2 滑坡位移分量預(yù)測(cè) 趨勢(shì)項(xiàng)位移反映滑坡長時(shí)間以來的變形趨勢(shì)，在目前的研究中該位移多使用多項(xiàng)式函數(shù)進(jìn)行預(yù)測(cè)。對(duì)于上節(jié)中分解得到的趨勢(shì)項(xiàng)數(shù)據(jù)，利用最小二乘法進(jìn)行預(yù)測(cè)。選取圖5中前79組數(shù)據(jù)先進(jìn)行擬合，采用二次多項(xiàng)式作為擬合函數(shù)，擬合之后再將公式外推獲得后面測(cè)試數(shù)據(jù)集的預(yù)測(cè)結(jié)果，利用本案例的訓(xùn)練數(shù)據(jù)最后獲得的預(yù)測(cè)公式如下：

y=0.0955t2+6.0632t+324.59 (R2=0.998)

(6)

趨勢(shì)項(xiàng)位移預(yù)測(cè)結(jié)果公式(6)的擬合優(yōu)度達(dá)到了0.998，表明該公式的預(yù)測(cè)效果很好，能夠有效反映滑坡變形隨時(shí)間變化的整體趨勢(shì)。

對(duì)于周期項(xiàng)位移數(shù)據(jù)，它通常與滑坡的周期性影響因素有關(guān)，本文選取降雨和庫水位變化作為滑坡周期性的影響因素，再參考類似研究，還需要考慮滑坡變形對(duì)兩種誘發(fā)因素的響應(yīng)滯后時(shí)間。因此，最終選擇上月降雨量、當(dāng)月降雨量、上月庫水位平均高程和當(dāng)月庫水位平均高程作為影響因素。這些影響因素被作為輸入數(shù)據(jù)集，而當(dāng)月周期項(xiàng)位移被作為輸出數(shù)據(jù)集，然后基于自編代碼構(gòu)建蟻群算法優(yōu)化的極限學(xué)習(xí)機(jī)模型，進(jìn)行周期項(xiàng)位移的預(yù)測(cè)。以2012年2月至2013年1月為預(yù)測(cè)期，得出的滑坡趨勢(shì)項(xiàng)位移和周期項(xiàng)位移預(yù)測(cè)結(jié)果與實(shí)際值比較如圖6所示。由圖6可看出，周期項(xiàng)位移預(yù)測(cè)結(jié)果與實(shí)際較為吻合，表明預(yù)測(cè)效果也較為理想。

圖6 2012年2月至2013年1月滑坡趨勢(shì)項(xiàng)和周期項(xiàng)位移預(yù)測(cè)結(jié)果與實(shí)際值比較

4 討 論

將3.3節(jié)得到的趨勢(shì)項(xiàng)和周期項(xiàng)兩項(xiàng)分量相加，得到預(yù)測(cè)期蟻群算法優(yōu)化的極限學(xué)習(xí)機(jī)模型預(yù)測(cè)的最終累計(jì)位移結(jié)果，如圖7所示，同時(shí)，將單一的極限學(xué)習(xí)機(jī)模型預(yù)測(cè)結(jié)果以及采用兩種常用的優(yōu)化算法(粒子群算法和遺傳算法)對(duì)極限學(xué)習(xí)機(jī)模型進(jìn)行優(yōu)化后的預(yù)測(cè)結(jié)果一并繪入圖7，并與位移實(shí)際值進(jìn)行比較，以分析本文所采用的蟻群算法優(yōu)化的極限學(xué)習(xí)機(jī)模型的預(yù)測(cè)效果。上述各模型預(yù)測(cè)效果主要參數(shù)對(duì)比見表1。

圖7 各模型預(yù)測(cè)結(jié)果

分析圖7和表1可知：(1)在所有4種模型當(dāng)中，單一極限學(xué)習(xí)機(jī)模型的計(jì)算用時(shí)最短，但是其計(jì)算精度最差，平均誤差達(dá)到了23.5 mm，擬合優(yōu)度為0.973，此2項(xiàng)指標(biāo)均低于其他3種預(yù)測(cè)模型。然而，考慮到滑坡的累計(jì)位移已經(jīng)超過了1 000 mm，因此單一極限學(xué)習(xí)機(jī)模型的預(yù)測(cè)精度也是可以接受的，這說明該模型在階躍式滑坡位移曲線方面的預(yù)測(cè)是較為可靠的，也驗(yàn)證了極限學(xué)習(xí)機(jī)模型對(duì)于本文應(yīng)用實(shí)例的適用性。(2)對(duì)于優(yōu)化后的極限學(xué)習(xí)機(jī)模型，蟻群算法優(yōu)化的極限學(xué)習(xí)機(jī)計(jì)算用時(shí)為98 s，位移預(yù)測(cè)結(jié)果平均誤差為10.1 mm，預(yù)測(cè)結(jié)果擬合優(yōu)度為0.998，各指標(biāo)均優(yōu)于其他兩種算法優(yōu)化的極限學(xué)習(xí)機(jī)模型。因此，對(duì)于本文所示的案例而言，蟻群算法具有最佳的參數(shù)選取與優(yōu)化效果，蟻群算法優(yōu)化的極限學(xué)習(xí)機(jī)模型是進(jìn)行三峽庫區(qū)八字門滑坡累計(jì)位移預(yù)測(cè)的最優(yōu)模型，可在類似滑坡體的位移預(yù)測(cè)中進(jìn)一步研究與推廣。

表1 各模型預(yù)測(cè)效果主要指標(biāo)對(duì)比

5 結(jié)論與展望

本文探討了蟻群算法優(yōu)化的極限學(xué)習(xí)機(jī)模型在三峽庫區(qū)堆積層滑坡位移預(yù)測(cè)中的應(yīng)用，并以秭歸縣八字門滑坡為例進(jìn)行了預(yù)測(cè)效果分析，主要得到以下幾點(diǎn)結(jié)論：

(1)極限學(xué)習(xí)機(jī)模型能夠較好地進(jìn)行三峽庫區(qū)八字門滑坡的累計(jì)位移預(yù)測(cè)，單一極限學(xué)習(xí)機(jī)模型預(yù)測(cè)的平均誤差為23.5 mm，擬合優(yōu)度可達(dá)0.973，是一種可靠的位移預(yù)測(cè)方法。

(2)蟻群算法優(yōu)化的極限學(xué)習(xí)機(jī)模型比單一的極限學(xué)習(xí)機(jī)模型的預(yù)測(cè)精度更高，說明該算法的優(yōu)化性能較好；同時(shí)與其他2種傳統(tǒng)的優(yōu)化算法(粒子群算法、遺傳算法)相比，蟻群算法在計(jì)算用時(shí)和計(jì)算精度上更優(yōu)，其預(yù)測(cè)結(jié)果與實(shí)際結(jié)果間的平均誤差為10.1 mm，擬合優(yōu)度達(dá)0.998。該模型可在類似滑坡的位移預(yù)測(cè)研究中進(jìn)行推廣。

(3)在蟻群算法優(yōu)化參數(shù)選取的過程中，部分算法參數(shù)仍需進(jìn)行進(jìn)一步討論，算法參數(shù)值與計(jì)算結(jié)果的用時(shí)和精度間的敏感性分析將會(huì)是未來的一個(gè)研究方向。