基于改進蟻群算法的邊坡穩定分析方法研究

馬 俊，資 強，呂錄娜，陳麗剛

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

在巖土工程領域，邊坡穩定性分析是一個極為重要而古老的研究課題，占有十分重要的地位。長期以來，國內外眾多學者們在這方面做了大量的研究工作，得出了很多重要的研究成果。目前，極限平衡法和有限元分析方法[1- 2]是進行邊坡穩定分析的主要方法，其中有限元法包括有限元強度折減法和滑面應力分析法。極限平衡法假定土體為剛性體，不能考慮土體的真實應力應變關系，同時對條間力的一些假定使得坡體受力分析不夠準確；盡管有限元強度折減法克服了極限平衡法的這些缺點，但假定邊坡處于某一極限平衡狀態進行受力分析得出邊坡的安全系數，且邊坡失穩的判斷標準沒有達成統一的認識。

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

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

1 改進蟻群算法

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

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

對于連續空間優化問題，引入二進制對優化變量進行編碼。假定二進制串的長度為k位，則基于圖解的蟻群系統可以編制成如圖1所示的構造圖[9]。

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

圖1 改進蟻群算法解構造圖

在圖1中，任意2個總結點Ni和Nj中間有4條并行路徑，分別用來連接子結點Ni0→Nj0，Ni0→Nj1，Ni1→Nj0，Ni1→Nj1，其相應路徑上的信息素分別記為τ(i，j)1，τ(i，j)2，τ(i，j)3，τ(i，j)4。每只螞蟻都會根據τ(0，1)s，(s∈{1，2，3，4})的大小概率選擇前兩個子結點。如果螞蟻選擇了子結點Ni0，就需要根據τ(i，j)1和τ(i，j)2的大小按照一定狀態轉移規則選擇路徑Ni0→Nj0或者Ni0→Nj1；同理，如果選擇了子結點Ni1，就需要根據τ(i，j)3和τ(i，j)4的大小按照一定狀態轉移規則選擇路徑Ni1→Nj0或者Ni1→Nj1。依次按概率選擇下一個子節點，當螞蟻走完終止節點時，就建立了一條完整的路徑。

2 邊坡穩定分析數學模型的建立

2.1 邊坡安全系數的定義

在傳統的極限平衡法中，滑裂面的安全系數被定義為抗滑力與滑動力的比值，其表達式為：

(1)

式中，τf—土體的抗剪強度；τ—土的實際剪應力。

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

已知滑裂面上某一點的應力狀態為σx、σy、τxy，則滑裂面上的正應力σn和剪應力τn的計算公式如下：

(2)

(3)

式中，α—滑裂面上該點處的正應力與x正方向的夾角。

該點的土體抗剪強度為τf，其計算公式為：

τf=c+σntanφ

(4)

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

于是，整個滑裂面上的安全系數可定義為：

(5)

2.2 選取設計變量

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

圖2 邊坡滑裂面模型圖

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

(6)

(7)

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

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

(8)

在傳統各旅行商問題或者路徑規劃問題中，最優路徑對應的是最短距離，而在邊坡穩定分析問題中，最優路徑為其對應目標函數值即邊坡滑動安全系數最小的螞蟻行進的路徑。

2.3 確定約束條件

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

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

(2)滑裂面的滑入點必須在滑出點的左側，即xo

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

2.4 目標函數

目標函數取臨界滑裂面對應的最小安全系數Fs，用數學描述為：

FS=minFS(X)，X∈Ω

(9)

3 應用研究

3.1 計算模型和計算參數

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

圖3 邊坡模型幾何尺寸

圖4 邊坡有限元計算網格

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

3.2 計算成果分析

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

圖5 邊坡水平向應力分布圖

圖6 邊坡豎向應力分布圖

圖7 邊坡剪應力分布圖

圖8 邊坡等效塑性應變場

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

圖9 最小安全系數收斂過程

圖10 安全系數平均值收斂過程

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

本文蟻群優化算法及文獻[10]法搜索到的最危險滑裂面如圖11—12所示，最小安全系數見表1。

表1 計算結果對比

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

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

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

4 結語

本文在邊坡土體應力場的基礎上定義了潛在滑裂面的安全系數，綜合考慮了土體的真實受力狀態和應力應變關系，構建了基于改進蟻群算法的邊坡穩定分析新方法，實現了在連續解空間范圍內搜尋邊坡的最小安全系數和滑裂面。經驗證，該方法評價邊坡穩定性是可行的，計算結果是可靠的。本文提出的邊坡穩定分析方法尚未初步研究成果，通過簡單算例對其進行了驗證，下一步將對更加復雜的邊坡工程的應用以及非圓弧滑裂面的分析做深入的研究工作。