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

李生清

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

0 引言

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

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

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

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

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

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

z=f(x,y)

(1)

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

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

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

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

(2)

(3)

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

(4)

(5)

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

(6)

(7)

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

G=G1/G2

(8)

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

(9)

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

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

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

(10)

圖1 搜索示意圖

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

圖2 受拉裂縫深度示意圖

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

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

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

K=g(φn,γn)

(11)

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

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

(12)

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

(13)

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

圖3 變異操作示意圖

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

2 實驗與分析

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

2.1 實驗準(zhǔn)備

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

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

(14)

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

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

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