基于改進(jìn)螢火蟲優(yōu)化算法和考慮空間變異特性的邊坡最小安全系數(shù)確定研究

褚雪松，李 亮

(青島理工大學(xué) 土木工程學(xué)院，山東 青島 266033)

0 引 言

確定邊坡的最小安全系數(shù)(失穩(wěn)模式或者滑動面)是邊坡穩(wěn)定分析以及邊坡可靠度分析中非常重要的步驟。對于邊坡確定性分析而言，需要知道邊坡最可能的滑動面以便為邊坡設(shè)計與加固提供決策依據(jù)；對于邊坡可靠度分析而言，需要知道蒙特卡洛每個樣本下對應(yīng)的確定性分析中可能的最小安全系數(shù)(失穩(wěn)模式)，以便為邊坡失穩(wěn)的可靠度分析(風(fēng)險分析)提供定量支持。實際上，歸根結(jié)底，在輸入?yún)?shù)確定的前提下，如何確定邊坡最小安全系數(shù)是非常古老而又關(guān)鍵的問題。在邊坡穩(wěn)定確定性分析中，許多學(xué)者均開展過卓有成效的研究工作，陳祖煜等[1]采用單純形等方法確定邊坡的最小安全系數(shù)，鄒廣電[2]采用非線性規(guī)劃原理結(jié)合Bishop法對邊坡的最小安全系數(shù)進(jìn)行了搜索，香港理工大學(xué)鄭榕明博士自主研發(fā)了SLOPE2000邊坡穩(wěn)定分析軟件，集成了包括模擬退火算法[3]、和聲搜索算法[4-5]、粒子群算法[6-7]在內(nèi)的多種優(yōu)化算法，其他方法譬如復(fù)合形法[8]、禁忌魚群算法[9]、遺傳算法[10]、蟻群算法[11-12]、網(wǎng)魚算法[13-14]等都應(yīng)用于邊坡最小安全系數(shù)的搜索中，取得了一系列的研究成果。上述研究多采用均質(zhì)邊坡或者多層含軟弱層的邊坡為例，然而在考慮土體材料空間變異特性時，邊坡的材料參數(shù)一般采用隨機(jī)場模擬，對于每一個隨機(jī)場的實現(xiàn)，將其參數(shù)輸入到邊坡確定性分析模型，然后進(jìn)行最小安全系數(shù)的搜索，通過判斷該最小安全系數(shù)是否小于1來決定該樣本是否屬于失效樣本，該最小安全系數(shù)對應(yīng)的滑動面積或者體積作為邊坡失穩(wěn)風(fēng)險分析的重要參數(shù)，因此非常有必要研究土體材料參數(shù)隨機(jī)場模擬前提下如何才能有效地確定邊坡的最小安全系數(shù)及失穩(wěn)模式。

1 邊坡最小安全系數(shù)的確定

邊坡穩(wěn)定分析的確定性模型可以表述為如下的數(shù)學(xué)模型

minFS{Gx,Px,Fx,Sx}

(1)

式中：FS為安全系數(shù)；Gx為邊坡的幾何模型；Px為邊坡土層的材料參數(shù)，包括密度ρ、內(nèi)摩擦角φ以及黏聚力c；Fx為邊坡所承受的荷載，譬如自重，堆載等；Sx為邊坡失穩(wěn)模式相關(guān)的數(shù)學(xué)變量，譬如圓弧滑動面的圓心或者一系列點的坐標(biāo)等。

圖1為邊坡最小安全系數(shù)確定模型，其中，ρi，ci，φi分別為第i個土層土體材料的密度、黏聚力以及內(nèi)摩擦角。邊坡最小安全系數(shù)的確定就是在給定的邊坡幾何形狀Gx、給定的土層材料參數(shù)Px、給定的荷載Fx下，通過不斷地變換可能的失穩(wěn)模式Sx，以期尋找到最小的安全系數(shù)及其對應(yīng)的失穩(wěn)模式。上述介紹的所有優(yōu)化方法都是不同的變換Sx策略，失穩(wěn)模式可以假定為圓弧或者非圓弧，作為一般情況，本文假定失穩(wěn)模式為非圓弧，由一系列的點通過樣條函數(shù)插值得到，限于篇幅，本文不再詳細(xì)描述非圓弧滑動面的模擬步驟，讀者可參考文獻(xiàn)[15]，[16]。安全系數(shù)的計算方法采用Spencer法[17]。

2 空間變異特性描述

土體材料的空間變異特性是土體材料不確定性來源中的一種，該特性的形成與土體材料在漫長的地質(zhì)年代中所經(jīng)受的沉積作用及地質(zhì)環(huán)境有關(guān)，通常認(rèn)為即使是在宏觀均質(zhì)的邊坡土層內(nèi)，任意不同兩點之間的材料強(qiáng)度參數(shù)或多或少都存在一定的變異性。在常規(guī)的邊坡穩(wěn)定分析中，均質(zhì)土層假定實際上是假定了任意不同兩點之間的材料強(qiáng)度參數(shù)不存在變異性或者說是完全相關(guān)的。研究表明[18]，土體材料的空間變異特性是土體材料的基本屬性，在邊坡可靠度分析中忽略該屬性會得到失真的邊坡失效概率，在某些情況下，失真程度會很大。由于Phoon等[18]研究表明，邊坡土層的水平變異特性較弱，因此本文僅考慮了土層的垂直變異特性。描述空間變異特性最常用的方法首推Vanmarcke[19]提出的隨機(jī)場模擬方法，下面以圖2所示的邊坡為例，簡單介紹一維隨機(jī)場模擬的步驟。圖2所示為一宏觀均質(zhì)的邊坡，該均質(zhì)邊坡被離散為12個水平的小土層C1～C12，在每一個小土層內(nèi)，材料參數(shù)是沒有空間變異特性的，或者說是完全相關(guān)的，不同的小土層之間，材料參數(shù)存在變異特性或者說具有一定的相關(guān)性，隨機(jī)場認(rèn)為，這種相關(guān)性隨著小土層相互之間距離的增加而衰減，譬如C1與C2之間的相關(guān)性要比C1與C12之間的相關(guān)性要強(qiáng)。具體的衰減公式為[19]

(1)

式中：ρCiCj為Ci與Cj之間的相關(guān)系數(shù)；yi為土層Ci中心的y坐標(biāo)；yj為土層Cj中心的y坐標(biāo)；λ為描述土層空間變異特性的參數(shù)，稱之為波動范圍，λ值越大，則相應(yīng)土層的空間變異特性越弱，λ→+∞時，則變?yōu)榫|(zhì)土層。

在得到相關(guān)系數(shù)矩陣之后，可以利用Cholesky分解結(jié)合概率密度函數(shù)得到隨機(jī)場的實現(xiàn)，隨機(jī)場的每一個實現(xiàn)就是確定了一組C1～C12小土層的強(qiáng)度參數(shù)值，譬如表1中每一行即為一個隨機(jī)場的實現(xiàn)，記為Ri。

3 改進(jìn)螢火蟲算法

螢火蟲算法[20-23]是一種仿生算法，算法主要是模仿螢火蟲群體在獵食等過程中根據(jù)生物螢光體發(fā)光亮度的強(qiáng)弱互相吸引來調(diào)整螢火蟲個體之間的飛行路線，從而達(dá)到群體之間的信息交流，實現(xiàn)種群的和諧與發(fā)展。將螢火蟲個體比擬為優(yōu)化問題中的一個解向量，個體之間位置的調(diào)整比擬為算法的尋優(yōu)過程，個體發(fā)光的亮度比擬為優(yōu)化問題的目標(biāo)函數(shù)(譬如本文中的安全系數(shù)FS)，該仿生算法遵循以下規(guī)則：

(1)螢火蟲個體沒有性別之分。

(2)螢火蟲的亮度與分析問題的目標(biāo)函數(shù)有關(guān)，對于最小安全系數(shù)的搜索問題而言，安全系數(shù)越小，則該個體的亮度就越大。

(3)螢火蟲個體之間的吸引度隨著個體之間距離的增加而遞減。

(1)算法初始化，確定變量N,β,以及最大迭代次數(shù)It；隨機(jī)生成N個螢火蟲個體，計算各自的亮度Bi。

(2)對N個螢火蟲群體中的每一個個體Xi進(jìn)行如下操作。判斷是否存在比該個體亮度更高的個體，若有，則記為Xj，利用公式(2)調(diào)整Xi的位置并更新其亮度；如果群體中不存在比該個體更亮的個體，則采用公式(3)隨機(jī)調(diào)整Xi的位置并更新其亮度。

(2)

(3)

式中：r1～r5均為[0,1]之間的隨機(jī)數(shù)；uk，lk分別為描述數(shù)學(xué)變量Sx所需第k個變量的取值上、下限，k=1，2，…，n。

(3)重復(fù)執(zhí)行步驟(2)共計It次，然后算法停止，從螢火蟲群體中找出最亮的個體作為優(yōu)化問題的解。

需要說明的是，基本的螢火蟲算法在采用公式(2)，(3)更新完螢火蟲個體的位置之后，沒有判斷新位置的亮度是否優(yōu)于當(dāng)前位置的亮度，就從當(dāng)前位置移動到新位置。實際上，在優(yōu)化過程中，為了使螢火蟲個體的亮度快速增加，可以采取一種“短視”策略，即只有新位置的亮度比原來位置亮度增加時，螢火蟲個體才會移動至該新位置，否則保持原位置不動。在基本的螢火蟲算法中是沒有采取“短視”策略的。這種策略通常會引起算法陷入局部最優(yōu)，但是在算法的初期會提高算法的搜索效率。

表1 2種波動范圍下的隨機(jī)場實現(xiàn)Tab.1 Realization of Random Field at Two Fluctuation Ranges

4 算法影響參數(shù)分析

本算例考慮圖2所示的邊坡，坡高10 m，考慮不排水工況下邊坡的短期穩(wěn)定性，強(qiáng)度參數(shù)取為不排水強(qiáng)度cu，即考慮了黏聚力與內(nèi)摩擦角的綜合強(qiáng)度。假設(shè)cu符合對數(shù)正態(tài)分布，其均值為40 kPa，變異系數(shù)為0.3，即標(biāo)準(zhǔn)差為12 kPa，cu的空間變異特性用一維(垂直)隨機(jī)場模擬，采用厚度為1 m的水平土層離散整個邊坡，同時考慮離散尺寸為1 m時的方差折減系數(shù)以及對應(yīng)的相關(guān)系數(shù)，具體公式請參考文獻(xiàn)[19]。

4.1 “短視”策略及β對算法的影響

首先假定波動范圍λ=1 m，隨機(jī)抽取10個隨機(jī)場的實現(xiàn)樣本Ri，其對應(yīng)的12個小土層的cu值如表1所示，將每個隨機(jī)場實現(xiàn)確定的12個cu值視為邊坡確定性分析中的輸入?yún)?shù)，利用螢火蟲算法搜索其相應(yīng)的最小安全系數(shù)。由于螢火蟲算法為隨機(jī)搜索算法，因此隨機(jī)生成10組不同的螢火蟲群體，利用10組不同的螢火蟲群體計算得到的結(jié)果統(tǒng)計算法的均值和標(biāo)準(zhǔn)差，利用均值與標(biāo)準(zhǔn)差的來衡量算法的優(yōu)越性，一般來看，均值越小，算法的尋優(yōu)能力越強(qiáng)；在平均值相等的前提下，標(biāo)準(zhǔn)差越小，則算法的尋優(yōu)能力越強(qiáng)。

螢火蟲算法參數(shù)設(shè)置：個體數(shù)目N=30，It=50，β=0.0,0.2,0.5,1.0,2.0,5.0,10.0,20.0,50.0,100.0,200.0,500.0,1 000.0，分別結(jié)合是否采取“短視”策略，探討了本文提出的“短視”策略對計算結(jié)果的影響，并給出了13個不同的β值對每一個隨機(jī)場實現(xiàn)下邊坡最小安全系數(shù)的影響(圖3)。由圖3可見，在不同的β取值下，算法得到的安全系數(shù)平均值均會有不同程度的起伏，在β>50以后，算法結(jié)果趨于穩(wěn)定。從安全系數(shù)平均值來看，對于第1，2，6，9個隨機(jī)場實現(xiàn)而言，β=0.2時算法得到了最小的安全系數(shù)平均值；對于第3，4，5個隨機(jī)場實現(xiàn)而言，β=0.5時算法得到了最小的安全系數(shù)平均值；對于第7，10個隨機(jī)場實現(xiàn)而言，β=0.0時算法得到了最小的安全系數(shù)平均值；對于第8個隨機(jī)場實現(xiàn)而言，β=2.0時算法得到了最小的安全系數(shù)平均值。由此可見，對于考慮空間變異特性下不同的隨機(jī)場實現(xiàn)而言，必須要變化算法的參數(shù)設(shè)置以便能夠找到最小的安全系數(shù)，一般不存在對所有隨機(jī)場都能找到最小安全系數(shù)的參數(shù)組合。β的合理取值范圍建議為[0,2]。圖4為“短視”策略對算法的影響。從圖4(a)～(c)可知，基本的螢火蟲算法沒有采取“短視”策略，安全系數(shù)的平均值較大，說明采取“短視”策略能夠增強(qiáng)螢火蟲算法的尋優(yōu)能力。從圖4(d)～(e)可知，采取“短視”策略后，螢火蟲算法結(jié)果的標(biāo)準(zhǔn)差較未采取“短視”策略的結(jié)果標(biāo)準(zhǔn)差要小(β=0.2時，第6個隨機(jī)場實現(xiàn)及β=0.5時第4個隨機(jī)場實現(xiàn)除外)，再次證明了采取“短視”策略能夠提高算法的尋優(yōu)能力。本文還比較了采取“短視”策略下不同β取值對計算結(jié)果的影響。

4.2 N對算法的影響

設(shè)置參數(shù)如下：It=50，N=10,30,50，β=0.5，采取“短視”策略，計算結(jié)果的比較如圖5所示。由圖5(a)可知，隨著群體中螢火蟲個體數(shù)目的增加，安全系數(shù)的平均值減小，這說明隨著N的增加，螢火蟲算法的尋優(yōu)能力提高。從圖5(b)中的標(biāo)準(zhǔn)差比較可知，有8個隨機(jī)場實現(xiàn)的最小標(biāo)準(zhǔn)差出現(xiàn)在N=50，另各有1個出現(xiàn)在N=10和N=30，綜合比較來看，隨著N的增加，算法的尋優(yōu)能力增強(qiáng)。圖5(c)還給出了算法在不同N時的計算時間。圖5(c)中計算時間指的是算法對10種隨機(jī)場實現(xiàn)全部的耗時，由圖5(c)可見，計算時間T與N之間符合乘冪函數(shù)關(guān)系，T與N2.1呈正比關(guān)系，所以N的取值不宜過大。必須要說明的是，算法的計算時間與計算機(jī)的配置有關(guān)，本文算法均在Intel(R) CoreTMi7-2600 CPU@3.40GHz,內(nèi)存16.0G的臺式機(jī)上運行。

4.3 λ對算法的影響

假定波動范圍λ=1 000 m，重新計算了螢火蟲算法(N=30,It=50,β=0.0,0.2,0.5)在表1中所示10個隨機(jī)場實現(xiàn)中的表現(xiàn)，統(tǒng)計了安全系數(shù)的均值和標(biāo)準(zhǔn)差。由表1可知，λ=1 000 m時10個隨機(jī)場的實現(xiàn)與λ=1 m時的10個隨機(jī)場實現(xiàn)不同，所以其安全系數(shù)平均值不具有可比性?？梢栽O(shè)想，算法參數(shù)設(shè)置相同時，算法應(yīng)該會具有相同的尋優(yōu)能力，如果算法所得安全系數(shù)的標(biāo)準(zhǔn)差較小，則間接反映了搜索空間的復(fù)雜度較小。圖6為λ=1 m和1 000 m時安全系數(shù)標(biāo)準(zhǔn)差比較。由圖6可見，當(dāng)移動步長β=0.0，λ=1 000 m時，螢火蟲算法得到的標(biāo)準(zhǔn)差有8次小于λ=1 m時的標(biāo)準(zhǔn)差，β=0.2時，有6次小于λ=1 m時的標(biāo)準(zhǔn)差，β=0.5時，有8次小于λ=1 m時的標(biāo)準(zhǔn)差。由此可見，在波動范圍增加后，邊坡最小安全系數(shù)這一優(yōu)化問題對應(yīng)的搜索空間變得簡單。

5 結(jié) 語

在基本螢火蟲算法基礎(chǔ)上，引入了“短視”策略，通過不排水邊坡的短期穩(wěn)定問題，探討了不同算法參數(shù)下螢火蟲算法在考慮土體材料空間變異特性時的尋優(yōu)能力。計算結(jié)果表明，引入“短視”策略可以大大提高算法的尋優(yōu)能力，增加螢火蟲個體數(shù)目N也可以提高算法的尋優(yōu)能力，但是其計算時間與N2.1成正比，為節(jié)省計算時間，N取值不宜過大。移動步長β對結(jié)果的影響較大，對于不同的隨機(jī)場實現(xiàn)而言，不同的β取值可以得到最小的安全系數(shù)平均值，因此，在考慮空間變異特性時的邊坡穩(wěn)定分析中，必須要變換算法的參數(shù)以便能夠保證搜索到最小的安全系數(shù)。當(dāng)土體材料空間變異特性變?nèi)?或波動范圍增加)時，邊坡最小安全系數(shù)這一優(yōu)化問題的搜索空間復(fù)雜度變?nèi)酢?/p>