李生清

(山東省地質(zhì)礦產(chǎn)勘查開(kāi)發(fā)局八○一水文地質(zhì)工程地質(zhì)大隊(duì)(山東省地礦工程勘察院)， 山東濟(jì)南 250014)

0 引言

自然環(huán)境變遷以及人類(lèi)社會(huì)過(guò)度開(kāi)發(fā)，給土地、山體以及坡地帶來(lái)較大的負(fù)面影響，不穩(wěn)定的土壤條件、松散的地層以及活躍的大陸板塊，使坡地出現(xiàn)了嚴(yán)重程度不同的滑裂面，影響自然生態(tài)安全。在過(guò)去的若干年間，國(guó)內(nèi)外大量學(xué)者對(duì)于邊坡穩(wěn)定問(wèn)題的研究一直很關(guān)注(祝志恒，2009；佘小年等，2011；樊成等，2020；鄧宇軒等，2021；宋二祥，2021)。目前應(yīng)用較多的方法為剛體極限平衡法，即通過(guò)簡(jiǎn)單計(jì)算將邊坡穩(wěn)定問(wèn)題轉(zhuǎn)化為二維問(wèn)題。盡管一些方法能夠快速獲得滑裂面，但對(duì)于滑裂面的危險(xiǎn)程度卻難以區(qū)分。從20世紀(jì)中后期開(kāi)始，越來(lái)越多的學(xué)者將邊坡穩(wěn)定性問(wèn)題轉(zhuǎn)化為三維問(wèn)題，通過(guò)三維分析模型以及智能化、自動(dòng)化系統(tǒng)搜索最不利滑裂面，但三維滑裂面本身較為復(fù)雜，因此該項(xiàng)工作仍存在較大困難(溫樹(shù)杰等，2018；韓林和蔡永昌，2020；石峰等，2021)。

針對(duì)幾種常規(guī)搜索方法的不足，吳志軒等(2019)中提出利用Bishop方法和有限元方法搜索邊坡最不利滑裂面，從無(wú)限長(zhǎng)邊坡基-填界面開(kāi)挖臺(tái)階情況出發(fā)，利用Bishop 方法推導(dǎo)不同破壞模式下的安全系數(shù)公示，提出了理論解法并解釋了穩(wěn)定性強(qiáng)化機(jī)理。方宏偉等(2019)基于SLFT(滑移線場(chǎng)理論)設(shè)計(jì)邊坡雙折減系數(shù)強(qiáng)度失穩(wěn)判據(jù)，根據(jù)極限狀態(tài)下的邊坡數(shù)據(jù)分析結(jié)果，討論邊坡局部位置是否開(kāi)始出現(xiàn)裂紋。該方法利用整體安全系數(shù)定位發(fā)生問(wèn)題的位置，實(shí)現(xiàn)對(duì)邊坡穩(wěn)定性、安全性的基本判斷，并且相較于已有方法，提出了失穩(wěn)判據(jù)的客觀標(biāo)準(zhǔn)，排除了人為主觀因素的影響。余鵬等(2022)基于XFEM構(gòu)造模擬土質(zhì)邊坡失穩(wěn)滑裂過(guò)程的方法，將邊坡失穩(wěn)過(guò)程概化為滑裂面的發(fā)生、擴(kuò)展和貫通的過(guò)程?？椎轮镜?2022)基于極限平衡理論，建立目標(biāo)函數(shù)與變量的關(guān)系表達(dá)式,通過(guò)求駐點(diǎn)的方法求解各條塊的坐標(biāo)值確定其相應(yīng)的穩(wěn)定性系數(shù)和滑裂面。盧坤林等(2009)基于滑面正應(yīng)力修正模式的假定，根據(jù)邊坡的平衡條件，導(dǎo)出含安全系數(shù)的平衡方程組,通過(guò)反復(fù)迭代求解安全系數(shù)。研究結(jié)果表明：該方法能夠得到更小的安全系數(shù)，最危險(xiǎn)滑面所對(duì)應(yīng)的冪次隨著滑體長(zhǎng)度的增大先略微降低后逐漸增大。然而，上述傳統(tǒng)方法在選取潛在滑裂面關(guān)鍵點(diǎn)時(shí)，未綜合考量剪切滑裂點(diǎn)、受拉滑裂點(diǎn)對(duì)搜索過(guò)程的影響程度，導(dǎo)致搜索過(guò)程易陷入局部最優(yōu)，從而丟失其他同等級(jí)搜索數(shù)據(jù)。

針對(duì)傳統(tǒng)方法存在的不足，本研究針對(duì)GA-Sarma算法的遺傳性和概率性特點(diǎn)，進(jìn)一步優(yōu)化邊坡滑裂面搜索方法，通過(guò)GA-Sarma算法搜索邊坡最不利滑裂面，為邊坡支護(hù)、加固及改造等地質(zhì)工作提供更加可靠的治理技術(shù)。

1 基于GA-Sarma算法的邊坡最不利滑裂面搜索方法

1.1 計(jì)算穩(wěn)定性系數(shù)提取潛在滑裂面

根據(jù)一般滑裂面的基本結(jié)構(gòu)建立三維勻質(zhì)邊坡模型，以其中坡腳線中點(diǎn)為原點(diǎn)，設(shè)邊坡水平面為x軸，與邊坡水平面垂直的坡腳線設(shè)為y軸，過(guò)原點(diǎn)垂直于兩軸的坐標(biāo)軸為z軸。設(shè)置滑面方程，即存在：

z=f(x,y)

(1)

設(shè)置邊坡土體粘聚力為a，內(nèi)摩擦角為θ，土體重度為b，則可以做出三組假設(shè)：

(1)假設(shè)邊坡發(fā)生變形時(shí)滑體為剛體(楊山奇等，2018；劉鋒濤等，2019)；

(2)假設(shè)求解滑裂面存儲(chǔ)勢(shì)能時(shí)，變形后的微面變形大小可利用剛度系數(shù)dij模擬，則剛度系數(shù)與微面面積之間，存在正比例關(guān)系；

(3)假設(shè)外力作用產(chǎn)生虛位移能夠減少滑裂面存儲(chǔ)勢(shì)能，則外力與虛位移分別為：

(2)

(3)

式中：λ表示隨機(jī)微面單元外法線的單位向量，該值可通過(guò)公式(4)獲得。由于每一個(gè)微面單元不一定是完全平整的曲面，為了便于計(jì)算與研究，默認(rèn)微面單元為光滑平面，則根據(jù)微面坐標(biāo)軸顯示的坐標(biāo)可得到該參數(shù)的取值：

(4)

(5)

確定土體單元受到的法向力計(jì)算剪應(yīng)力S1(李珂等，2018)。由于該值的計(jì)算過(guò)程與上述彈性勢(shì)能高度近似，因此不再重復(fù)說(shuō)明，默認(rèn)剪切勢(shì)能為S2，則求解獲得的總勢(shì)能為：

(6)

(7)

式中：u、v、w表示函數(shù)求解參數(shù)。設(shè)置抗滑力參數(shù)為G1，設(shè)置邊坡滑力參數(shù)為G2，則穩(wěn)定性系數(shù)的一般計(jì)算表達(dá)式為：

G=G1/G2

(8)

G1和G2的值兩組參數(shù)的計(jì)算表達(dá)式為：

(9)

式中：li表示第i處的活動(dòng)距離；c表示滑動(dòng)系數(shù)；α表示剪切力與虛位移之間夾角(孫加平等，2017)。公式(9)是在滑動(dòng)系數(shù)的控制下，根據(jù)不同方向上的作用效果和邊坡活動(dòng)距離，對(duì)抗滑力和滑力參數(shù)展開(kāi)計(jì)算。將公式(9)代入公式(8)中，求出邊坡穩(wěn)定性系數(shù)，然后可根據(jù)該系數(shù)提取潛在滑裂面。

1.2 選取潛在滑裂面關(guān)鍵點(diǎn)

隨著大位移發(fā)生形成潛在滑裂面，無(wú)論在靜力還是動(dòng)力條件下，失穩(wěn)時(shí)的坡體存在大面積破壞的情況(黃詩(shī)淵等，2018；馮晴楓等，2020)。已知巖土內(nèi)部任意平面的剪應(yīng)力超過(guò)極限剪應(yīng)力，致使剪切應(yīng)變?cè)隽堪l(fā)生突變現(xiàn)象，因此假設(shè)無(wú)限小的時(shí)間dt中，剪切應(yīng)變?cè)隽繛椋?/p>

(10)

圖1 搜索示意圖

圖1中的縱坐標(biāo)表示坡高。將層數(shù)相同的單元合并為一組，利用單元ID標(biāo)記每個(gè)單元。完成剪切滑點(diǎn)的選取工作后，選取受拉滑裂點(diǎn)。環(huán)境動(dòng)態(tài)變化下，邊坡受到水平慣性力，坡頂位置由于受拉而出現(xiàn)破裂(丁鑫品等，2020)。采用簡(jiǎn)化方法假設(shè)受拉裂縫為直向的細(xì)長(zhǎng)裂縫，當(dāng)坡頂出現(xiàn)連續(xù)受拉塑性單元時(shí)，疊加所有單元獲得高度值，并將該值記為裂縫最大深度。受拉塑性單元與剪切塑性單元存在交互關(guān)系，針對(duì)這一研究結(jié)果定義受拉塑性單元高度為受拉裂縫深度，圖2為受拉裂縫深度示意圖。

圖2 受拉裂縫深度示意圖

圖2中縱坐標(biāo)表示裂縫深度，黑色區(qū)域表示拉塑性單元區(qū)，h表示受拉塑性單元與剪切塑性單元之間的深度(m)，也就是受拉裂縫深度，豎直提取該值以上單元坐標(biāo)，作為受拉滑裂關(guān)鍵單元坐標(biāo)。綜合上述兩組關(guān)鍵點(diǎn)，搜索邊坡最不利滑裂面。

1.3 GA-Sarma算法搜索邊坡滑裂面

GA-Sarma算法能夠解決以折線形為滑面形態(tài)、滑裂路徑可追蹤、坡向不連續(xù)的巖質(zhì)邊坡最危險(xiǎn)滑面的全局最優(yōu)問(wèn)題(張劍波等，2013；王紫杰等，2020)，其將滑裂面看作折線形，利用Sarma極限平衡法確定穩(wěn)定性系數(shù)，通過(guò)遺傳算法優(yōu)化搜索最不利滑裂面的位置。很明顯，GA-Sarma算法是具有搜索全局最優(yōu)問(wèn)題能力的遺傳算法和傳統(tǒng)力學(xué)模型的融合，為快速、合理得搜索邊坡滑裂面提供了技術(shù)支持(毛謙等，2008；林永生和陳勝宏，2013；麻官亮等，2016)。根據(jù)預(yù)先模擬的邊坡模型，該算法首先建立目標(biāo)函數(shù)。當(dāng)滑裂面由隨機(jī)一點(diǎn)向坡頂擴(kuò)展時(shí)，可能存在無(wú)數(shù)條變化路徑，因此假設(shè)滑裂路徑的方向?yàn)棣?。?dāng)坡體內(nèi)存在順坡向不連續(xù)結(jié)構(gòu)面時(shí)，則路徑沿不連續(xù)結(jié)構(gòu)面變化(王章瓊等，2018)。因此根據(jù)Sarma算法得到：

K=g(φn,γn)

(11)

式中：φn、γn表示方向?yàn)棣聲r(shí)，選取的剪切滑裂點(diǎn)和受拉滑裂點(diǎn)。決策變量就是上述公式中參數(shù)的數(shù)量，因此將函數(shù)優(yōu)化問(wèn)題轉(zhuǎn)換成算法搜索問(wèn)題。根據(jù)其他方法中適應(yīng)度函數(shù)的使用方法可知，該方法可以度量個(gè)體最優(yōu)解。將求解穩(wěn)定性系數(shù)作為滑裂面全局搜索的一個(gè)最小值問(wèn)題，則得到適應(yīng)度函數(shù)的計(jì)算公式：

h(φn,γn)=max(g(φn,γn))-g(φn,γn)

(12)

公式(12)的物理意義，描述了穩(wěn)定性系數(shù)最小的路徑適應(yīng)度最大，當(dāng)搜索過(guò)程中出現(xiàn)變異問(wèn)題時(shí)，有更大機(jī)率被保存下來(lái)。算法利用選擇算子評(píng)價(jià)個(gè)體，根據(jù)得到的適應(yīng)度評(píng)價(jià)結(jié)果，強(qiáng)化算法的收斂性能(竇勇，2019)。已知隨機(jī)采樣方式下，適應(yīng)度實(shí)際值與個(gè)體選中概率呈正相關(guān)，因此假設(shè)N表示整個(gè)搜索群體集合，用i表示個(gè)體排序情況，通過(guò)上述公式(12)獲得，則概率結(jié)果為：

(13)

式(13)中：i∈N表示單個(gè)個(gè)體和個(gè)體總數(shù)量。交叉算子作為有效計(jì)算工具，能夠產(chǎn)生新個(gè)體，提高整個(gè)群體中所有個(gè)體的多樣性和動(dòng)態(tài)變化特性，防止部分個(gè)體結(jié)果，提前成熟而影響最終的位置搜索(張子映等，2018)。為了防止破壞大量現(xiàn)有的較好模式，要求變異概率Pi的值要滿足最小值要求，則變異操作可按照?qǐng)D3所示的方法執(zhí)行。

圖3 變異操作示意圖

為擺脫個(gè)體多樣性喪失、搜索進(jìn)程陷入僵局的問(wèn)題，GA-Sarma算法通過(guò)多次增大變異概率Pi的值改善這些現(xiàn)象，同時(shí)為了進(jìn)一步加強(qiáng)搜索結(jié)果，還引入了災(zāi)變策略(CS)，大規(guī)模消減處理群體，并進(jìn)行新后代生成操作，強(qiáng)化搜索最優(yōu)這一目的。至此基于GA-Sarma算法，實(shí)現(xiàn)邊坡最不利滑裂面搜索方法。

2 實(shí)驗(yàn)與分析

為驗(yàn)證上述設(shè)計(jì)的基于GA-Sarma算法的邊坡最不利滑裂面搜索方法的應(yīng)用性能，設(shè)計(jì)如下實(shí)驗(yàn)。

2.1 實(shí)驗(yàn)準(zhǔn)備

通過(guò)地形測(cè)量技術(shù)、地質(zhì)測(cè)量技術(shù)、物探技術(shù)、鉆探技術(shù)以及巖礦分析等技術(shù)手段，對(duì)西北某地區(qū)中存在潛在滑裂面的峽谷、山坡進(jìn)行勘察。此次勘察工作依據(jù)8項(xiàng)不同的城市治理規(guī)范、土地使用規(guī)范、滑坡防治規(guī)范以及工程勘察方案而進(jìn)行。

實(shí)驗(yàn)的工程測(cè)量工作中，采用2000國(guó)家大地坐標(biāo)系，建立仿真環(huán)境測(cè)試模型，高模型的基準(zhǔn)為1985國(guó)家高程基準(zhǔn)，設(shè)置比例尺為1:500，范圍為滑坡體及周?chē)獾?0 m范圍內(nèi)，測(cè)量面積0.2 km2。利用無(wú)人機(jī)對(duì)整個(gè)邊坡進(jìn)行掃描，單次掃描面積設(shè)置為0.30 km2。采用高密度電法進(jìn)行區(qū)域物探。在選擇的勘測(cè)區(qū)域內(nèi)設(shè)置260個(gè)物理點(diǎn)、15個(gè)質(zhì)檢點(diǎn)。計(jì)算各點(diǎn)均方相對(duì)誤差，計(jì)算公式為：

(14)

圖4 Ⅱ區(qū)、Ⅲ區(qū)山體滑動(dòng)裂縫圖(部分)

根據(jù)圖4顯示的實(shí)際調(diào)查結(jié)果可知，在實(shí)驗(yàn)選擇的勘測(cè)區(qū)域中，存在大量威脅性較大的滑裂面以及潛在滑裂面。為此，使用不同的勘察手段分析測(cè)試區(qū)域內(nèi)的巖土類(lèi)型。由于篇幅的限制，本研究主要對(duì)Ⅲ區(qū)展開(kāi)研究。表1為獲得的Ⅲ區(qū)各個(gè)巖土層力學(xué)指標(biāo)。

表1 測(cè)試區(qū)域巖土層主要力學(xué)指標(biāo)勘測(cè)結(jié)果

為了檢驗(yàn)本文方法的實(shí)際應(yīng)用性能，設(shè)計(jì)對(duì)比應(yīng)用實(shí)驗(yàn)。將本文方法作為實(shí)驗(yàn)組，將兩組傳統(tǒng)搜索方法-Bishop方法及基于SLFT(滑移線場(chǎng)理論)設(shè)計(jì)邊坡雙折減系數(shù)強(qiáng)度的失穩(wěn)判據(jù)，然后根據(jù)極限狀態(tài)下的邊坡數(shù)據(jù)分析結(jié)果，討論邊坡局部位置是否開(kāi)始出現(xiàn)裂紋方法(分別是方宏偉等(2019)和呈志軒等(2019)所得的方法)分別作為對(duì)照A組、對(duì)照B組，結(jié)合圖4與表1顯示的數(shù)據(jù)，對(duì)整個(gè)勘察地區(qū)進(jìn)行邊坡最不利滑面搜索工作，通過(guò)比較搜索結(jié)果是否最優(yōu)，來(lái)確定不同方法之間的差異性。

2.2 最優(yōu)解測(cè)試結(jié)果分析

利用仿真軟件、CAD等軟件模擬整個(gè)搜索區(qū)域的三維影像，已知仿真場(chǎng)景中存在3處最不利滑裂面、5處剛剛形成且危險(xiǎn)程度偏低的滑裂面以及14處潛在滑裂面。因此，分別以最不利滑裂面、剛剛形成且危險(xiǎn)程度偏低的滑裂面(一般滑裂面)和潛在滑裂面為搜索對(duì)象，利用不同方法對(duì)其展開(kāi)搜索，根據(jù)搜索結(jié)果驗(yàn)證不同方法的實(shí)際性能，其比較準(zhǔn)則為不同方法對(duì)仿真場(chǎng)景中的不同滑裂面搜索的吻合度和數(shù)量。圖5為三組不同方法的最優(yōu)滑裂面搜索結(jié)果。

圖5 最優(yōu)搜索結(jié)果對(duì)比

根據(jù)圖5顯示的測(cè)試結(jié)果可知，面對(duì)同樣的搜索場(chǎng)景，實(shí)驗(yàn)組獲得了所有最不利滑裂面的所在位置。對(duì)照組A遺漏了1處最不利滑裂面所在位置，同時(shí)將3處剛剛形成且危險(xiǎn)程度偏低的滑裂面，當(dāng)成了最不利滑裂面。對(duì)照B組只獲得了一組最不利滑裂面所在位置，剩余的2處最不利滑裂面直至實(shí)驗(yàn)測(cè)試結(jié)束，也沒(méi)有被發(fā)現(xiàn)，同時(shí)將4處剛剛形成且危險(xiǎn)程度偏低的滑裂面，當(dāng)成了最不利滑裂面。綜合上述實(shí)驗(yàn)測(cè)試結(jié)果可以看出，兩組常規(guī)邊坡最不利滑裂面搜索方法中，使用的算法性能較差，過(guò)早收斂導(dǎo)致獲得的結(jié)果不是最優(yōu)。

在此基礎(chǔ)上，為了進(jìn)一步驗(yàn)證不同方法的應(yīng)用性能、保證實(shí)驗(yàn)測(cè)試結(jié)果的可靠性，保證基本測(cè)試條件不變，利用3種方法對(duì)選擇的測(cè)試場(chǎng)地中的最不利滑裂面展開(kāi)20輪搜索，以此來(lái)更直觀地觀察不同方法的搜索效果。表2為20輪測(cè)試過(guò)程中得到的最優(yōu)解測(cè)試結(jié)果。

表2 最不利滑裂面搜索最優(yōu)解統(tǒng)計(jì)表

表2匯總結(jié)束后，發(fā)現(xiàn)對(duì)照B組從第11輪開(kāi)始，一直到第17輪的測(cè)試結(jié)果均為0。為了保證實(shí)驗(yàn)測(cè)試數(shù)據(jù)完全可信，利用事先安裝的監(jiān)控軟件，查詢對(duì)照B組儀器設(shè)備以及軟件的運(yùn)行狀態(tài)，發(fā)現(xiàn)對(duì)照B組正常完成了20輪搜索工作，并不存在任務(wù)系統(tǒng)或軟硬件故障問(wèn)題，因此第11輪至第17輪的測(cè)試結(jié)果可信。根據(jù)表2中的數(shù)據(jù)可以直觀看到，實(shí)驗(yàn)組方法在20輪測(cè)試中，全部得到了3組最不利滑裂面；對(duì)照A組每次搜索最不利滑裂面時(shí)，總是存在遺漏位置的情況；對(duì)照B組則有35%的概率，難以獲得最不利滑裂面的搜索結(jié)果。

為了方便比較，計(jì)算三組方法20輪測(cè)試的最優(yōu)解平均值，分別為3處、1.4處以及0.8處。根據(jù)該計(jì)算結(jié)果可以得出：在同樣的測(cè)試背景環(huán)境下，本文方法完全能夠獲得最優(yōu)解，而兩組傳統(tǒng)方法獲得最優(yōu)解的平均值不足2處和1處，即面對(duì)3處最不利滑裂面所在位置，兩組方法每次搜索時(shí)只能獲得1.4處以及0.8處最不利滑裂面。由此可見(jiàn)，本文設(shè)計(jì)的基于GA-Sarma算法的邊坡最不利滑裂面搜索方法更加滿足設(shè)計(jì)預(yù)期。

3 結(jié)論

針對(duì)實(shí)際實(shí)驗(yàn)結(jié)果，本研究在選取潛在滑裂面關(guān)鍵點(diǎn)的基礎(chǔ)上，引入了GA-Sarma算法，通過(guò)加強(qiáng)算法自身的計(jì)算能力、降低其他因子對(duì)算法的影響，并建立目標(biāo)函數(shù)和適應(yīng)度函數(shù)，在考慮異變的前提下搜索邊坡最不利滑裂面。具體結(jié)論如下：

(1)工程測(cè)量工作中，采用2000國(guó)家大地坐標(biāo)系，建立仿真環(huán)境測(cè)試模型，利用無(wú)人機(jī)對(duì)整個(gè)邊坡進(jìn)行掃描且采用高密度電法進(jìn)行區(qū)域物探。勘測(cè)區(qū)域質(zhì)量檢查總均方相對(duì)誤差在±3.42%之間，滿足上述8項(xiàng)規(guī)范中不大于±5%的要求。

(2)本文引入GA-Sarma算法中，利用遺傳算法對(duì)邊坡最不利滑裂面進(jìn)行搜索，顯示出強(qiáng)大的全局最優(yōu)搜索能力。遺傳算法和Sarma算法結(jié)合是具有搜索全局最優(yōu)問(wèn)題能力的遺傳算法和傳統(tǒng)力學(xué)模型的融合，為快速、可靠得搜索邊坡滑裂面提供了技術(shù)支持。

(3)在20輪測(cè)試中，對(duì)照A組每次搜索最不利滑裂面時(shí)，總是存在遺漏位置的情況；對(duì)照B組則有35%的概率難以搜索到最不利滑裂面；實(shí)驗(yàn)組方法全部得到了3個(gè)最不利滑裂面，取得了最好的搜索效果。

此外，在今后研究工作中，基于GA-Sarma算法進(jìn)一步擴(kuò)大搜索范圍，對(duì)多種不同地質(zhì)進(jìn)行系統(tǒng)研究，優(yōu)化特征提取或智能化處理步驟是值得研究的課題，為土地保護(hù)與治理提供更加快速、可靠的技術(shù)支持。