基于模擬退火法尋優的公路邊坡抗剪強度參數反分析研究

張曉軍, 曹 娟

(1.湖南省交通科學研究院有限公司， 湖南 長沙 410015; 2.湖南金沙路橋建設有限公司， 湖南 長沙 410100)

0 引言

隨著我國公路建設事業的快速發展，公路技術人員積累了大量的工程建設資料和經驗，在此前提下，為削減勘探和土工試驗成本，對許多邊坡工程采用經驗法或者類比法來估算邊坡的抗剪強度參數[1]。此類方法盡管操作簡單方便，但受人員工程建設經驗、巖土體所處地域差異等不確定因素影響，存在一定局限性。反分析法是獲取邊坡巖土體參數的一種間接方法，即以現場監測信息為基礎，通過反分析模型或算法反推得到巖土體的實際材料參數[2]。

國外巖土工程反分析研究始于上世紀70年代，目前已得到了長足的應用與發展。Finno等[3]利用與誤差函數拓撲相關的梯度方法和遺傳算法，根據巖體原位測試數據反推Mohr-Coulomb模型的參數。Deb等[4]統計了大量巖體位移數據，利用BP神經網絡和模糊系統間接估計巖體的變形模量。Yazdani等[5]利用數值仿真計算替代現場監測，通過位移計算值對Siah Bisheh電站附近巖體參數、應變率以及節理參數進行反演，發現將數值計算結果作為反分析的依據具有可行性。國內方面，傳統的邊坡參數反分析往往基于剩余推力法[6]、極限分析法[7]和傳遞系數法[8]進行，而計算科學的進步促進反分析的手段逐漸向智能算法方向發展。張志增等[9]設計了正交試驗表，利用有限差分和人工神經網絡進行邊坡巖體力學參數的反分析。李南生等[10]為克服常規最優化方法在解決非線性問題中的缺陷，利用MATLAB進行了基于遺傳算法的邊坡反分析研究。張曉詠等[11]將ABAQUS有限元分析與強度折減法進行結合，對穩定滲流作用下的邊坡進行了抗剪切強度參數反分析。

實際邊坡工程中的反分析問題大多是非線性問題，若采用傳統優化方法存在速率低、無解性、多解性、穩定性差等缺陷。本文利用模擬退火法(SA)的優勢，將其引入到邊坡巖體抗剪強度參數反分析中，形成結合模擬退火法與多點位移優化進行反分析的方法，并以湖南省某高速公路2處邊坡為例，編寫了模擬退火法反分析程序，探討和驗證了該方法的可行性，研究結果可為邊坡分析和設計提供一定的參考。

1 模擬退火原理

在冶金學中，將材料加熱后，原子的能量會增大，導致原子容易離開原來的位置隨機移動而處于不穩定狀態；但溫度降低，原子的能量會減弱，移動能力會被削減，如果冷卻速率選擇恰當，原子就很有可能停留在能量更低的位置，這個過程就是退火[12]。如果將該原理引入尋優過程中，就形成了模擬退火法(SA)。在模擬退火法中，最優化問題的目標函數可被視為原子能量，為了獲取目標函數的全局極小值，可以先假定一個初始搜索模型，這就相當于退火中原子能量的初始值，采用合適的降溫速度，那么目標函數就有可能向降低的方向移動，使搜索參數逐漸收斂。

根據統計熱力學定律，分子在溫度T時，處于某種狀態Ei滿足Boltzmann率分布：

(1)

式中：T為絕對溫度；kb為玻爾茲曼常數；Z(T)為分配函數，即各個狀態相對幾率的總和。

物體從狀態i躍遷到狀態j的概率為：

(2)

如果，將模型參數向量m視為物體的某種狀態，將目標函數E(m)等效為物體的能量函數，利用控制參數T模擬物體的溫度，就可以得到Metropolis接受準則[13]：

P(mi→mj)=

(3)

從式(3)可以看出，當新模型的目標函數較小時，新模型100%將被接受，從而保證搜索向最優模型的方向移動；而新模型的目標函數較大時，也不會完全被拒絕，而是按照一定的概率接受該模型，這樣可以防止過快收斂而陷入局部極值，能夠有一定的概率跳出局部極值的陷阱，體現了全局化尋優的特點。

2 基于SA的邊坡巖體抗剪強度參數反演

在所研究邊坡巖體抗剪強度參數未知的情況下，利用模擬退火法(SA)進行參數反演的過程如下：

1)首先建立有限元模型，在計算能力滿足的情況下盡可能多設置幾組位移數據提取節點。對有限元模型輸入初始的抗剪強度參數c0和φ0，設定多組強度折減系數，進行有限元分析，然后提取各個節點在不同強度折減系數下的位移。

2)建立各個節點的位移-折減系數關系曲線，并進行函數擬合。

3)實際邊坡開挖后，在上述各個節點進行位移監測，獲取邊坡開挖后上述各個節點的監測位移穩定值。

4)根據計算數據和監測數據設定目標函數，編寫模擬退火法尋優程序，搜索獲取最優強度折減系數。

5)將初始抗剪強度參數按最優強度折減系數進行折減，得到反分析結果，進行誤差分析。

3 工程案例

3.1 工程概況

以湖南省某公路的2處邊坡為案例，對抗剪強度參數反演方法進行驗證。進行反分析的2處邊坡概況如下,現場如圖1所示。

a) A邊坡

b) B邊坡

1) K43+700～K43+800右邊坡(簡稱A邊坡)，3級坡，基巖為中風化砂質板巖，褐黃-青灰色，變余砂質結構，板狀構造，巖體較為完整，邊坡1級坡坡比為1∶0.75，2～3級坡坡比為1∶1。

2)K50+600～K50+675右邊坡(簡稱B邊坡)，2級坡，邊坡覆蓋層較薄，基巖為強風化砂質板巖，呈灰黃色，巖石結構構造可見，巖石呈砂狀及塊狀。邊坡1級坡坡比為1∶0.5，2級坡坡比為1∶0.75。

3.2 有限元模型分析

利用有限元分析軟件MIDAS NX建立了邊坡的初始有限元模型(見圖2)，考慮到初始強度參數要進行一系列折減，因此初始值的選取盡量大于真實強度參數取值，初始材料參數如表1所示，其中A邊坡巖體采用中風化板巖參數，B邊坡巖體采用強風化板巖參數。僅對粘聚力和內摩擦角進行折減，忽略坡面的框架防護結構，具體建模過程在此不做詳述。

對整個模型首先進行地應力平衡，然后去除開挖部分單元，再計算開挖后邊坡的位移。為了使擬合結果更好地反映邊坡整體變形的演變過程，在邊坡穩定的前提下進行了8次抗剪強度參數折減，然后選擇如圖2所示的5個節點提取位移數據，作為反分析的依據。

a) A邊坡

b) B邊坡

本文中強度折減法僅對模型上部的砂質板巖進行強度折減，強度折減系數與抗剪強度參數對應關系如表2、表3所示。

表1 初始材料參數材料重度/kN·m-3 泊松比彈性模量/MPa粘聚力/kPa內摩擦角/(°)強風化砂質板巖23.00.3580.0100.045.0中風化砂質板巖24.50.35200.0210.0100.0地基20.00.3860.080.035.0

表2 A邊坡在不同折減系數下的抗剪強度參數折減系數粘聚力/kPa內摩擦角/(°)1.4150.0071.431.6131.2562.501.8116.6755.562.0105.0050.002.295.4545.452.487.5041.672.680.7738.462.875.0035.71

表3 B邊坡在不同折減系數下的抗剪強度參數折減系數粘聚力/kPa內摩擦角/(°)1.190.9140.911.376.9234.621.566.6730.001.758.8226.471.952.6323.682.147.6221.432.343.4819.572.540.0018.00

3.3 位移-折減系數曲線

根據表2和表3中的數據設置有限元模型材料參數，對2個邊坡分別進行8組折減系數下的有限元分析，提取5個節點(如圖2所示)監測開挖后的位移，得到節點位移-強度折減系數曲線如圖3和圖4所示。

從圖3和圖4可以看出，各個節點的位移隨著折減系數上升，首先緩慢增加，而后增加速率越來越大，此趨勢對邊坡下部的節點更加明顯。采用了多種函數對2個邊坡的節點位移di-強度折減系數Fa關系進行了擬合，要求R2達到0.95以上，結果發現對于邊坡下部節點，采用4次多項式函數的擬合效果較好，而對于上部節點，采用3次多項式函數即可達到較好的擬合效果，表4顯示了對A邊坡5個節點的擬合情況。

圖3 A邊坡的節點位移-強度折減系數曲線

圖4 B邊坡的節點位移-強度折減系數曲線

表4 A邊坡5個節點的位移—折減系數擬合情況節點號擬合函數R2①d1=6.544 9 F3a-35.727 F2a+64.695 Fa-39.5210.965 5②d2=9.663 1 F3a-53.219 F2a+98.453 Fa-60.1450.977 8③d3=12.603 F3a-68.445 F2a+125.14 Fa-74.4340.954 7④d4=31.024 F4a-239.11 F3a+685.83 F2a-862.06 Fa+402.210.969 7⑤d5=44.494 F4a-342.71 F3a+981.64 F2a-1 233 Fa+576.230.988 7

3.4 利用模擬退火法對折減系數進行尋優

以圖2中各個位置開挖后的表面監測位移為依據，進行參數反分析，驗證利用模擬退火法對折減系數進行尋優的可靠性。其中監測數據來自對于觀測樁的定期觀測，觀測樁位于提取數據的節點處，開挖后進行位移觀測，直至位移趨于穩定，認為所發生的位移為開挖引起的位移。A邊坡和B邊坡5個節點的監測位移如表5所示。

表5 2個邊坡的位移監測值 mm邊坡①②③④⑤A0.3473.215.817.439.23B-0.7350.4361.783.614.97

建立如式(4)的目標函數C，為使C達到最小值，利用MATLAB軟件編寫模擬退火法程序進行Fa的全局化尋優，該程序的流程圖如圖5所示。

(4)

式中：dmi為節點i的位移監測值；di為節點i的位移計算值。

圖5 模擬退火程序流程圖

對于最優強度折減系數的搜索過程如圖6所示，最終得到2個邊坡的反分析最優強度折減系數分別為2.38和1.55；將初始的抗剪強度參數除以折減系數，則得到A邊坡的反分析抗剪強度參數為88.24 kPa(c)和42.02°(φ)，B邊坡的反分析抗剪強度參數為64.5 kPa(c)和29.03°(φ)。將反分析抗剪強度參數代入有限元模型中，再次計算5個節點的位移，并與監測值進行對比(見圖7和圖8)，可以看出，5個節點的數據點都比較靠近斜率為1的虛線段(監測位移=計算位移)，這表明利用反分析抗剪強度參數計算得出的位移可以較為準確地反映監測位移情況。

圖6 對于最優強度折減系數的搜索

圖7 A邊坡位移對比

圖8 B邊坡計算對比

4 結論

1)準確確定公路邊坡巖體的抗剪強度參數對邊坡設計具有重要意義。本文將模擬退火法(SA)引入到邊坡巖體抗剪強度參數反分析中，根據計算數據和監測數據設定目標函數，編寫模擬退火法尋優程序，搜索獲取最優的強度折減系數，再進一步得出反分析抗剪強度參數。

2)針對所述方法開展了2處邊坡的工程案例研究，將采用模擬退火法反分析得出的抗剪強度參數代入有限元模型中進行計算，各個節點的計算位移都比較接近監測位移，從而在一定程度上驗證了該方法的可靠性和有效性。