張曉軍, 曹 娟

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

0 引言

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

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

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

1 模擬退火原理

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

根據(jù)統(tǒng)計熱力學(xué)定律，分子在溫度T時，處于某種狀態(tài)Ei滿足Boltzmann率分布：

(1)

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

物體從狀態(tài)i躍遷到狀態(tài)j的概率為：

(2)

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

P(mi→mj)=

(3)

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

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

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

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

2)建立各個節(jié)點的位移-折減系數(shù)關(guān)系曲線，并進行函數(shù)擬合。

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

4)根據(jù)計算數(shù)據(jù)和監(jiān)測數(shù)據(jù)設(shè)定目標(biāo)函數(shù)，編寫模擬退火法尋優(yōu)程序，搜索獲取最優(yōu)強度折減系數(shù)。

5)將初始抗剪強度參數(shù)按最優(yōu)強度折減系數(shù)進行折減，得到反分析結(jié)果，進行誤差分析。

3 工程案例

3.1 工程概況

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

a) A邊坡

b) B邊坡

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

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

3.2 有限元模型分析

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

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

a) A邊坡

b) B邊坡

本文中強度折減法僅對模型上部的砂質(zhì)板巖進行強度折減，強度折減系數(shù)與抗剪強度參數(shù)對應(yīng)關(guān)系如表2、表3所示。

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

表2 A邊坡在不同折減系數(shù)下的抗剪強度參數(shù)折減系數(shù)粘聚力/kPa內(nèi)摩擦角/(°)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邊坡在不同折減系數(shù)下的抗剪強度參數(shù)折減系數(shù)粘聚力/kPa內(nèi)摩擦角/(°)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 位移-折減系數(shù)曲線

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

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

圖3 A邊坡的節(jié)點位移-強度折減系數(shù)曲線

圖4 B邊坡的節(jié)點位移-強度折減系數(shù)曲線

表4 A邊坡5個節(jié)點的位移—折減系數(shù)擬合情況節(jié)點號擬合函數(shù)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 利用模擬退火法對折減系數(shù)進行尋優(yōu)

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

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

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

(4)

式中：dmi為節(jié)點i的位移監(jiān)測值；di為節(jié)點i的位移計算值。

圖5 模擬退火程序流程圖

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

圖6 對于最優(yōu)強度折減系數(shù)的搜索

圖7 A邊坡位移對比

圖8 B邊坡計算對比

4 結(jié)論

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

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