張曉軍, 曹 娟

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

0 引言

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

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

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

1 模擬退火原理

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

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

(1)

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

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

(2)

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

Pmi→mj=

(3)

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

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

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

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

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

3)實(shí)際邊坡開(kāi)挖后，在上述各個(gè)節(jié)點(diǎn)進(jìn)行位移監(jiān)測(cè)，獲取邊坡開(kāi)挖后上述各個(gè)節(jié)點(diǎn)的監(jiān)測(cè)位移穩(wěn)定值。

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

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

3 工程案例

3.1 工程概況

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

a) A邊坡

b) B邊坡

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

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

3.2 有限元模型分析

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

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

a) A邊坡

b) B邊坡

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

表1 初始材料參數(shù)材料重度/kN·m-3 泊松比彈性模量/MPa粘聚力/kPa內(nèi)摩擦角/(°)強(qiáng)風(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ù)下的抗剪強(qiáng)度參數(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ù)下的抗剪強(qiáng)度參數(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ù)，對(duì)2個(gè)邊坡分別進(jìn)行8組折減系數(shù)下的有限元分析，提取5個(gè)節(jié)點(diǎn)(如圖2所示)監(jiān)測(cè)開(kāi)挖后的位移，得到節(jié)點(diǎn)位移-強(qiáng)度折減系數(shù)曲線如圖3和圖4所示。

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

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

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

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

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

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

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

圖5 模擬退火程序流程圖

(4)

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

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

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

圖7 A邊坡位移對(duì)比

圖8 B邊坡計(jì)算對(duì)比

4 結(jié)論

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

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