基于智能算法的水工隧洞圍巖穩定可靠度分析

0 引言

在實際工程和理論研究中常將圍巖的變形（或應變）作為判定其是否穩定的重要指標[1]。在考慮巖土力學參數不確定性的基礎上，以水工隧洞圍巖的變形（或應變）建立可靠度計算功能函數并進行可靠度分析，可為水工隧洞圍巖的安全穩定評判提供重要依據。

收斂約束法是考慮圍巖與支護相互作用下計算隧洞圍巖變形的實用設計方法之一[2]。但該方法中圍巖的變形隱含在關于巖土力學參數的非線性方程組中，這使得以隧洞圍巖的變形（或應變）建立的可靠度計算功能函數為隱式非線性函數。傳統的可靠度計算方法在處理隱式非線性功能函數時，或表現為計算效率低下（如蒙特卡洛法），或因無法求解關于隨機變量的偏導數而不適用（如一次二階矩法），或計算精度不高（如以二次多項式構建響應面的傳統響應面法），因而需探索更高效準確的可靠度計算方法。

可靠指標是衡量可靠性的數量指標，根據其幾何含義，可靠度計算問題可轉化為約束優化問題進行求解[3]。對于實際工程往往需要建立數值模型來求解可靠度計算的功能函數，而一次數值模型的計算將耗費大量的計算機時。為此，大量元模型方法被引人到可靠度分析領域[4]，這些方法以少量的真實樣本建立訓練集并構造元模型，以替代原本復雜的可靠度計算功能函數，從而減少真實功能函數的調用次數。高斯過程回歸[5]是一種監督學習算法，因其具有建模過程易實現、超參數可自適應獲得等優點，已被廣泛應用于可靠度分析領域[6]。本文選用該算法建立元模型以代替真實功能函數，并將其作為約束優化問題中的約束函數。智能優化算法相比于傳統優化算法，可克服其易陷入局部最優解、要求目標函數必連續可微等缺點，從而更適合于隱式非線性目標函數的優化求解[7]。成吉思汗鯊魚優化器[8]是一種較新的群體智能優化算法，具有預設參數少、優化效率高的優點，可作為可靠度計算中約束優化問題的尋優方法。綜上，本文提出高斯過程回歸-成吉思汗鯊魚優化器智能算法，并將其應用于由收斂約束法構建的水工隧洞圍巖穩定可靠度計算問題中，以期為水工隧洞圍巖穩定分析提供一種精確、高效的可靠指標計算方法。

1 收斂約束法

收斂約束法常用于隧洞（道）及其支護的初步設計階段。該方法是一種理論分析方法，其最初的提出基于一系列基本假定[9]。但隨著理論研究的深入和數值計算技術的發展，該方法逐漸打破了最初的一些假定，將更多現實因素考慮在內，且與數值計算結合得更加緊密，以更好地應用于實際工程。目前，收斂約束法仍是隧洞（道）初步設計的主流方法之一。

收斂約束法包含3條基本曲線，見圖1，分別是縱斷面變形曲線（LDP）、圍巖特性曲線（GRC）和支護特性曲線（SCC）。這3條曲線簡單明了地解釋了隧洞（道）圍巖與支護間的相互作用，并可獲得圍巖的最終變形。

以圍巖應變建立可靠度計算功能函數，如式（1）所示，其中圍巖變形由圍巖特性曲線和支護特性曲線的交點求得，即圍巖變形的求解隱含在圍巖特性曲線和支護特性曲線的非線性方程組中，因此功能函數為隱式非線性函數。

式中： ε_max 是圍巖允許最大應變； u 是圍巖徑向變形； r₀

是隧洞半徑。

2基于高斯過程回歸－成吉思汗鯊魚優化器響應面法的圍巖穩定可靠度計算方法

2.1 可靠指標幾何含義

在標準正態空間中，設計驗算點 ）是位于極限狀態面 g（X）=0 上且距離坐標原點最近的點，該最短距離即為可靠指標 β 。根據此幾何意義，可靠指標 β 的求解可轉化成式（2）所示約束優化模型：

式中： 為隨機變量 X 的協方差矩陣。

則最終尋優結果點 P（X₁，X₂，…，X_n） 即為設計驗算點 P^* ，點 P^* 到原點的距離即為可靠指標 β 。

2.2 高斯過程回歸

高斯過程回歸是一種適用于處理小樣本、非線性回歸問題的基于貝葉斯理論和統計學習理論的非參數機器學習算法[10]。它基于高斯過程理論，通過學習訓練集的輸人-輸出關系，來估計未知函數的概率分布及給出對預測結果的不確定性估計。

假定訓練集為 ，其中 為第 i 個多維輸入向量， y_i 為對應的第 i 個輸出值。高斯過程可完全由均值函數 和協方差函數 k（x， x^′ ）確定，即：

基于貝葉斯理論，由訓練集輸出 y 和預測樣本 x^* 的預測輸出 y^* 的聯合先驗分布，可得預測值 y^* 的后驗分布為

其中，

式中： 是協方差矩陣；高斯白噪聲 ε～N（0，σ_n²）;I 為單位矩陣。

據此，當給定訓練集 T 及預測樣本 x^* 時，其預測輸出 y^* 就可根據上述理論獲得，

為減少真實功能函數的調用次數，本文用高斯過程回歸元模型g代替真實功能函數，即將式（2）中的約束方程 替換為

2.3 成吉思汗鯊魚優化器

成吉思汗鯊魚優化器可模擬成吉思汗鯊魚的捕食行為，是一種較新的群體智能優化算法[8]。該算法模擬了成吉思汗鯊魚捕食行為的4個階段，分別是狩獵階段、移動階段、覓食階段和自我保護階段。其中自我保護階段是一個特殊的階段，可幫助算法跳脫局部最優。成吉思汗鯊魚優化器在一次迭代中依次執行上述4種行為，以不斷更新種群個體的最優位置，直到達到迭代停止條件為止。具體計算公式可參見文獻[8]。

相比于其他優化算法，該優化器穩定性強，效率高，且除需設定迭代次數和種群個體數外，只需再另設參數 m 即可完成所有參數的設定。而已有研究表明，m =1.5時算法取得最優結果[8]

2.4 可靠度計算方法流程

根據可靠指標的幾何含義，可將水工隧洞圍巖穩定可靠度計算轉化為約束優化問題的求解。目標函數為可靠指標 β 的數學表達式，即式（2）中的第一式；約束方程為近似極限狀態面，即式（7），且該近似極限狀態面會隨著訓練集的更新而不斷更新；優化方法為成吉思汗鯊魚優化器。具體執行步驟為：

（1）用拉丁超立方抽樣選取 5d（d 為問題維度）個初始樣本，并求得其真實功能函數值，從而構成初始訓練集，記所需樣本量為 N_noE=5d 。

（2）基于訓練集，構建高斯過程回歸元模型。

（3）將高斯過程回歸元模型代替式（2）中的等式約束條件（真實功能函數），并用外點罰函數法將等式約束轉化為無約束優化問題，從而用成吉思汗鯊魚優化器進行優化求解。為減少波動性，重復使用5次成吉思汗鯊魚優化器進行尋優，5次結果中的最優解即為新得到的計算樣本。

（4）將新得到的計算樣本代人真實功能函數，得到其真實功能函數值，并加入到訓練集中以更新訓練集， N_DOE=N_DOE+1 。

（5）重復步驟（2）＼～（4），直至達到收斂標準，即lβnow-βbefore|/βbefore lt;10-3

本文所提方法用高斯過程回歸元模型替代水工隧洞圍巖穩定可靠度計算功能函數，并將其融入到成吉思汗鯊魚優化器的算法框架內，計算流程圖如圖2所示。

3 算例分析

如圖3所示典型城門洞型地下洞室，巖體為均質理想彈塑性體，服從Mohr-Coulomb屈服準則。巖體容重 γ=25kN/m³ ，黏聚力 c=0.5MPa ，內摩擦角 φ= 38^° 。彈性模量 E 、泊松比 u 和均勻分布的地表荷載 p 為隨機變量，且 ） ，ν～

LN（0.3，0.015）_?p～LN（0.6MPa，0.12MPa） 。確定性分析下，圍巖沉降位移等值線圖見圖3。其中洞頂沉降位移為 16.1mm ，與文獻[11]中所得 16.2mm 極為接近。以洞頂沉降位移為分析對象，功能函數為式（8）。

式中： Z_max=25mm 為洞頂沉降極限位移值。

運用本文所提智能算法進行可靠指標的計算。首先用拉丁超立方抽樣選取15（ d=3 ）個初始樣本，代入模型計算其真實響應值以組成初始訓練集，從而構建初始高斯過程回歸元模型。接著，用該元模型代替真實功能函數，并根據可靠指標的幾何含義，將可靠度問題轉化為優化問題。然后用成吉思汗鯊魚優化器進行優化求解，且為了減少波動性，用成吉思汗鯊魚優化器求解5次并選出最優解，該最優解即為新選取的樣本點。之后將新樣本點及其真實響應值加入到訓練集中，以修正高斯過程回歸元模型，再用成吉思汗鯊魚優化器進行優化求解，重復此步驟，直到計算所得可靠指標收斂為止。可靠指標計算結果見表1。

由表1可看出，本文方法適用于隱式功能函數的求解。以蒙特卡洛法計算結果為基準解，本文方法相比傳統響應面法和蒙特卡洛法，真實功能函數調用次數更少，即較少調用數值軟件進行計算。由于運行一次成吉思汗鯊魚優化器需消耗約35s（電腦配置為：Intel（R）Core（TM）i7-9700 CPU @ 3.00GHz ，RAM為8GB，下同），即尋找新樣本點所需時間約 175s ，且該算例模型簡單，一次模型調用及計算結果返回需約6s，因此本文所提方法的計算耗時大于傳統響應面法。但傳統響應面法分析復雜功能函數時計算精度往往難以保證，而高斯過程回歸具有很強的適應性[10]；且對于復雜工程結構，消耗約 3min 的時間尋找新樣本點以參與計算，比大量運行一次就需幾十分鐘甚至幾小時的數值模型更為節約計算機時。因此，本文方法對于大型復雜工程結構的可靠度分析，能很好地達到計算效率和準確率之間的平衡。

4 工程實例

位于雅襲江干流上的錦屏二級引水式地下電站，共有4條引水隧洞。隧洞穿越綠泥石片巖地層時，遭遇了軟巖大變形。為此，本文以穿越綠泥石片巖地層的洞段為研究對象，并以圍巖應變建立可靠度計算功能函數，計算其圍巖穩定可靠指標。

根據該洞段所受初始地應力情況[12]及開挖尺寸，建立三維數值模型如圖4所示，

為簡化計算，假設隧洞一次開挖成型，開挖進尺為1m 。開挖巖體為理想彈塑性體，服從Mohr-Coulomb屈服準則。隨機變量統計特征見表2。

用收斂約束法求解拱頂變形。結合數值模型及式（9）[13]，計算圍巖的縱斷面變形曲線。采用應力釋放

法計算圍巖特性曲線，并用式（10）對模型計算結果進行擬合。將支護看作理想彈塑性材料，支護特性曲線見式（11），初期支護參數見表3[14]

式中： u 為圍巖變形； 是圍巖的最大徑向變形； λ 是位移完成率； λ₀=u₀∣u_r^max∣，u₀ 是掌子面達到監測斷面時監測點圍巖的變形； M_?1=λ₀/（1-λ₀）M，M 和 M_?1 是未知參數，可通過Levenberg-Marguard和通用全局優化算法得到[15]； L_s 是支護安裝距掌子面的距離。

式中： p_i 為圍巖壓力； k，f，a，b，c，d 是未知參數，可通過Levenberg-Marguard和通用全局優化算法得到； p_i^cr 是彈塑性分界區臨界圍巖壓力。

式中： p_s 為支護力； k_tot 為總的支護剛度； u_r^D 為支護（或圍巖）的最終變形； u_rⁱⁿ 為支護安裝前圍巖先期變形；p_s^max 為支護的最大承載力。 k_tot 和 p_s^max 可由支護參數求得[14] O

性曲線如圖5所示。根據收斂約束法，假設在距掌子面 1m 后安裝初期支護，由縱斷面變形曲線可得圍巖在支護安裝前的先期變形為 0.11m 。由圍巖特性曲線和支護特性曲線的交點可得，圍巖的最終變形為0.21m ，應變為 ε=u/r₀=0.21/6.90=3.04% 。由于施工期可按收斂應變的大小對軟巖變形程度進行初步預測評估[16] .2.5%lt;ε?5.0% 屬于中等擠壓變形，穩定性差，圍巖自穩時間很短。此外，本洞段采用 4% 的應變來判定圍巖是否會產生流變變形[17]， 3.04% lt;4% ，則在確定性分析下，拱頂圍巖為中等擠壓變形，且不會產生流變變形。

在考慮巖土力學參數不確定性的情況下，以圍巖應變建立功能函數，即式（1）， ε_max 取不同值時，可靠度計算結果見表4。由表4可知，本文方法在較少調用數值模型（真實功能函數）的情況下即可得到計算結果。由于一次數值分析需消耗 5.5h 左右，因此本文方法大大減少了計算機時。計算所得可靠指標隨著ε_max 取值的增大而增大，失效概率則隨之而減少。當（20 ε_max=2.5% 時，拱頂圍巖的失效概率為 83.81% ，即拱頂產生中等及以上程度擠壓變形的概率為 83.81% ;當 ε_max=5% 時，拱頂圍巖的失效概率為 3.92% ，即拱頂產生嚴重及以上程度擠壓變形的概率為 3.92% （2 （5.0%lt;ε?10.0% 為嚴重擠壓變形[16]）。因此，可靠度分析可為拱頂圍巖的擠壓變形程度賦予一定的概率意義。此外，當 ε_max=4% 時，即從圍巖流變變形角度分析，拱頂圍巖的失效概率為 14.46% ，這在工程上通常不可接受，即可靠度分析結果認為拱頂是不安全的；而確定性分析結果表明，拱頂圍巖不會產生流變變形。可見，可靠度分析結果與確定性分析結果不一致，可靠度分析結果更為保守。

此外，本文所提方法還可求出驗算點。由于標準正態空間中的驗算點可反映功能函數對各變量的敏感程度[18]，由表4可知，對于該工程，當 ε_max=2.5% 時，功能函數對參數 φ 最為敏感，其次為參數 E 和 ∣c∣ ；當ε_max=4% ？ 5% 時，功能函數對參數 E 和 c 最為敏感，對參數 φ 最不敏感。可見，當 ε_max 取不同值時，功能函數對不同變量的敏感性會有所改變。因此， ε_max 的取值極為重要，它不僅會決定可靠指標（或失效概率）的取值從而影響圍巖安全穩定性的評估結果，還會影響功能函數對不同變量的敏感程度，從而在一定程度上影響隧洞及其支護的設計和施工方案。

綜上，本文方法不僅適用于隱式非線性功能函數可靠度的求解，從可靠性的角度來分析和指導工程的建設施工；還可求得驗算點信息，以分析各因素對功能函數的敏感性，從而對重點影響因素進行勘測、采集、整理、分析，以獲得更準確的數據信息，更好地服務于研究對象。

5結論

本文基于收斂約束法，在考慮巖土力學參數不確定性的基礎上，以水工隧洞圍巖應變建立可靠度計算功能函數。由于功能函數為隱式非線性函數，本文根據可靠指標的幾何含義，將水工隧洞圍巖穩定可靠度計算轉化為約束優化問題的求解，并提出了基于高斯過程回歸-成吉思汗鯊魚優化器的水工隧洞圍巖穩定可靠度智能計算方法。通過算例和工程實例的計算分析，得到如下結論：

（1）算例分析表明，本文方法適用于隱式非線性功能函數的求解，且相比于傳統響應面法，具有更高的計算效率。本文方法與蒙特卡洛法的計算結果精度相當，計算精度可滿足工程需要。

（2）工程實例分析表明，本文方法只需調用少量數值模型（真實功能函數）即可獲得可靠指標及失效概率，計算效率高；且可以求得驗算點，分析各因素對所研究功能函數的敏感性，從而對重點影響因素進行重點分析，以更好地服務于實際工程的設計和安全穩定性分析。

（3）根據工程實例的確定性分析結果，在施工期，拱頂圍巖發生中等擠壓變形，且不會產生流變變形。

而根據可靠度分析結果，拱頂圍巖產生中等及以上程度的擠壓變形概率為 83.81% ，產生嚴重及以上程度擠壓變形的概率為 3.92% ，且有 14.46% 的概率會產生流變變形。可見，可靠度分析能在考慮不確定性因素的前提下對確定性分析結果賦予相應的概率。對于本工程實例中拱頂圍巖的流變變形穩定性分析，可靠度分析結果與確定性分析結果不一致，可靠度分析結果更為保守。因此，可靠度分析應與確定性分析結合使用，以為實際工程的安全穩定性提供更為全面和綜合的評價，從而為工程的安全建設提供指導和建議。

參考文獻：

[1]侯欽禮，張雨霆，孫海清.弱膠結軟巖水工隧洞力學參數反演及穩 定性評價[J].人民長江，2023，54（增2）：202-207.

[2]LV Q，CHAN C L，LOW B K. Probabilistic evaluation of ground - support interaction for deep rock excavation using artificial neural network and uniform design[J]. Tunnelling and Underground Space Technology，2012，32：1-18.

[3]王剛，秦凈凈，管莉莉.基于遺傳算法的重力壩多滑面穩定可靠度 分析[J].巖石力學與工程學報，2016，35（增1）：3153-3161.

[4]姜廣倫，孔存芝，南驍聰，等.土質邊坡位移概率反分析與失穩概 率預測[J].人民長江，2022，53（9）：155-162.

[5]RASMUSSEN C E，WILLIAMS C K I. Gaussian processes for machine learning[M].Cambridge，MA：MIT Press，2005.

[6]LI M，WANG G，QIAN L，et al. ASS-GPR： adaptive sequential sampling method based on gaussian process regression for reliability analysis of complex geotechnical engineering[J].International Journal of Geomechanics，2021，21（10） ：04021192.

[7］王彪龍，劉曉，郭將.基于雜交粒子群響應面的邊坡穩定可靠性算 法［J].人民長江，2018，49（16）：97-105.

[8]HU G，GUO Y，WEI G，et al. Genghis Khan shark optimizer：a novel nature- inspired algorithm for engineering optimization[J].Advanced Engineering Informatics，2023，58：102210.

[9]ORESTE P.The convergence - confinement method ： roles and limits in modern geomechanical tunnel design[J].American Journal of Applied Sciences，2009，6（4） ：757 -771.

[10］肖義龍，蘇國韶.結構可靠度分析的高斯過程響應面方法[J].人 民長江，2016，47（23）：82-85.

[11]LI D Q，JIANG S H，CHEN Y F，et al. Reliability analysis ofserviceability performance for an underground cavern using a non- intrusive stochastic method[J]. Environmental Earth Sciences，2014，71（3）： 1169 -1182.

[12]呂品.錦屏水電站綠片巖段擴挖及落底開挖穩定研究[D].大連： 大連理工大學，2011.

[13]ZHAO D，JIA L，WANG M，et al. Displacement prediction of tunnels based ona generalised Kelvin constitutivemodel and its applicationin a subsea tunnel[J]. Tunnelling and Underground Space Technology， 2016，54：29 -36.

[14]ORESTE P P.Analysis of structural interaction in tunnels using the convergence -confinement approach[J]. Tunnelling and Underground Space Technology，2003，18（4）：347-363.

「15]張妍珺，蘇凱，周利，等.基于收斂-約束法的隧洞縱向變形演化規律研究與支護時機估算[J].巖土力學，2017，38（增1）：471-478.

[16]彭土標.水力發電工程地質手冊［M].北京：中國水利水電出版社，2011.

[17］張春生.深埋隧洞巖石力學問題與實踐[M].北京：中國水利水電出版社，2016.

[18]張明.結構可靠度分析：方法與程序[M].北京：科學出版社，2009.

（編輯：鄭毅）

Stability reliability analysis of surrounding rock of hydraulic tunnels based on intelligent algorithm

LI Mengyao，WANG Gang，MA Zhenyue，KANG Fei （School of Infrastructure Engineering，Dalian University of Technology，Dalian 116O24，China）

Abstract：Theconvergence-confinement methodismostlyusedinthepreliminarydesignstageof tunnels，buttheformulasfor thedeformationofthesuroundingockinthismethdareimplicitnonlinarfuctions，hichmakesitdiffcult tosolvethmbyusingtraditionalreliabilitycalculationmethod.Iviewofthis，thestabilityreliabilitycalculationofthesuroundingockofhdau lictunnelswastransformedintoaconstrainedoptimizationproblembasedonthegeometric meaningofthereliabilityindex.Then， theGaussianprocessregressionwasusedtofttherealperformancefunction，andtheGenghisKhanSharkOptimizerwasusedto solvetheoptimizationproblem.Asaresult，newmethdforcalculatingthestabityreliabilityofthesuroundingrockofhydraulic tuels，basedontheGauianprocessregresionandtheGenghisKhan SharkOptimizer，asproposed.ThenalysisofacalculationexapleshowedtatteproposedmethdhasstrongplabilityighomputationalfiencyndacuacyItealculation，thereliabilityindexandfailureprobabilitycouldbeobtainedbyafewtimesofoperatingthenumericalmodel，andatthe sametime，thedesign points were given，whichcouldbeusedtofurtheranalyzethesensitivityofperformancefunctiontoeach factor.Theexampleofapracticalhydraulictunnelshowedthattheprobabilityof producing moderateorhigherdegreeofextrusion deformation at the tunnel crown was 83.81% and the probability of producing rheological deformation was 14.46% . Compared withthedeterministicanalysis，thereliabilityanalysiscangiveprobabilisticsignificancetothedeterministicanalysis.Withrgard totherheologicaldeformation，theresultofthereliabilityanalysiswasmoreconservativethanthatofthedeterministicanalysis.

Keywords：hydraulictunel;stabilityofsurroundingrock；reliabity；convergence-confinement method；Gausianprocess regression；Genghis Khan Shark Optimizer