苗曉燕，陳昌富

(湖南大學(xué) 土木工程學(xué)院，湖南 長沙 410082)

1 引言

臨界滑動(dòng)面的合理確定是邊坡穩(wěn)定性分析的關(guān)鍵。傳統(tǒng)滑動(dòng)面形狀常人為假設(shè)為直線形、圓弧形、對(duì)數(shù)螺旋形、二次拋物線形等，而工程實(shí)踐中受邊坡幾何形狀各異、材料特性多變的影響，滑動(dòng)面形狀并不完全規(guī)則。因此，采用任意形狀滑動(dòng)面更具客觀性與廣泛適用性，且因滑動(dòng)面形狀所受約束較少，往往能搜索到更接近實(shí)際的臨界滑動(dòng)面和安全系數(shù)值。任意形狀滑動(dòng)面的構(gòu)造一般采用若干點(diǎn)將滑動(dòng)面曲線離散，點(diǎn)間通過直線或曲線相連，點(diǎn)的坐標(biāo)作為優(yōu)化變量。顯然，離散點(diǎn)越多滑動(dòng)面模擬精度越高，但相同模擬精度要求下，折線模式所需離散點(diǎn)個(gè)數(shù)要大于曲線模式，應(yīng)用曲線模擬任意形狀滑動(dòng)面相比折線具有優(yōu)化變量少、搜索效率高的優(yōu)勢(shì)。陳祖煜、李亮等，采用三次樣條插值函數(shù)較好地模擬了邊坡任意形狀滑動(dòng)面。但是，插值樣條要求曲線必須通過全部控制點(diǎn)，其合理性較大依賴于控制點(diǎn)本身的精確度。而控制點(diǎn)坐標(biāo)作為優(yōu)化變量，結(jié)果具有一定的隨機(jī)性，更多地代表的是一種趨勢(shì)和走向，此時(shí)，要求擬合曲線嚴(yán)格通過控制點(diǎn)意義不大，滑動(dòng)面的非凸性也難以保證，且計(jì)算公式較復(fù)雜。更合理的做法應(yīng)是，利用三次B樣條構(gòu)造一條曲線使之整體上最大程度地逼近控制點(diǎn)構(gòu)成的特征多邊形。B樣條是當(dāng)前工業(yè)幾何設(shè)計(jì)和形狀數(shù)學(xué)描述中表示和設(shè)計(jì)自由曲線的主流方法。三次B樣條具有三次插值樣條的光滑度，同時(shí)由于其曲線包含在控制多邊形內(nèi)，便于將控制多邊形作為邊界來約束滑動(dòng)面形狀。當(dāng)控制多邊形為凹形時(shí)，滑動(dòng)面滿足非凸性；當(dāng)相鄰4個(gè)控制點(diǎn)共線時(shí)，其對(duì)應(yīng)曲線段退化為一直線段，可模擬直曲組合滑動(dòng)面。其中，三次均勻B樣條因幾何性質(zhì)優(yōu)良、公式簡明易行應(yīng)用廣泛。

20世紀(jì)70年代以來，大量學(xué)者分別引入各種優(yōu)化方法搜索邊坡臨界滑動(dòng)面。它們相比傳統(tǒng)手動(dòng)試算或經(jīng)驗(yàn)方法，人工計(jì)算量大大減輕，成果更為準(zhǔn)確可靠。蟻群算法、粒子群算法、遺傳算法等啟發(fā)式算法以其強(qiáng)大的全局搜索能力、廣泛的適用性成為近年來的研究熱點(diǎn)。2009年，Rashedi等提出一種新的啟發(fā)式算法——萬有引力搜索算法(Gravitational search Algorithm，GSA)。該算法通過模擬物理學(xué)萬有引力定律實(shí)現(xiàn)種群優(yōu)化，它擁有強(qiáng)大的全局搜索能力，性能優(yōu)于常見的粒子群算法和遺傳算法。蔣建國等研究發(fā)現(xiàn)，GSA在邊坡臨界滑動(dòng)面搜索中具有很好的應(yīng)用前景，有待深入研究。GSA作為一種全局優(yōu)化算法，也存在易陷入局部最小的通病，結(jié)合其他算法思想和優(yōu)化策略進(jìn)行組合改進(jìn)是提升算法優(yōu)化性能的有效途徑。該文采用三次均勻B樣條函數(shù)模擬邊坡任意形狀滑動(dòng)面，采用Morgenstern-Price法(M-P法)計(jì)算滑動(dòng)面對(duì)應(yīng)安全系數(shù)。然后將啟發(fā)式算法中的GSA與PSO相融合，引入雁陣效應(yīng)和非線性慣性權(quán)重進(jìn)行改進(jìn)，提出一種新型混合啟發(fā)式搜索算法(Hybrid Heuristic Algorithm，HHA)，以安全系數(shù)為目標(biāo)函數(shù)，對(duì)邊坡臨界滑動(dòng)面進(jìn)行搜索。最后通過四道標(biāo)準(zhǔn)考題來驗(yàn)證該文方法的可靠性。

2 邊坡任意滑動(dòng)面模擬方法

2.1 控制點(diǎn)序列與優(yōu)化變量

將滑動(dòng)面用N個(gè)點(diǎn)P1、P2,…,PN離散(滑入點(diǎn)為P1，滑出點(diǎn)為PN)，作為控制點(diǎn)序列，則在平面內(nèi)，控制點(diǎn)的位置矢量共含2N個(gè)未知變量(xP1，yP1，…，xPN，yPN)。若直接將所有未知變量作為優(yōu)化變量，存在工作量大、效率低的問題。為精簡優(yōu)化變量，以坡腳為原點(diǎn)建立直角坐標(biāo)系，將控制點(diǎn)等橫間距布置，則當(dāng)x1、xN確定后，xi(i=2，3，…，N-1)均可知，且有y1=0，yN=H+(xN-H/tanβ)tanε，優(yōu)化變量降至N個(gè)，即(xP1，yP2，yP3，…，yPN-1，xPN)。

2.2 三次均勻B樣條曲線的生成

當(dāng)控制點(diǎn)坐標(biāo)確定后，代入三次均勻B樣條表達(dá)式即可生成樣條曲線模擬的滑動(dòng)面。三次均勻B樣條曲線由多個(gè)三次多項(xiàng)式構(gòu)成的曲線段分段拼接而成，每個(gè)曲線段由基函數(shù)與相鄰4個(gè)控制點(diǎn)的線性組合唯一確定。該表達(dá)式不含導(dǎo)數(shù)，較插值樣條函數(shù)更加簡單、便于實(shí)現(xiàn)。第i段的表達(dá)式為：

(1)

式中：Si(u)為第i段B樣條曲線上任一點(diǎn)的位置矢量，即[xi(u)，yi(u)]；Pi為控制點(diǎn)序列中第i點(diǎn)的位置矢量，即[xPi，yPi]；u為[0,1]區(qū)間內(nèi)參變量，以Δu的速度從0增至1，該文取Δu=0.001，即令每個(gè)曲線段由1 000個(gè)點(diǎn)組成。

2.3 曲線端點(diǎn)的處理

由于三次均勻B樣條曲線始于第二個(gè)控制點(diǎn)附近，終于倒數(shù)第二個(gè)控制點(diǎn)附近[圖1(a)]，故為保證曲線能完整擬合滑動(dòng)面，在離散點(diǎn)首尾分別補(bǔ)充一個(gè)邊界點(diǎn)P0、PN+1，控制點(diǎn)總數(shù)變?yōu)镹+2個(gè)。邊界點(diǎn)坐標(biāo)一般參照離散點(diǎn)的間距和方向設(shè)置，不影響優(yōu)化變量個(gè)數(shù)。但常規(guī)的邊界點(diǎn)設(shè)置只能使曲線端點(diǎn)位于首尾離散點(diǎn)附近，無法使二者重合，導(dǎo)致樣條曲線的端點(diǎn)不能剛好落在坡底和坡頂平面上。為解決這一問題，同時(shí)保持曲線在起點(diǎn)和終點(diǎn)的光順性，根據(jù)文獻(xiàn)[11]，該文給出一種較好的端點(diǎn)處理方法。

B樣條曲線的起點(diǎn)處，有S0(0)=P1，即：

P0+4P1+P2=6P1

(2)

設(shè)P2=η0P0，解得P0的位置矢量為：

(3)

B樣條曲線的終點(diǎn)處，有SN-2(1)=PN，即：

PN-1+4PN+PN+1=6PN

(4)

設(shè)pN+1=η1PN-1，解得PN+1的位置矢量為：

(5)

上述處理后，B樣條曲線的起、終點(diǎn)分別與首末離散點(diǎn)重合，且曲線保持光滑[圖1(b)]。

圖1 三次均勻B樣條曲線端點(diǎn)處理

2.4 滑動(dòng)面合理性控制

借鑒文獻(xiàn)[2]、[5]、[12]、[13]的滑面生成策略并結(jié)合B樣條曲線的凸包性，通過對(duì)控制點(diǎn)變量的初始化和邊界條件進(jìn)行約束，保證生成滑動(dòng)面的合理有效，無須二次修復(fù)或篩選。

如圖2所示，滑入點(diǎn)x1與滑出點(diǎn)xN按常規(guī)的常數(shù)上下限[x1min，x1max]、[xNmin，xNmax]約束。x1與xN在約束范圍內(nèi)隨機(jī)確定后，均分相鄰控制點(diǎn)的水平間距得Δx=(xN-x1)/(N-1)，則其余控制點(diǎn)的橫坐標(biāo)xi均可知。縱坐標(biāo)yi的約束條件較復(fù)雜，直觀來看，yi的上限不能超過坡面，下限不能低于基巖fs(xi)，若不存在基巖可人為設(shè)定一個(gè)較大負(fù)值，如-H等。設(shè)相鄰控制點(diǎn)間線段PiPi+1的水平傾角為αi(水平線下為負(fù))，為不生成鋸齒狀滑面，須保證αi≥αi-1，故tanαi≥tanαi-1，即(yi+1-yi)/Δx≥(yi-yi-1)/Δx，化簡得到y(tǒng)i+1≥2yi-yi-1。同時(shí)，為避免前面的搜索參數(shù)取某些值后，后面的搜索參數(shù)在約束條件下可行域?yàn)榭占那闆r，連接控制點(diǎn)Pi、PN交直線x=xi+1于Mi+1點(diǎn)，其縱坐標(biāo)為yMi+1，應(yīng)保證yi+1≤yMi+1，否則后面的滑動(dòng)面必然出現(xiàn)凹陷。根據(jù)相似三角形原理，(yMi+1-yi)/(yN-yi)=1/(N-i)，解得yi+1≤yi+(H-yi)/[N-(i-1)]。需要說明的是，由于已知y1為0，yi的搜索實(shí)際從i=2開始。此時(shí)，線段P1P2的水平傾角α1是滑動(dòng)面生成的初始方向，根據(jù)經(jīng)驗(yàn)，α1滿足-45°≤α1，所以-Δx≤y2，又因?yàn)榫€段P1PN交直線x=x2于x軸以上，故當(dāng)x2≤0時(shí)，只需不超過坡面y=0即可。綜上，yi的約束范圍如表1所示。

圖2 控制點(diǎn)取值范圍

表1 控制點(diǎn)縱坐標(biāo)yi的約束范圍

3 混合啟發(fā)式搜索算法(HHA)

3.1 混合啟發(fā)式搜索算法(HHA)思想

GSA算法中，粒子被視為一組在多維搜索空間無阻力運(yùn)動(dòng)的有質(zhì)量物體，粒子的位置代表優(yōu)化變量，每個(gè)位置對(duì)應(yīng)的優(yōu)化問題的目標(biāo)函數(shù)值稱為適應(yīng)值，粒子的質(zhì)量代表該位置適應(yīng)值的優(yōu)劣，適應(yīng)值越佳的粒子慣性質(zhì)量越大。粒子遵循萬有引力定律更新速度和位置，朝大質(zhì)量粒子趨近，逐漸逼近最優(yōu)位置對(duì)應(yīng)的最優(yōu)解。GSA具有較強(qiáng)的全局搜索能力和收斂速度，但也存在易陷入局部最小和解精度不高的問題。為此，該文結(jié)合粒子群算法PSO、雁陣效應(yīng)及非線性慣性權(quán)重對(duì)其進(jìn)行改進(jìn)，得到一種新型混合啟發(fā)式搜索算法(Hybrid Heuristic Algorithm，HHA)。

GSA每次迭代時(shí)，粒子的位置都會(huì)更新，粒子群到達(dá)過的較優(yōu)位置不能被保留。雖然迭代歷史上的種群最優(yōu)都會(huì)被記錄，對(duì)比后得出最終解，但在粒子移動(dòng)過程中，粒子的速度和位置更新僅依據(jù)當(dāng)前引力，不能學(xué)習(xí)粒子個(gè)體和群體的有益經(jīng)驗(yàn)，算法收斂速度較慢、解精度不高。而粒子群算法在速度的更新中，通過跟蹤個(gè)體最優(yōu)和群體最優(yōu)實(shí)現(xiàn)了對(duì)粒子記憶信息的有效利用。采用粒子群算法的記憶功能改進(jìn)萬有引力算法，粒子的速度公式為：

(6)

然而改進(jìn)后，由于受全局最優(yōu)的強(qiáng)烈吸引，粒子群在搜索后期趨向同一化。要想為群體輸入新鮮信息，就要增加種群數(shù)量或降低對(duì)全局最優(yōu)的追逐。增加種群數(shù)量會(huì)加大計(jì)算量，而減弱對(duì)全局最優(yōu)的追逐又使算法不易收斂。為在保持粒子多樣性的同時(shí)減少代價(jià)，受自然界雁陣效應(yīng)啟發(fā)，將粒子按個(gè)體最優(yōu)適應(yīng)值排序，最佳的粒子作為頭雁，依靠其自身飛行經(jīng)驗(yàn)決定下一步的運(yùn)動(dòng)，其他粒子依次跟隨其前方較優(yōu)粒子飛行。每次迭代后，更新粒子群的歷史最優(yōu)適應(yīng)值，并重新排序。這樣就使粒子群不止朝一個(gè)方向運(yùn)動(dòng)，保持了粒子的多樣性。綜合雁陣效應(yīng)和粒子群記憶功能，最終可得混合啟發(fā)式搜索算法(HHA)的速度公式為：

(7)

由式(7)可知，粒子當(dāng)下速度以r0的比例繼承于先前速度，這種繼承能力的大小影響算法的全局探索和局部開發(fā)能力。當(dāng)r0較大時(shí)，算法全局探索能力增強(qiáng)，但局部開發(fā)能力減弱，相同迭代次數(shù)的最優(yōu)解精度較差；當(dāng)r0較小時(shí)，算法局部開發(fā)能力增強(qiáng)，但全局探索能力減弱，在優(yōu)化前期易陷入局部最優(yōu)。為更好地平衡算法的全局探索能力和局部開發(fā)能力，引入非線性慣性權(quán)重函數(shù)代替隨機(jī)數(shù)來體現(xiàn)繼承能力：

r0=ws-(ws-we)(t/T)2

(8)

式中：ws為初始慣性權(quán)重；we為最大迭代次數(shù)的慣性權(quán)重，通常取ws=0.9、we=0.4時(shí)算法優(yōu)化性能最好；t為當(dāng)前迭代次數(shù)；T為最大迭代次數(shù)。由式(8)可知，r0隨迭代次數(shù)逐漸下降，且前期下降較慢，后期下降較快。由此可保證算法在迭代前期著重于全局搜索，而后期著重于局部開發(fā)。

3.2 混合啟發(fā)式搜索算法(HHA)流程

(1)確定搜索空間，生成初始種群。設(shè)在d維搜索空間中存在N個(gè)粒子，第i個(gè)粒子的位置為：

(9)

(3)計(jì)算并更新粒子質(zhì)量。定義第i個(gè)粒子t時(shí)刻的慣性質(zhì)量為：

(10)

式中：fiti(t)為第i個(gè)粒子t時(shí)刻的適應(yīng)值；best(t)、worst(t)分別為t時(shí)刻種群最優(yōu)和最差適應(yīng)值。

(4)計(jì)算并更新萬有引力。由萬有引力定律，t時(shí)刻粒子i受粒子j在第d維上的引力為：

(11)

Ri,j(t)=‖xi(t)，xj(t)‖2

(12)

G(t)為t時(shí)刻萬有引力常數(shù)，其值隨宇宙實(shí)際年齡的增大而變小，表達(dá)式為：

G(t)=G0×e-αt/T

(13)

式中：G0為t0時(shí)刻的G值，通常取100；α為衰減系數(shù)，通常取20；T為最大迭代次數(shù)。

(5)計(jì)算合力。t時(shí)刻，慣性質(zhì)量最大的前k個(gè)粒子集合作用于粒子i第d維上的合力為：

(14)

式中：randj為[0,1]之間的隨機(jī)數(shù)，可使搜索具有一定的隨機(jī)性而更加合理；kbest為慣性質(zhì)量最大的前k個(gè)粒子集合，其初值為N，隨時(shí)間t線性減至1。

(6)更新粒子的加速度、速度和位置。由牛頓第二定律，t時(shí)刻粒子i在第d維產(chǎn)生的加速度為：

(15)

速度更新公式見式(7)、(8)，位置更新公式為：

(16)

3.3 數(shù)值試驗(yàn)驗(yàn)證

為評(píng)估HHA的優(yōu)化性能與效率，對(duì)比文獻(xiàn)[10]中8個(gè)經(jīng)典測(cè)試函數(shù)進(jìn)行數(shù)值試驗(yàn)，結(jié)果如表2所示。MGSA為文獻(xiàn)[10]改進(jìn)GSA算法。表中，F(xiàn)1、F2、F3、F4為單峰函數(shù)，其余為多峰函數(shù)。函數(shù)維數(shù)均為30，全局最小真值均為0。參數(shù)G0=100，N=50，α=40，itmax=1 000。每個(gè)函數(shù)計(jì)算30次。

表2 測(cè)試函數(shù)不同算法優(yōu)化性能對(duì)比

注：表中GSA、MGSA的計(jì)算結(jié)果源自文獻(xiàn)[10]。

對(duì)比表2均值可知：HHA的優(yōu)化精度較GSA有數(shù)量級(jí)的提升，改進(jìn)效果優(yōu)于MGSA。說明HHA能有效跳出局部最優(yōu)，具有更好的收斂精度。再對(duì)比表2標(biāo)準(zhǔn)差可知，HHA的標(biāo)準(zhǔn)差最小，穩(wěn)定性優(yōu)于GSA與MGSA。

由于HHA相比GSA增加了計(jì)算步驟，為評(píng)估HHA算法的效率，以F1函數(shù)為例，比較等精度要求下的算法耗時(shí)。設(shè)誤差要求控制在5%以內(nèi)，HHA只需迭代50次，耗時(shí)0.96 s，而GSA需迭代500次以上，耗時(shí)大于3.6 s，可見等精度要求下HHA優(yōu)化效率更高。

綜上所述，HHA能有效克服GSA易陷入局部最小的不足，具有更好的優(yōu)化精度、效率與穩(wěn)定性。

3.4 基于HHA的邊坡穩(wěn)定性分析

由三次均勻B樣條函數(shù)生成有效的任意形狀滑動(dòng)面后，采用計(jì)算嚴(yán)密、收斂性好的M-P法計(jì)算滑面對(duì)應(yīng)的安全系數(shù)Fs。再以安全系數(shù)為目標(biāo)函數(shù)，根據(jù)3.2節(jié)流程采用HHA算法搜索得到邊坡臨界滑動(dòng)面及對(duì)應(yīng)的最小安全系數(shù)。

M-P法中邊坡的條分方式和條塊受力如圖3所示。目標(biāo)函數(shù)Fs的計(jì)算公式為：

圖3 邊坡土條受力分析

(17)

式中：

(18)

其中比例因子λ按下式計(jì)算：

λ=

(19)

而Ei通過以下遞推公式得到：

(20)

4 算例驗(yàn)證

選取澳大利亞計(jì)算機(jī)應(yīng)用協(xié)會(huì)(ACADS)的4道邊坡穩(wěn)定標(biāo)準(zhǔn)考題，應(yīng)用樣條逼近生成滑動(dòng)面并采用M-P法計(jì)算安全系數(shù)，再以安全系數(shù)為目標(biāo)函數(shù)，分別采用GSA與HHA搜索臨界滑動(dòng)面，與推薦答案進(jìn)行對(duì)比。其中，考題1(a)為均質(zhì)邊坡，考題1(c)為非均質(zhì)邊坡，考題1(d)在考題1(c)基礎(chǔ)上土體水平方向增加一系數(shù)為0.15g的地震加速度，為考慮地震效應(yīng)邊坡，考題3(a)為含軟弱夾層邊坡。

題中邊坡的幾何形狀與剖面尺寸見圖4，材料特性如表3所示。為平衡搜索效率與滑動(dòng)面模擬精度，考題1(a)中離散點(diǎn)取7個(gè)，考題1(c)、1(d)中離散點(diǎn)取9個(gè)，考題3(a)中離散點(diǎn)取11個(gè)，條分?jǐn)?shù)均為20。GSA與HHA中種群數(shù)為50，總迭代次數(shù)為200，HHA中c1=c2=c3=0.5。計(jì)算結(jié)果及推薦答案列于表4。

由表4可知：考題1(a)、1(c)、1(d)中，該文方法的計(jì)算安全系數(shù)均與文獻(xiàn)推薦答案較為接近，相應(yīng)的臨界滑動(dòng)面與文獻(xiàn)圓弧滑動(dòng)面形狀相似、位置相近，證明三次均勻B樣條能合理逼近圓弧滑動(dòng)面，GSA與HHA對(duì)于均質(zhì)邊坡、非均質(zhì)邊坡及考慮地震效應(yīng)邊坡的滑動(dòng)面搜索均適用。且HHA的優(yōu)化結(jié)果比GSA更小，優(yōu)化性能更佳。

圖4 臨界滑動(dòng)面對(duì)比

表3 考題材料參數(shù)

考題3(a)為典型的滑動(dòng)面沿軟弱夾層呈非圓弧狀的邊坡，由于該文不對(duì)滑動(dòng)面直線段預(yù)先假定，且土體的軟弱夾層很薄(0.5 m)，控制點(diǎn)個(gè)數(shù)又多(11個(gè))，故臨界滑裂面的搜索非常困難。此時(shí)，GSA顯然難以勝任，最優(yōu)解停留在1.6左右，遠(yuǎn)大于推薦答案，滑動(dòng)面位置也相距較遠(yuǎn)。而HHA的安全系數(shù)依然接近并略小于推薦答案，二者滑動(dòng)面位置相近，形狀吻合，特別是在軟弱土層交界處搜索得到較長距離的直線段，證明樣條曲線能合理逼近非圓弧狀滑動(dòng)面。HHA對(duì)含軟弱夾層邊坡和多控制點(diǎn)的復(fù)雜滑動(dòng)面搜索依然有效。

表4 考題答案對(duì)比

5 結(jié)論

(1)基于三次均勻B樣條曲線構(gòu)造邊坡任意形狀滑動(dòng)面，探討了曲線端點(diǎn)的處理方法與有效滑動(dòng)面生成策略，有效地減少了優(yōu)化變量，提高了搜索效率，保證了生成滑動(dòng)面的形狀自由與合理有效。

(2)將萬有引力算法與粒子群算法相融合，并引入雁陣效應(yīng)、非線性慣性權(quán)重進(jìn)行改進(jìn)，提出一種新的混合啟發(fā)式搜索算法(HHA)，通過與其他算法進(jìn)行對(duì)比，驗(yàn)證了HHA的優(yōu)良性能。

(3)通過4道邊坡穩(wěn)定標(biāo)準(zhǔn)考題，分析驗(yàn)證了該文提出的樣條滑動(dòng)面和HHA搜索算法在邊坡穩(wěn)定分析應(yīng)用中的可靠性與有效性，發(fā)現(xiàn)三次均勻B樣條曲線能合理逼近任意形狀臨界滑動(dòng)面；HHA對(duì)臨界滑動(dòng)面的搜索優(yōu)于GSA，對(duì)非均質(zhì)、含軟弱夾層及地震作用下的復(fù)雜邊坡均具有適用性。