基于改進(jìn)蟻群算法的邊坡穩(wěn)定分析方法研究

馬 俊，資 強(qiáng)，呂錄娜，陳麗剛

(1.黃河勘測規(guī)劃設(shè)計(jì)有限公司，河南 鄭州 450003；2.中國市政工程華北設(shè)計(jì)研究總院有限公司，天津 300074)

在巖土工程領(lǐng)域，邊坡穩(wěn)定性分析是一個(gè)極為重要而古老的研究課題，占有十分重要的地位。長期以來，國內(nèi)外眾多學(xué)者們?cè)谶@方面做了大量的研究工作，得出了很多重要的研究成果。目前，極限平衡法和有限元分析方法[1- 2]是進(jìn)行邊坡穩(wěn)定分析的主要方法，其中有限元法包括有限元強(qiáng)度折減法和滑面應(yīng)力分析法。極限平衡法假定土體為剛性體，不能考慮土體的真實(shí)應(yīng)力應(yīng)變關(guān)系，同時(shí)對(duì)條間力的一些假定使得坡體受力分析不夠準(zhǔn)確；盡管有限元強(qiáng)度折減法克服了極限平衡法的這些缺點(diǎn)，但假定邊坡處于某一極限平衡狀態(tài)進(jìn)行受力分析得出邊坡的安全系數(shù)，且邊坡失穩(wěn)的判斷標(biāo)準(zhǔn)沒有達(dá)成統(tǒng)一的認(rèn)識(shí)。

邊坡穩(wěn)定分析的最終目的是尋找最小安全系數(shù)，鑒于巖土邊坡土層組成和邊界條件的復(fù)雜性，其安全系數(shù)是空間和時(shí)間的復(fù)雜函數(shù)，很難給出顯式計(jì)算公式。蟻群算法作為一種啟發(fā)式人工智能型優(yōu)化算法，具有魯棒性好，對(duì)目標(biāo)函數(shù)(安全系數(shù))無特殊要求，逐漸被學(xué)者們接受，并在邊坡穩(wěn)定分析問題中得到初步的應(yīng)用。陳昌富等[3]通過在蟻群算法中引入混沌擾動(dòng)算子，提出了具有混沌擾動(dòng)算子的啟發(fā)式蟻群算法，計(jì)算結(jié)果表明該方法可以有效的搜索出邊坡的最小安全系數(shù)和相應(yīng)的滑裂面。王成華等[4- 5]利用邊坡體的非線性有限元數(shù)值分析結(jié)果定義安全系數(shù)，分別研究了邊坡臨界滑裂面的蟻群算法和遺傳算法搜索技術(shù)，經(jīng)實(shí)例驗(yàn)證，其提出的新算法可以用來分析邊坡穩(wěn)定性。石露等[6]在邊坡應(yīng)力場的基礎(chǔ)上，提出基于邊坡應(yīng)力場的矢量和法求解滑裂面的安全系數(shù)，在露天采場邊坡工程中聯(lián)合蟻群算法與遺傳算法搜尋邊坡的臨界滑裂面和最小安全系數(shù)。

綜上，盡管很多學(xué)者將蟻群算法成功的應(yīng)用到了邊坡穩(wěn)定分析中，但在計(jì)算滑裂面安全系數(shù)時(shí)，采用條分法，不能充分考慮坡體的真實(shí)受力狀態(tài)；同時(shí)，在搜索滑裂面時(shí)，其搜索思想還是基于基本蟻群算法的有向圖原理，其計(jì)算精度受邊坡離散點(diǎn)的多少制約，而且螞蟻只能經(jīng)過預(yù)先設(shè)定的點(diǎn)，存在滑裂面多樣性不足的缺點(diǎn)。本文在利用非線性有限元應(yīng)力場定義滑裂面安全系數(shù)的基礎(chǔ)上，提出了基于改進(jìn)蟻群算法的邊坡穩(wěn)定性分析新方法，以期解決傳統(tǒng)條分法計(jì)算安全系數(shù)和離散蟻群算法搜索滑裂面的不足。

1 改進(jìn)蟻群算法

受蟻群覓食行為的啟發(fā)，1991年，意大利學(xué)者M(jìn).Dorigo[7]首次在他的博士論文中系統(tǒng)給出了基于蟻群的新型優(yōu)化算法——蟻群算法(Ant Colony Algorithm，簡稱ACA)，成功解決了許多組合優(yōu)化問題。隨后，很多學(xué)者對(duì)蟻群算法做了大量的研究和改進(jìn)，用于解決實(shí)際工程問題。如求解旅行商問題(TSP)、二次分配問題(QAP)、車輛路線問題(VRP)等離散域路徑優(yōu)化難題，取得了不錯(cuò)的效果。

由于蟻群算法是建立在離散域基礎(chǔ)上的，因此其能夠很好的解決復(fù)雜的組合優(yōu)化問題，但是不能有效地處理連續(xù)對(duì)象的優(yōu)化問題。通常在求解連續(xù)對(duì)象優(yōu)化問題時(shí)，需要對(duì)標(biāo)準(zhǔn)的蟻群算法做一定的改進(jìn)。奧地利學(xué)者，W.J.Gutjahr教授(1999年)給出了求解組合優(yōu)化問題的基于構(gòu)造圖概念的通用框架。通過行程對(duì)待優(yōu)化問題的可行解進(jìn)行編碼，并把蟻群啟發(fā)算法應(yīng)用到構(gòu)造的框架里，稱之為基于圖解的蟻群系統(tǒng)[8]。與基本蟻群算法相比較，其對(duì)連續(xù)對(duì)象的優(yōu)化問題具有更強(qiáng)的魯棒性。

對(duì)于連續(xù)空間優(yōu)化問題，引入二進(jìn)制對(duì)優(yōu)化變量進(jìn)行編碼。假定二進(jìn)制串的長度為k位，則基于圖解的蟻群系統(tǒng)可以編制成如圖1所示的構(gòu)造圖[9]。

該構(gòu)造圖共有k+1個(gè)總結(jié)點(diǎn)Ni(i=1，2，…，k+1)，每個(gè)總結(jié)點(diǎn)又分成兩個(gè)子結(jié)點(diǎn)Ni0和Ni1，在編程中分別表示二進(jìn)制的0和1。根據(jù)蟻群算法中螞蟻的特性知，螞蟻在行進(jìn)過程中只能選擇每個(gè)總結(jié)點(diǎn)中的某一個(gè)子結(jié)點(diǎn)。當(dāng)螞蟻從起始節(jié)點(diǎn)N00或N01出發(fā)，依次經(jīng)過結(jié)點(diǎn)N1t，N2t，…，Nkt(t∈{0，1})中某一個(gè)子結(jié)點(diǎn)，這樣最終可構(gòu)成一條完整路徑表示的二進(jìn)制解。

圖1 改進(jìn)蟻群算法解構(gòu)造圖

在圖1中，任意2個(gè)總結(jié)點(diǎn)Ni和Nj中間有4條并行路徑，分別用來連接子結(jié)點(diǎn)Ni0→Nj0，Ni0→Nj1，Ni1→Nj0，Ni1→Nj1，其相應(yīng)路徑上的信息素分別記為τ(i，j)1，τ(i，j)2，τ(i，j)3，τ(i，j)4。每只螞蟻都會(huì)根據(jù)τ(0，1)s，(s∈{1，2，3，4})的大小概率選擇前兩個(gè)子結(jié)點(diǎn)。如果螞蟻選擇了子結(jié)點(diǎn)Ni0，就需要根據(jù)τ(i，j)1和τ(i，j)2的大小按照一定狀態(tài)轉(zhuǎn)移規(guī)則選擇路徑Ni0→Nj0或者Ni0→Nj1；同理，如果選擇了子結(jié)點(diǎn)Ni1，就需要根據(jù)τ(i，j)3和τ(i，j)4的大小按照一定狀態(tài)轉(zhuǎn)移規(guī)則選擇路徑Ni1→Nj0或者Ni1→Nj1。依次按概率選擇下一個(gè)子節(jié)點(diǎn)，當(dāng)螞蟻?zhàn)咄杲K止節(jié)點(diǎn)時(shí)，就建立了一條完整的路徑。

2 邊坡穩(wěn)定分析數(shù)學(xué)模型的建立

2.1 邊坡安全系數(shù)的定義

在傳統(tǒng)的極限平衡法中，滑裂面的安全系數(shù)被定義為抗滑力與滑動(dòng)力的比值，其表達(dá)式為：

(1)

式中，τf—土體的抗剪強(qiáng)度；τ—土的實(shí)際剪應(yīng)力。

為了應(yīng)用方便以及具有可比性，安全系數(shù)的定義與傳統(tǒng)極限平衡法中的安全系數(shù)概念一致。本文采用基于剪應(yīng)力的定義方法定義安全系數(shù)。

已知滑裂面上某一點(diǎn)的應(yīng)力狀態(tài)為σx、σy、τxy，則滑裂面上的正應(yīng)力σn和剪應(yīng)力τn的計(jì)算公式如下：

(2)

(3)

式中，α—滑裂面上該點(diǎn)處的正應(yīng)力與x正方向的夾角。

該點(diǎn)的土體抗剪強(qiáng)度為τf，其計(jì)算公式為：

τf=c+σntanφ

(4)

式中，c、φ—土體材料的黏聚力和內(nèi)摩擦角。

于是，整個(gè)滑裂面上的安全系數(shù)可定義為：

(5)

2.2 選取設(shè)計(jì)變量

對(duì)于圓弧滑裂面，通常選用圓弧的圓心坐標(biāo)(xo，yo)和半徑R(或者滑弧深度R-yo)作為優(yōu)化問題的設(shè)計(jì)變量。在實(shí)際工程中，巖土邊坡的坡面形狀通常較為復(fù)雜，因此很難直接計(jì)算出滑裂面與滑弧入口坡面和滑弧出口處坡面的交點(diǎn)坐標(biāo)(如圖2中A點(diǎn)和B點(diǎn))。所以本文對(duì)滑裂面設(shè)計(jì)變量的選取進(jìn)行改進(jìn)，采用滑裂面的滑入點(diǎn)A、滑出點(diǎn)B以及圓心O的x坐標(biāo)xo，xL，xM為優(yōu)化設(shè)計(jì)變量，這樣可以很容易的確定圓弧的幾何控制參數(shù)。

圖2 邊坡滑裂面模型圖

如圖2所示，在已知滑入點(diǎn)A和滑出點(diǎn)B的σx坐標(biāo)后，就可以根據(jù)邊坡的坡面幾何形狀確定相應(yīng)的yM、yL，然后再根據(jù)圓弧滑裂面的幾何關(guān)系可求得圓心的y坐標(biāo)yo，其計(jì)算公式為：

(6)

(7)

若令設(shè)計(jì)變量為X=[Xo，XL，XL]T，則安全系數(shù)可表示為Fs=Fs(X)，X∈Ω，Ω?R3。其中，Ω為滑裂面的可行域，且是三維歐式空間R3的一個(gè)子集。

確定了設(shè)計(jì)優(yōu)化變量以后，就可以對(duì)其進(jìn)行二進(jìn)制編碼。針對(duì)邊坡穩(wěn)定分析問題中的設(shè)計(jì)變量，其可行解可表示為{x1，x2，…，xm}，其中xi用長度為l的二進(jìn)制串表示為{blbl-1…b2b1}，其中bj∈{0，1}，對(duì)應(yīng)于改進(jìn)蟻群算法結(jié)構(gòu)圖1中的一段路徑(NL+1，t，NL+2，t，…，NL+l，t)。假設(shè)xi變量的取值范圍為[xmin，xmax]，則其對(duì)應(yīng)的十進(jìn)制數(shù)值計(jì)算公式為：

(8)

在傳統(tǒng)各旅行商問題或者路徑規(guī)劃問題中，最優(yōu)路徑對(duì)應(yīng)的是最短距離，而在邊坡穩(wěn)定分析問題中，最優(yōu)路徑為其對(duì)應(yīng)目標(biāo)函數(shù)值即邊坡滑動(dòng)安全系數(shù)最小的螞蟻行進(jìn)的路徑。

2.3 確定約束條件

由于計(jì)算滑裂面的安全系數(shù)時(shí)必須保證滑裂面不能脫離坡體，因此選取的設(shè)計(jì)變量要滿足一定的幾何約束條件，以使其能夠形成合理的、符合實(shí)際的滑裂面。如圖2所示，設(shè)計(jì)變量需要滿足以下約束條件：

(1)滑裂面的滑入點(diǎn)、滑出點(diǎn)位于邊坡表面上(即滑入與坡面的交點(diǎn))。

(2)滑裂面的滑入點(diǎn)必須在滑出點(diǎn)的左側(cè)，即xo

(3)滑裂面最低點(diǎn)高于硬土層(基巖)，則R-|yo-yL|≤H，H為邊坡坡腳到硬土層的垂直距離，如圖2所示。

2.4 目標(biāo)函數(shù)

目標(biāo)函數(shù)取臨界滑裂面對(duì)應(yīng)的最小安全系數(shù)Fs，用數(shù)學(xué)描述為：

FS=minFS(X)，X∈Ω

(9)

3 應(yīng)用研究

3.1 計(jì)算模型和計(jì)算參數(shù)

本算例摘自文獻(xiàn)[10]，如圖3所示，某均質(zhì)土坡，坡腳到坡頂高度為20m，坡面傾角為45(xM，xL，xO)，土體重度為20Fs，土體黏聚力42m，內(nèi)摩擦角為17k。有限元網(wǎng)格剖分情況如圖4所示，采用四邊形單元對(duì)邊進(jìn)行離散，共剖分為1912個(gè)單元，2011個(gè)節(jié)點(diǎn)。

圖3 邊坡模型幾何尺寸

圖4 邊坡有限元計(jì)算網(wǎng)格

邊界條件：在計(jì)算模型的左側(cè)豎向邊界和右側(cè)豎向邊界施加水平約束，基礎(chǔ)底部施加固定約束，限制其水平位移和豎向位移，其它為自由邊界。

3.2 計(jì)算成果分析

本算例中土體本構(gòu)關(guān)系采用莫爾-庫倫理想彈塑性模型，使用ABAQUS進(jìn)行有限元數(shù)值計(jì)算，得出邊坡的應(yīng)力場數(shù)據(jù)如圖5—7所示。邊坡有限元非線性數(shù)值分析給出的等效塑性應(yīng)變場分布如圖8所示。

圖5 邊坡水平向應(yīng)力分布圖

圖6 邊坡豎向應(yīng)力分布圖

圖7 邊坡剪應(yīng)力分布圖

圖8 邊坡等效塑性應(yīng)變場

由以上云圖可見，邊坡的應(yīng)力各分量分布基本符合一般規(guī)律，坡腳處存在應(yīng)力集中現(xiàn)象。在以上有限元應(yīng)力場的基礎(chǔ)上，采用本文構(gòu)建的邊坡穩(wěn)定分析數(shù)學(xué)模型搜索邊坡的臨界滑裂面和相應(yīng)的最小安全系數(shù)。蟻群規(guī)模為50只螞蟻，終止迭代次數(shù)為200代，優(yōu)化過程如圖9—10所示。

圖9 最小安全系數(shù)收斂過程

圖10 安全系數(shù)平均值收斂過程

圖9給出了蟻群優(yōu)化算法的全局最小安全系數(shù)的收斂過程。在第25次迭代時(shí)全局最小安全系數(shù)從1.795降低到了1.277，隨后在第152次迭代時(shí)又逐漸降低為1.276，最后趨于穩(wěn)定，全局最小安全系數(shù)僅減小了0.001。由此可見，蟻群優(yōu)化算法先期主要進(jìn)行全局搜索，后期主要是在當(dāng)前最優(yōu)解附近的可行域內(nèi)進(jìn)行細(xì)致的局部搜索。圖10給出了當(dāng)前迭代步安全系數(shù)的平均值與迭代次數(shù)的變化系，由圖可知，蟻群優(yōu)化算法能夠快速的收斂到全局最優(yōu)解附近。另外，在全局最小安全系數(shù)趨于穩(wěn)定的過程中，安全系數(shù)的平均值仍存在一定的波動(dòng)性，這恰恰體現(xiàn)了蟻群算法能夠不斷搜索新解的能力。

本文蟻群優(yōu)化算法及文獻(xiàn)[10]法搜索到的最危險(xiǎn)滑裂面如圖11—12所示，最小安全系數(shù)見表1。

表1 計(jì)算結(jié)果對(duì)比

圖11 收斂后的試算滑裂面分布

圖12 不同方法搜索的臨界滑裂面

由表1比較可知，本文算法得到的最小安全系數(shù)要比文獻(xiàn)中給出的小一些，說明本文算法比文獻(xiàn)[10]算法的全局搜索能力更強(qiáng)。但是由于都是在用有限元給出的應(yīng)力場進(jìn)行安全系數(shù)的計(jì)算，所以二者得到的全局最優(yōu)值很接近，誤差為0.85%。另外，這兩種方法給出的結(jié)果都比簡化Bishop法給出的大，主要是因?yàn)閭鹘y(tǒng)極限平衡法不能考慮邊坡真實(shí)的應(yīng)力變形關(guān)系以及引入的條塊間力分布的假定等，而有限元法能夠提供真實(shí)的邊坡應(yīng)力場，滿足所有的平衡關(guān)系。由圖12可見，盡管各種方法給出的安全系數(shù)略有差異(基本在5%以內(nèi))，但是計(jì)算得到的臨界滑裂面形狀和位置都非常接近，而且滑裂面的滑出點(diǎn)基本都在圖8中的等效塑性應(yīng)變區(qū)內(nèi)。同時(shí)，圖11給出了本文算法收斂后的一系列臨界滑裂面的分布情況，基本都落在一個(gè)很窄的區(qū)域，構(gòu)成了一個(gè)與最優(yōu)滑動(dòng)面一致的滑裂帶，這與邊坡的實(shí)際破壞情況基本吻合，而且滑裂帶的滑出部分和坡腳的塑性區(qū)一致。

4 結(jié)語

本文在邊坡土體應(yīng)力場的基礎(chǔ)上定義了潛在滑裂面的安全系數(shù)，綜合考慮了土體的真實(shí)受力狀態(tài)和應(yīng)力應(yīng)變關(guān)系，構(gòu)建了基于改進(jìn)蟻群算法的邊坡穩(wěn)定分析新方法，實(shí)現(xiàn)了在連續(xù)解空間范圍內(nèi)搜尋邊坡的最小安全系數(shù)和滑裂面。經(jīng)驗(yàn)證，該方法評(píng)價(jià)邊坡穩(wěn)定性是可行的，計(jì)算結(jié)果是可靠的。本文提出的邊坡穩(wěn)定分析方法尚未初步研究成果，通過簡單算例對(duì)其進(jìn)行了驗(yàn)證，下一步將對(duì)更加復(fù)雜的邊坡工程的應(yīng)用以及非圓弧滑裂面的分析做深入的研究工作。