基于元胞自動機(jī)模型的庫岸邊坡失穩(wěn)模擬方法

黃梓莘， 陳 波

(1.中國電建集團(tuán)中南勘測設(shè)計(jì)研究院有限公司， 湖南 長沙 410014； 2.河海大學(xué) 水利水電學(xué)院， 江蘇 南京 210098)

1 研究背景

庫岸邊坡自身地質(zhì)結(jié)構(gòu)復(fù)雜，且外部荷載條件多變，其穩(wěn)定性分析一直是水利工程治理和邊坡失穩(wěn)災(zāi)害評估的研究重點(diǎn)[1-2]?？紤]到水庫在運(yùn)行過程中庫岸邊坡時(shí)刻受到外部環(huán)境的影響，導(dǎo)致其不可避免地呈現(xiàn)動態(tài)變化的特征，有學(xué)者引入直觀可靠的實(shí)測資料，通過關(guān)注特征點(diǎn)位移變化分析庫岸邊坡的穩(wěn)定性[3-5]。安可君等[6]結(jié)合監(jiān)測資料分析重點(diǎn)岸坡的變形性態(tài)，總結(jié)了岸坡變形破壞的發(fā)展演化規(guī)律以及其對周邊環(huán)境的影響。裴向軍等[7]結(jié)合庫水位變化曲線分析了岸坡的變形監(jiān)測數(shù)據(jù)，發(fā)現(xiàn)了岸坡變形與庫水位之間存在的強(qiáng)相關(guān)性。王國強(qiáng)等[8]基于InSAR技術(shù)獲取大范圍、高精度的庫岸變形監(jiān)測信息，開展了流域庫岸的地質(zhì)災(zāi)害研究。然而，水庫蓄水打破了庫岸邊坡的自然演化過程，導(dǎo)致蓄水后庫岸邊坡呈現(xiàn)出與蓄水前不同的滑坡發(fā)育機(jī)制，僅僅憑借工程地質(zhì)勘查資料或宏觀變形監(jiān)測資料，往往難以準(zhǔn)確判斷其運(yùn)行穩(wěn)定情況和未來發(fā)展趨勢[9-10]。

針對復(fù)雜的庫岸邊坡失穩(wěn)分析問題，有學(xué)者引入數(shù)值模擬方法加以研究。晏鄂川等[11]重點(diǎn)考慮了庫水的動力作用因素，通過FLAC3D軟件構(gòu)建了滑坡的數(shù)值-力學(xué)模型。徐寅等[12]采用離散單元法研究了庫岸邊坡的動態(tài)破壞過程，為邊坡風(fēng)險(xiǎn)評價(jià)提供了全新思路。張坤勇等[13]采用有限元強(qiáng)度折減法對邊坡失穩(wěn)的漸進(jìn)破壞過程進(jìn)行了解釋和模擬。鄭宇豪等[14]考慮庫水位和降雨的影響，利用Geo-studio對滑坡進(jìn)行失穩(wěn)破壞概率計(jì)算。高道路等[15]結(jié)合Linear Cohesion接觸模型研究了含水堆積體的運(yùn)動特性，并采用離散元方法計(jì)算堆積體的失穩(wěn)-滑動-堆積全過程。然而，目前的數(shù)值模擬方法大多將庫岸邊坡設(shè)定為均質(zhì)體，而實(shí)際庫岸高邊坡在巖體風(fēng)化、剝蝕、沉積等外力作用下往往形成上部為堆積體、下部為基巖體的二元結(jié)構(gòu)邊坡[16]，無論單獨(dú)從巖質(zhì)邊坡還是單獨(dú)從土質(zhì)邊坡的角度展開分析，均無法全面準(zhǔn)確地描述庫岸高邊坡的運(yùn)行性狀，同時(shí)庫岸邊坡的物質(zhì)組分具有強(qiáng)烈的隨機(jī)性和不均勻性，導(dǎo)致實(shí)際模擬過程中的巖土力學(xué)參數(shù)取值也存在較大困難[17]。

因此，針對傳統(tǒng)邊坡失穩(wěn)分析方法的巖土力學(xué)參數(shù)取值困難和動態(tài)變化特征考慮不足的弊端，本文提出一種基于元胞自動機(jī)模型的庫岸邊坡失穩(wěn)模擬方法，考慮庫岸邊坡不同區(qū)域的物理力學(xué)性質(zhì)，在變形實(shí)測資料驅(qū)動下動態(tài)辨識庫岸邊坡失穩(wěn)路徑，從而模擬庫岸邊坡發(fā)生失穩(wěn)破壞的演化全過程。

2 研究方法

2.1 元胞自動機(jī)模型

元胞自動機(jī)是在離散動力系統(tǒng)中構(gòu)造微觀個體間的相互作用規(guī)則，以描述局部個體在強(qiáng)烈非線性作用下造成整體系統(tǒng)實(shí)現(xiàn)自組織演化的過程[18]，其模型具有建模思路簡單、模擬功能強(qiáng)大、約束條件寬松等優(yōu)點(diǎn)?，F(xiàn)將元胞自動機(jī)模型的主要構(gòu)成元素介紹如下：

(1)元胞：元胞通常以規(guī)則散落在元胞空間的網(wǎng)格形式表示，是模型中的最基礎(chǔ)的組成部分。

(2)元胞狀態(tài)：元胞狀態(tài)多以布爾型變量進(jìn)行表示，即采用0或1分別代表元胞發(fā)生破壞和元胞維持完好的兩種對立狀態(tài)。

(3)元胞空間：元胞空間是指元胞在空間中分布的網(wǎng)格集合，以二維元胞空間為例，主要的網(wǎng)格劃分形式如圖1所示。

圖1 元胞空間主要網(wǎng)格劃分形式

(4)元胞鄰居：元胞鄰居直接影響了元胞自動機(jī)模型模擬的準(zhǔn)確度和復(fù)雜性，二維空間局部范圍內(nèi)的元胞鄰居模型如圖2所示，圖中黑色網(wǎng)格均代表主元胞，而灰色網(wǎng)格則代表元胞鄰居。

圖2 二維空間局部范圍內(nèi)的元胞鄰居模型

(5)元胞時(shí)間：元胞時(shí)間是一個概念抽象的離散尺度，是一系列無量綱的整數(shù)值，通常假設(shè)當(dāng)前時(shí)刻為t，時(shí)間步長Δt為1，則下一時(shí)刻表示為t+1。

(6)演化規(guī)則：根據(jù)當(dāng)前時(shí)刻元胞狀態(tài)以及鄰近范圍的元胞鄰居狀態(tài)，共同確定下一時(shí)刻元胞狀態(tài)的函數(shù)表達(dá)式f，即元胞自動機(jī)的演化規(guī)則，存在：

S1t+1=f(S1t,SNt+1)

(1)

式中：S1t為主元胞t時(shí)刻的狀態(tài);SNt+1為元胞鄰居t+1時(shí)刻的狀態(tài)。

(7)邊界條件：理論上元胞空間是呈無限延伸的狀態(tài)，但由于目前的計(jì)算模擬條件尚無法有效處理無限網(wǎng)格，因此需要對模型設(shè)置邊界條件。

2.2 邊坡失穩(wěn)模擬的元胞自動機(jī)模型

對于庫岸邊坡而言，在空間尺度上，基本單元之間存在相互作用的關(guān)系，基本單元的狀態(tài)變化會受到鄰近單元狀態(tài)的影響；在時(shí)間尺度上，基本單元的狀態(tài)變化則會受到前一時(shí)刻狀態(tài)的影響，這種時(shí)空局部性與元胞自動機(jī)模型的基本原理相契合，故本文引入元胞自動機(jī)模型模擬庫岸邊坡發(fā)生失穩(wěn)破壞的演化過程?，F(xiàn)結(jié)合圖3所示的邊坡失穩(wěn)元胞自動機(jī)模型流程圖，將方法主要步驟介紹如下：

圖3 邊坡失穩(wěn)元胞自動機(jī)模型流程圖

(1)構(gòu)建元胞空間。選定邊坡典型斷面作為分析對象，將三維邊坡簡化為二維平面，然后根據(jù)邊坡的形態(tài)尺寸，將邊坡離散成規(guī)則的正方形網(wǎng)格，設(shè)定包括鄰居型式、時(shí)間步長、邊界條件、能量傳遞系數(shù)等模型參數(shù)，構(gòu)建用于計(jì)算模擬的元胞空間。其中，元胞i的坐標(biāo)為(xi,yi)，i=1,2,…,M，M為元胞空間內(nèi)的元胞總數(shù)?？紤]到邊坡能量傳遞方向具有隨機(jī)性，采用圖2(b)所示的摩爾型鄰居型式，以元胞i為中心，定義局部范圍內(nèi)的元胞為元胞鄰居j(j=1,2,…,8)。

(2)賦予初始能量。邊坡在形成過程中存儲有一定的彈性應(yīng)變能，因此在邊坡失穩(wěn)的元胞自動機(jī)模型中需要考慮元胞初始能量E(xi,yi, 0)?？紤]到邊坡監(jiān)測信息是邊坡運(yùn)行過程中最直觀的表現(xiàn)形式，其變化規(guī)律可以反映邊坡自身的結(jié)構(gòu)性態(tài)及外部環(huán)境的荷載特性和變化規(guī)律，此處結(jié)合邊坡工程的實(shí)時(shí)監(jiān)測資料對各個邊坡元胞進(jìn)行初始能量賦值，其典型力-變形-能量曲線如圖4所示[19]。

圖4 邊坡元胞典型力-變形-能量關(guān)系曲線

將位移測點(diǎn)映射到邊坡典型斷面，選定同一時(shí)間斷面的監(jiān)測資料對不同位移測點(diǎn)位進(jìn)行賦值。在根據(jù)邊坡斷面輪廓尺寸確定的邊界條件約束下，利用克里格法對輸入的測點(diǎn)監(jiān)測信息進(jìn)行插值處理，形成監(jiān)測信息的插值網(wǎng)格δi(xi,yi)，i= 1, 2, …,M，分別計(jì)算各個元胞應(yīng)變能作為邊坡失穩(wěn)模擬的初始能量，其計(jì)算式為：

(2)

式中：Ui為元胞i的應(yīng)變能;δi為元胞i的變形值;σi為元胞i的應(yīng)力值。

對圖4中曲線三角形OAB進(jìn)行簡化考慮，近似計(jì)算為U(xi,yi, 0)=0.5·Eiδi2，其中Ei為元胞i的楊氏模量，由此形成元胞空間的初始狀態(tài)分布。

(3)設(shè)置破壞閾值。由于邊坡的巖土材料具有各向異性，同時(shí)堆積體區(qū)和巖體區(qū)的材料參數(shù)也存在顯著差異，因此在考慮邊坡物理力學(xué)性質(zhì)的情況下，此處采用兩參數(shù)的Weibull概率分布對邊坡的巖土材料參數(shù)進(jìn)行分區(qū)模擬[20]，設(shè)置各個邊坡元胞的破壞閾值：

(3)

式中：x為邊坡巖土的材料參數(shù)；μi為元胞i的彈模特征值；mi為元胞i的材料均質(zhì)度；λ為破壞系數(shù)。

(4)確定失穩(wěn)起始位置。搜索監(jiān)測信息插值網(wǎng)格確定邊坡位移測值的極大值，將該極值點(diǎn)作為模擬邊坡失穩(wěn)破壞的起始位置，并將起始位置所對應(yīng)的元胞狀態(tài)調(diào)整為0，顯示為黑色網(wǎng)格，代表元胞發(fā)生破壞；其余元胞狀態(tài)仍設(shè)置為1，顯示為白色網(wǎng)格，代表元胞維持完好。

(5)設(shè)定演化規(guī)則。在元胞空間中，所有單位元胞按照一定的演化規(guī)則進(jìn)行狀態(tài)轉(zhuǎn)移，這也是邊坡失穩(wěn)路徑的發(fā)展條件。圖5為邊坡元胞空間演化規(guī)則示意圖。

結(jié)合圖5所示對元胞空間的演化規(guī)則具體介紹如下：

① 能量積累。當(dāng)元胞存儲能量E(xi,yi,t)小于破壞閾值Ecr(xi,yi)時(shí)，元胞維持完好，元胞狀態(tài)為1，此時(shí)在外部荷載作用下，元胞將持續(xù)積累能量，則元胞i在時(shí)刻t的能量E(xi,yi,t)表達(dá)式為：

E(xi,yi,t)=E(xi,yi,0)+e(t)

(4)

式中：E(xi,yi,t)為t時(shí)刻元胞i的能量，kJ；E(xi,yi,0)為元胞i的初始能量，kJ；e(t)為外界環(huán)境的輸入能量，kJ，此處采用e(t)=0.05t·Rnd模擬能量輸入方式，其中Rnd為范圍在[0, 1]的隨機(jī)數(shù)，表征邊坡元胞能量積累的非均勻性質(zhì)。

②能量傳遞。當(dāng)元胞存儲能量E(xi,yi,t)大于破壞閾值Ecr(xi,yi)時(shí)，元胞發(fā)生破壞，狀態(tài)由1調(diào)整成0。根據(jù)能量守恒規(guī)則，破壞元胞i所存儲的能量中有一部分被耗散，另一部分按能量傳遞概率向元胞鄰居j傳遞，此時(shí)元胞鄰居j在下一時(shí)刻的存儲能量演化為：

E(xi,yi,t)=Ecr(xi,yi)→E(xi,yi,t)=0

(5)

ΔE(xi,yi)=(1-α)E(xi,yit,t)

(6)

E(xj,yj,t)→E(xj,yj,t+1)=E(xj,yj,t)+

w(i,j)ΔE(xi,yi)

(7)

式中：E(xj,yj,t+1)為t+1時(shí)刻元胞鄰居j的能量積累，kJ；E(xj,yj,t)為t時(shí)刻元胞鄰居j的能量積累，kJ；α為能量耗散系數(shù)，此處根據(jù)元胞鄰居的破壞個數(shù)將一定能量作為耗散能量處理； ΔE(xi,yi)為破壞元胞i所傳遞的能量，kJ；w(i,j)為元胞i與元胞鄰居j之間的能量傳遞概率，同時(shí)表示元胞鄰居吸收能量的能力，故取值參考各元胞鄰居的破壞閾值占所有元胞鄰居閾值之和的比例，即：

(8)

③發(fā)生破壞。在發(fā)生能量傳遞后，若每個元胞的存儲能量均滿足E(xi,yi,t)≤Ecr(xi,yi)，則元胞維持完好狀態(tài)，并根據(jù)公式(4)再次獲取能量；若某個元胞的存儲能量滿足E(xi,yi,t)>Ecr(xi,yi)，則元胞發(fā)生破壞，根據(jù)公式(5)～(7)發(fā)生能量傳遞。

(6)輸出模擬結(jié)果。重復(fù)第(4)步的演化過程，分別記錄元胞空間在不同時(shí)刻的破壞情況，即可模擬邊坡發(fā)生失穩(wěn)破壞的演化過程，以達(dá)到動態(tài)辨識邊坡失穩(wěn)路徑的目的，當(dāng)元胞空間不再發(fā)生能量傳遞，或元胞破壞數(shù)目最終恒定時(shí)，即代表演化過程停止，判定模型計(jì)算結(jié)束，并輸出最終的邊坡失穩(wěn)區(qū)域分布。同時(shí)，為了直觀展現(xiàn)邊坡元胞空間的破壞情況，此處引入邊坡失穩(wěn)破壞指標(biāo)I(t)=N/M(N為t時(shí)刻發(fā)生失穩(wěn)破壞的元胞個數(shù))，從而可定量計(jì)算不同時(shí)刻破壞元胞面積占整個邊坡元胞空間面積的比例。

圖5 邊坡元胞空間演化規(guī)則示意圖

3 工程實(shí)例模擬與分析

3.1 工程概況

本文所選取的研究對象為某水庫左岸邊坡一古老滑坡堆積區(qū)，該邊坡上覆巖體在河谷不斷下切過程中，受重力作用曾多次出現(xiàn)順層滑動的現(xiàn)象，且目前仍然存在失穩(wěn)的可能性，情況嚴(yán)重時(shí)甚至可能危及整個上部邊坡堆積體的穩(wěn)定性。該水庫岸邊坡滑坡體(I區(qū))的位置及工程地質(zhì)剖面圖分別如圖6(a)、6(b)所示。

3.2 元胞自動機(jī)模型參數(shù)設(shè)置

采用元胞自動機(jī)模型對所選庫岸邊坡進(jìn)行失穩(wěn)路徑辨識和失穩(wěn)破壞模擬，首先對邊坡失穩(wěn)元胞自動機(jī)模型的各項(xiàng)參數(shù)進(jìn)行設(shè)置：

(1)元胞空間：在圖6(b)所示的邊坡工程地址剖面圖的基礎(chǔ)上，根據(jù)邊坡形態(tài)尺寸將邊坡典型斷面離散成單位長度為1 m的正方形網(wǎng)格，所構(gòu)建的元胞空間為900×460的元胞列陣，如圖7所示。

(2)邊界條件：選取斷面輪廓作為定值型邊界，在失穩(wěn)路徑辨識圖中以黑色網(wǎng)格表示。

(3)元胞狀態(tài)：各個邊坡元胞被賦予兩種基本狀態(tài){0,1}，其中0表示元胞發(fā)生破壞，網(wǎng)格以黑色網(wǎng)格表示；1表示元胞維持完好狀態(tài)，網(wǎng)格以白色網(wǎng)格表示。

(4)鄰居型式：采用標(biāo)準(zhǔn)摩爾型定義元胞的鄰居型式，并設(shè)定鄰居半徑為1，即定義與主元胞鄰近的8個元胞為元胞鄰居。

(5)元胞破壞閾值：工程實(shí)例邊坡在河谷不斷下切的過程中暴露出假整合面，上覆巖體在重力作用下發(fā)生多次順層滑動，因此在考慮邊坡地質(zhì)構(gòu)造和演化歷史的情況下，參考《地基與基礎(chǔ)》[21]分別選取堆積體和巖體的巖土材料參數(shù)，并采用兩參數(shù)的Weibull概率分布計(jì)算各邊坡元胞的破壞閾值，如圖8所示。

圖6 工程實(shí)例的庫岸邊坡滑坡體位置及地質(zhì)剖面圖

圖7 邊坡典型斷面的元胞空間示意圖 圖8 邊坡元胞破壞閾值分布圖

(6)元胞初始能量：將滑坡體(Ⅰ區(qū))位移測點(diǎn)投影到圖6(b)所示的邊坡典型斷面上，基于各測點(diǎn)2020年全年的位移監(jiān)測數(shù)據(jù)，利用克里格法對輸入監(jiān)測信息進(jìn)行插值處理，初步形成圖9所示的監(jiān)測信息插值網(wǎng)格，以此為基礎(chǔ)分別計(jì)算各個元胞的應(yīng)變能作為邊坡失穩(wěn)模擬的初始能量，進(jìn)一步形成元胞空間的初始狀態(tài)分布。

圖9 克里格插值結(jié)果及邊坡失穩(wěn)模擬起始位置

(7)失穩(wěn)起始位置：為模擬邊坡發(fā)生失穩(wěn)破壞的變形演化過程，此處通過搜索克里格插值結(jié)果中的極大值，假設(shè)邊坡失穩(wěn)元胞自動機(jī)模型的起始位置，選取的邊坡失穩(wěn)模擬起始位置如圖9所示。

匯總邊坡失穩(wěn)的元胞自動機(jī)模型參數(shù)如表1所示，由此構(gòu)建邊坡失穩(wěn)元胞自動機(jī)模型的元胞空間。

表1 邊坡失穩(wěn)的元胞自動機(jī)模型參數(shù)

3.3 邊坡失穩(wěn)破壞模擬結(jié)果

在設(shè)定模型參數(shù)后，構(gòu)建邊坡失穩(wěn)元胞自動機(jī)模型的元胞空間，按照設(shè)定的演化規(guī)律進(jìn)行元胞空間的能量傳遞，并隨時(shí)間步長增加記錄不同時(shí)刻的元胞破壞狀態(tài)。圖10展示了t在10、50、55、60、70和90時(shí)步時(shí)，采用邊坡失穩(wěn)元胞自動機(jī)模型動態(tài)辨識邊坡失穩(wěn)路徑的模擬結(jié)果，圖11展示了邊坡失穩(wěn)破壞指標(biāo)I(t)隨時(shí)步長增加的判別結(jié)果。

圖10 基于元胞自動機(jī)的邊坡失穩(wěn)路徑辨識

圖11 邊坡失穩(wěn)破壞指標(biāo)判別結(jié)果

結(jié)合圖10和11，分時(shí)步分析邊坡失穩(wěn)破壞的模擬過程：

(1)在t=1時(shí)步至t=50時(shí)步失穩(wěn)初期時(shí)，邊坡元胞主要處于能量積累的時(shí)期，破壞元胞的數(shù)目基本沒有增加；

(2)在t=50時(shí)步至t=55時(shí)步時(shí)，邊坡元胞在經(jīng)歷長時(shí)間的能量積累后，從邊坡失穩(wěn)起始位置開始，突然出現(xiàn)了劇烈的大規(guī)模的失穩(wěn)破壞，破壞速度達(dá)到整個模擬時(shí)段的頂峰，破壞方向自高高程向低高程發(fā)展，此時(shí)的邊坡失穩(wěn)區(qū)域主要集中在滑坡堆積區(qū)的地表位置；

(3)在t=55時(shí)步至t=65時(shí)步時(shí)，隨著時(shí)間步長的增加，幾乎每時(shí)步都存在元胞狀態(tài)的改變，但元胞破壞數(shù)目增長速度明顯變慢，元胞失穩(wěn)區(qū)域的擴(kuò)張速度也逐漸減緩，此時(shí)的邊坡失穩(wěn)主要表現(xiàn)為高高程位置的局部坍滑現(xiàn)象；

(4)在t=70時(shí)步時(shí)，邊坡元胞破壞數(shù)目已經(jīng)基本保持穩(wěn)定，邊坡失穩(wěn)破壞區(qū)域不再發(fā)生擴(kuò)張，此時(shí)的失穩(wěn)區(qū)域主要集中在邊坡沿地層界限分布的地滑堆積區(qū)；

(5) 最終，在t=100時(shí)步時(shí)，元胞空間內(nèi)不再發(fā)生能量傳遞，由此判斷基于元胞自動機(jī)的邊坡失穩(wěn)破壞模擬過程結(jié)束，輸出最終的失穩(wěn)區(qū)域分布，此時(shí)失穩(wěn)破壞指標(biāo)I(t)的計(jì)算結(jié)果為21.95%。

綜上所述，庫岸邊坡失穩(wěn)破壞呈現(xiàn)出漸變到突變的累進(jìn)發(fā)展過程，同時(shí)失穩(wěn)破壞區(qū)域主要集中在沿地層界限分布的地滑堆積區(qū)。

4 結(jié) 論

針對傳統(tǒng)庫岸邊坡失穩(wěn)分析方法的巖土力學(xué)參數(shù)取值困難、動態(tài)變化特征考慮不足的弊端，本文提出了基于元胞自動機(jī)模型的庫岸邊坡失穩(wěn)模擬方法，并結(jié)合工程實(shí)例驗(yàn)證了方法的可行性，現(xiàn)得出主要結(jié)論如下：

(1)基于元胞自動機(jī)模型的庫岸邊坡失穩(wěn)模擬方法，在位移實(shí)測數(shù)據(jù)的驅(qū)動下動態(tài)辨識了邊坡失穩(wěn)的發(fā)展路徑，模擬了庫岸邊坡失穩(wěn)破壞的演化過程，驗(yàn)證了庫岸邊坡失穩(wěn)破壞需要經(jīng)歷漸變到突變的累進(jìn)發(fā)展過程，并得出了邊坡失穩(wěn)區(qū)域主要集中在沿地層界限分布的地滑堆積區(qū)的結(jié)論。

(2)對于特征復(fù)雜、組分混合的庫岸邊坡，本文所提出方法在失穩(wěn)路徑辨識上具有一定的應(yīng)用效果和推廣價(jià)值，但在邊坡失穩(wěn)起始位置的選取上比較困難，未來有必要綜合考慮邊坡實(shí)測資料和深入結(jié)合邊坡風(fēng)險(xiǎn)預(yù)警理論，合理確定邊坡失穩(wěn)起始位置以提高模擬結(jié)果的準(zhǔn)確程度。