余秀玲， 申翃， 鄭萬(wàn)棟

(武漢理工大學(xué) 土木工程與建筑學(xué)院， 湖北 武漢 430070)

邊坡穩(wěn)定性分析是巖土工程中的重要研究課題之一，其研究方法包括：工程地質(zhì)類比法、極限平衡法、極限分析法、數(shù)值分析法等。其中極限平衡法的應(yīng)用最為廣泛，但是其必須依賴研究者的經(jīng)驗(yàn)，試算一系列滑面尋找最小安全系數(shù)及其對(duì)應(yīng)的臨界滑面位置。

近年來(lái)，很多學(xué)者致力于臨界滑面自動(dòng)搜索技術(shù)的研究，提出了各種適用于圓弧滑面或非圓弧滑面的搜索方法，如蒙特卡羅法與有限元結(jié)合、臨界化方法、遺傳進(jìn)化算法、模式搜索法、相遇蟻群算法、神經(jīng)網(wǎng)絡(luò)與多元回歸等。與早期使用的窮舉法、二分法相比，上述方法的精度和收斂性提高了很多，對(duì)邊坡穩(wěn)定研究和實(shí)際工程應(yīng)用起到了積極的作用，但有時(shí)也會(huì)發(fā)生計(jì)算成果并非整體安全系數(shù)極值的情況。該文將果蠅優(yōu)化算法和模擬退火算法相融合，在Matlab中編程實(shí)現(xiàn)基于模擬退火機(jī)制的果蠅優(yōu)化算法，通過(guò)兩個(gè)土坡算例驗(yàn)證，并與其他方法進(jìn)行對(duì)比分析。以期為邊坡臨界滑面搜索提供新途徑。

1 邊坡臨界滑面搜索理論模型

簡(jiǎn)化Bishop法是目前常用的一種滑坡穩(wěn)定性分析方法，被認(rèn)為是圓弧滑面計(jì)算分析的最佳極限平衡法。其力學(xué)模型如圖1所示。假定條間力的方向水平，切向剪切力為零，根據(jù)滑面上土體的靜力平衡和整體的力矩平衡可以得到滑面的安全系數(shù)：

Fs=

(1)

式中：Fs為對(duì)應(yīng)滑面的安全系數(shù)；n為條分?jǐn)?shù)；ci、φi分別為滑動(dòng)體第i條塊土層的黏聚力和內(nèi)摩擦角；Wi為第i條塊自重；bi為第i條塊的寬度；θi為第i條塊中心線與圓弧滑面交點(diǎn)處圓弧切線的水平角。

圖1 簡(jiǎn)化Bishop法力學(xué)模型

若忽略條間力的影響，即為瑞典條分法，安全系數(shù)的表達(dá)式簡(jiǎn)化為：

(2)

式中：L為圓弧滑面的弧長(zhǎng)。其他參數(shù)同式(1)。

式(1)為隱式函數(shù)，求解安全系數(shù)時(shí)需要迭代，式(2)為顯式函數(shù)，可以直接求解。

當(dāng)邊坡的幾何參數(shù)和力學(xué)參數(shù)已知時(shí)，安全系數(shù)是滑弧圓心坐標(biāo)和半徑的函數(shù)，即：

Fs=f(x0,y0,R)

(3)

式中：x0、y0、R分別為圓弧滑面的圓心坐標(biāo)和半徑；f為安全系數(shù)Fs與可行搜索域參數(shù)的函數(shù)。

則式(3)轉(zhuǎn)化為一個(gè)優(yōu)化問(wèn)題：

(4)

式中:x(j)={x(1)=x0,x(2)=y0,x(3)=R}為優(yōu)化變量；[a(j),b(j) ]為第j個(gè)變量的取值區(qū)間。

2 基于模擬退火的果蠅優(yōu)化算法

2.1 基本果蠅優(yōu)化算法(FOA)

基本果蠅優(yōu)化算法(FOA)由臺(tái)灣學(xué)者潘文超基于果蠅的覓食行為而開(kāi)發(fā)，是一種新型的群體智能全局優(yōu)化算法，近年來(lái)已被研究人員廣泛運(yùn)用到諸多領(lǐng)域解決優(yōu)化問(wèn)題。果蠅優(yōu)化算法機(jī)理是基于果蠅種群覓食過(guò)程的模擬，利用果蠅良好的嗅覺(jué)和視覺(jué)，設(shè)計(jì)相應(yīng)的感官搜索環(huán)節(jié)，通過(guò)不斷迭代從而獲得最優(yōu)解。果蠅群體尋找食物的迭代過(guò)程如圖2所示。

果蠅優(yōu)化算法步驟如下：

(1) 果蠅種群相關(guān)參數(shù)的初始設(shè)定：位置X_axis、Y_axis，規(guī)模Sizepop，最大迭代次數(shù)Maxgen。

圖2 果蠅群體覓食迭代過(guò)程

(2) 更新果蠅個(gè)體覓食過(guò)程的隨機(jī)距離(RandomValue)與方位。

(5)

(3) 估量果蠅種群中個(gè)體與原點(diǎn)之間的距離。

(6)

(4) 將味道濃度判定值Si代入味道濃度判定函數(shù)，即適應(yīng)度函數(shù)(Fitness function)，得到果蠅個(gè)體的味道濃度Smell(i)。

Smell(i)=Function(Si)

(7)

(5) 在果蠅種群中找出味道濃度最佳的個(gè)體。

[bestSmellbestindex]=min[Smell(i)]

(8)

(6) 記錄并保存最佳味道濃度值bestSmell及其位置X和Y，此時(shí)其他果蠅飛向該位置。

(7) 重復(fù)步驟(2)～(5)，此次與前一次的最佳味道濃度做比較，若優(yōu)于前一次的味道濃度，并且迭代次數(shù)未達(dá)到Maxgen，則執(zhí)行步驟(6)，否則，結(jié)束算法。

基本果蠅優(yōu)化算法能夠快速高效地搜索出最優(yōu)解。但是，由于在利用適應(yīng)度函數(shù)的尋優(yōu)過(guò)程中極易陷入局部最優(yōu)解，導(dǎo)致果蠅種群停止迭代更新位置，因此得到的最優(yōu)解具有局限性。

2.2 模擬退火算法(SA)

模擬退火算法是根據(jù)物理學(xué)中固體退火原理而建立的一種全局最優(yōu)算法。固體加熱升溫時(shí)內(nèi)部粒子排序紊亂，內(nèi)能逐漸增大，降溫時(shí)，粒子排列趨于有序，內(nèi)能減小。假設(shè)在狀態(tài)A時(shí)，系統(tǒng)內(nèi)能為E(A)，受到外界擾動(dòng)后，變?yōu)闋顟B(tài)C，內(nèi)能變?yōu)镋(C)。基于Metropolis準(zhǔn)則，粒子在溫度T趨于平衡的概率為:

(9)

由式(9)可以看出：退火時(shí)，新?tīng)顟B(tài)使內(nèi)能函數(shù)減小，系統(tǒng)一定接受該狀態(tài)；新?tīng)顟B(tài)使內(nèi)能函數(shù)增加時(shí)，系統(tǒng)以一定概率接受該狀態(tài)，即在優(yōu)化進(jìn)程中，可以獲得最優(yōu)解和一定概率的差解。將內(nèi)能函數(shù)E(X)演變?yōu)槟繕?biāo)函數(shù)，退火溫度T演變?yōu)榭刂茀?shù)，隨著T值的變化，算法執(zhí)行“新解-判斷-舍棄或接受”的進(jìn)化，直至求得問(wèn)題的最優(yōu)解。

模擬退火算法的步驟如下：

(1) 隨機(jī)賦予初始狀態(tài)，設(shè)置初始解x0、初始溫度T0以及最大迭代次數(shù)Maxgen，計(jì)算目標(biāo)函數(shù)初始值E(x0)。

(2) 判斷所賦初始值是否滿足終止條件，若是，該解即為最優(yōu)解，結(jié)束算法；否則，在x0鄰域內(nèi)隨機(jī)選取x*重新賦值，計(jì)算內(nèi)能增量ΔE=E(x*)-E(x0)。

(3) 若ΔE<0，則x=x*；否則，以概率exp(-ΔE/T)作為x*的解。

(4) 進(jìn)入迭代尋優(yōu)過(guò)程，重復(fù)步驟(2)和(3)，直至尋得最優(yōu)解或達(dá)到最大迭代次數(shù)。

2.3 基于模擬退火的果蠅優(yōu)化算法(SA-FOA)

如前所述，基本果蠅優(yōu)化算法(FOA)參數(shù)少，易調(diào)節(jié)，效率高，卻易陷入局部最優(yōu)，而模擬退火算法(SA)可以有效地跳出局部最優(yōu)，因此該文將兩種算法融合，形成基于模擬退火的改進(jìn)果蠅優(yōu)化算法。其步驟為：

(1) 執(zhí)行基本果蠅算法的步驟(1)～(6)。

(2) 模擬退火步驟，確定初始溫度T0

T0=Smellbest/log(5)

(10)

(3) 計(jì)算當(dāng)下溫度狀態(tài)果蠅味道濃度判定值Si的適應(yīng)值：

TF(Si)=

(11)

(4) 計(jì)算隨機(jī)概率:

p=rand()

(12)

(5) 采用輪盤賭法策略，從Si中計(jì)算得到全局最佳的某個(gè)值(bestX，bestY)作為下一代果蠅的搜索中心位置。

(13)

(7) 計(jì)算果蠅個(gè)體的味道濃度Smelli*，若Smelli*

(14)

(8) 進(jìn)行退溫操作：

T=0.5T0

(15)

(9) 進(jìn)入迭代尋優(yōu)，重復(fù)執(zhí)行(3)～(9)，直至尋得最優(yōu)解或達(dá)到最大迭代次數(shù)。

將基于模擬退火機(jī)制的果蠅優(yōu)化算法(SA-FOA)應(yīng)用于邊坡的臨界滑面搜索，其流程如圖3所示。

待優(yōu)化的圓弧滑面圓心坐標(biāo)和半徑的取值區(qū)間通過(guò)地質(zhì)勘查和工程類比相關(guān)資料確定。根據(jù)式(16)、(17)進(jìn)行滑弧圓心坐標(biāo)和半徑的初始化：

(16)

(17)

圖3 基于SA-FOA的邊坡臨界滑面搜索流程圖

3 算例分析

基于該文提出的SA-FOA算法，結(jié)合瑞典條分法/簡(jiǎn)化Bishop法，在Matlab中編寫了相應(yīng)的邊坡臨界滑面搜索程序，并對(duì)兩個(gè)算例進(jìn)行了驗(yàn)證分析。

3.1 算例1：雙層土體邊坡

算例1是文獻(xiàn)[8]中的非均質(zhì)土坡，安全系數(shù)計(jì)算采用瑞典條分法。該文分別采用FOA和SA-FOA算法搜索該土坡的臨界滑面，并與文獻(xiàn)[8]的SA算法進(jìn)行比較，參見(jiàn)表1及圖4。3種算法中，F(xiàn)OA算法安全系數(shù)較大，SA-FOA算法安全系數(shù)最小，即基于模擬退火的果蠅優(yōu)化算法尋優(yōu)能力最強(qiáng)。

表1 不同算法的臨界滑面搜索結(jié)果對(duì)比

圖4 不同算法臨界滑面形狀和位置比較

在Matlab中，F(xiàn)OA、SA-FOA搜索最小安全系數(shù)的收斂過(guò)程如圖5所示，對(duì)比可見(jiàn)，SA-FOA迭代曲線的斜率比FOA更大，即SA-FOA收斂速度更快；SA-FOA迭代200次后的安全系數(shù)穩(wěn)定在0.857處，而FOA迭代200次后的安全系數(shù)還未達(dá)到0.93，因此SA-FOA的收斂穩(wěn)定性更好，搜索能力較FOA顯著提高。

圖5 最小安全系數(shù)收斂進(jìn)程

算例1證明將SA-FOA算法用于邊坡臨界滑面的搜索不僅是可行的，而且更高效、準(zhǔn)確。

3.2 算例2：復(fù)雜的3層土體邊坡

算例2取自澳大利亞計(jì)算機(jī)應(yīng)用協(xié)會(huì)(ACADS)發(fā)布的考題，土層剖面如圖6所示，土層參數(shù)見(jiàn)表2。

圖6 ACADS考核題剖面圖

該邊坡土層厚度不均、材料參數(shù)相差較大，是一個(gè)復(fù)雜邊坡，在最危險(xiǎn)滑面搜索過(guò)程中會(huì)出現(xiàn)多個(gè)極值。SA-FOA算法結(jié)合簡(jiǎn)化Bishop 法，假設(shè)一個(gè)初始滑面，則局部危險(xiǎn)滑面及最危險(xiǎn)滑面的搜索結(jié)果列于表3。SA-FOA算法在初始滑面與最危險(xiǎn)滑面(臨界滑面)相差較大的情況下，能夠跨越眾多局部危險(xiǎn)滑面，搜索到整體最危險(xiǎn)滑面，證明了SA-FOA具有較強(qiáng)的全局搜索能力。

表2 各土層材料參數(shù)

表3 不同圓弧滑面的安全系數(shù)

針對(duì)該考題各國(guó)學(xué)者以及裁判推薦答案同時(shí)列于表4，各個(gè)答案均基于簡(jiǎn)化Bishop法。由表4可見(jiàn)，該文SA-FOA計(jì)算結(jié)果與評(píng)判員答案吻合度很高，且臨界滑面的安全系數(shù)較其他答案略低，進(jìn)一步驗(yàn)證了該文SA-FOA算法用于邊坡臨界滑面搜索中是行之有效的。

表4 該文計(jì)算結(jié)果與評(píng)判員答案比較

4 結(jié)論

(1) 該文所提出的SA-FOA算法與瑞典條分法結(jié)合，在簡(jiǎn)單土坡中能夠快速、高效地搜索出臨界滑面，與單純采用SA算法或FOA算法相比，收斂速度更快，尋優(yōu)結(jié)果更準(zhǔn)確。

(2) SA-FOA算法與簡(jiǎn)化Bishop法結(jié)合，也適用于復(fù)雜邊坡的臨界滑面搜索，顯示出較強(qiáng)的全局搜索能力，得到的最小安全系數(shù)和臨界滑面更可靠。

(3) 將模擬退火算法與基本果蠅優(yōu)化算法相結(jié)合，克服了單一算法的局限性，經(jīng)驗(yàn)證具有較強(qiáng)的全局搜索能力和較快的收斂速度，為邊坡臨界滑面的搜索提供了一條高效的新途徑。