基于粒子群—有限差分耦合算法的隧道巖體力學(xué)參數(shù)反演①

孫 斌

(中鐵二十二局集團(tuán)有限公司，北京 100043)

由于巖土參數(shù)的不確定性和離散性，在數(shù)值正向計(jì)算過程中采用的初勘巖土力學(xué)參數(shù)勢(shì)必帶來一定的盲目性。反演分析作為一種基礎(chǔ)信息來自工程現(xiàn)場(chǎng)的“函數(shù)擬合”方法，其目的就是要依據(jù)獲得的量測(cè)信息來更好地認(rèn)識(shí)模型和描述模型的參數(shù)，以期獲得與系統(tǒng)輸出信息高度吻合的系統(tǒng)輸入信息，最終目標(biāo)仍然是更好地為工程設(shè)計(jì)、施工服務(wù)。

目前，用于工程領(lǐng)域的巖體力學(xué)及初始應(yīng)力場(chǎng)參數(shù)的反演方法較多基于有限元、有限差分和邊界元等數(shù)值方法的優(yōu)化反演。即利用正分析的過程和格式，通過優(yōu)化算法迭代最小誤差函數(shù)，逐次修正迭代過程中待反演參數(shù)的試算值，直至獲得“最佳值”。可見，問題實(shí)質(zhì)轉(zhuǎn)換為一個(gè)函數(shù)逼近問題，算法的優(yōu)劣性在一定程度上影響著“最佳值”。常用的優(yōu)化算法有:單純形法、復(fù)合型法、變量輪換法、共軛梯度法、罰函數(shù)法及鮑威爾法等。Gioda等［1－2］分別提出將單純形法、Rosenbrock 法、逆梯度法和鮑威爾法4種優(yōu)化算法應(yīng)用于求解巖體彈性及彈塑性力學(xué)參數(shù)，詳盡闡述了各方法在反分析中的適用性及不足;Arai［3］采用二次梯度法提出求解彈性模量和泊松比的新方法。然而，傳統(tǒng)的優(yōu)化方法由于計(jì)算模型的選擇而帶來待反演參數(shù)個(gè)數(shù)限制、尋優(yōu)過程過渡依賴初始值、計(jì)算效率低下等劣勢(shì)。為此，本文基于巖土數(shù)值仿真FLAC軟件的二次開發(fā)平臺(tái)，將FLAC2D平面有限差分與近些年來興起的一支仿生優(yōu)化算法—粒子群優(yōu)化算法相耦合，提出進(jìn)化的PSO－FLAC2D方法。該方法具有全局優(yōu)化和并行搜索的特點(diǎn)，能解決以往位移反演分析過程中由于計(jì)算模型的選擇而帶來的待反演參數(shù)個(gè)數(shù)的限制，同時(shí)在參數(shù)尋優(yōu)過程中，算法對(duì)初值取值不再依賴，融合算法計(jì)算效率高，收斂性好。最后，通過某隧道工程區(qū)巖體力學(xué)參數(shù)反演驗(yàn)證該方法的有效性。

1 粒子群優(yōu)化算法簡(jiǎn)介

粒子群算法(Particle Swarm optimization，簡(jiǎn)稱PSO)PSO算法就是從這種生物種群行為特性中得到啟發(fā)并用于求解優(yōu)化問題。在PSO中，每個(gè)優(yōu)化問題的潛在解都可以想象成d維搜索空間上的一個(gè)點(diǎn)，我們稱之為“粒子”(Particle)，所有的粒子都有一個(gè)被目標(biāo)函數(shù)決定的適應(yīng)度值(Fitness Value)，每個(gè)粒子還有一個(gè)速度決定他們飛翔的方向和距離，然后粒子們就追隨當(dāng)前的最優(yōu)粒子在解空間中搜索。Reynolds對(duì)鳥群飛行的研究發(fā)現(xiàn):鳥僅僅是追蹤它有限數(shù)量的鄰居但最終的整體結(jié)果是整個(gè)鳥群好像在一個(gè)中心的控制之下。即復(fù)雜的全局行為是由簡(jiǎn)單規(guī)則的相互作用引起的。

不失一般性，以最小化系統(tǒng)適應(yīng)度函數(shù)f為優(yōu)化目標(biāo)，標(biāo)準(zhǔn)粒子群算法數(shù)學(xué)原理可描述如下［4］。

粒子自身速度與位置更新:

所有粒子的全局極值選取:

式中:xj，vi，p besti分別表示第 i個(gè)粒子當(dāng)前位置、速度及個(gè)體歷史最優(yōu)位置;g best代表所有粒子經(jīng)歷的最好位置;w代表慣性權(quán)重，取較大值有利于算法在較大范圍內(nèi)進(jìn)行搜索，而取較小值有利于在局部進(jìn)行精細(xì)搜索;c1為個(gè)體進(jìn)化因子，c2為社會(huì)進(jìn)化因子，取值均在［0，2］;r1和r2是均服從U(0，1)分布的相互獨(dú)立隨機(jī)向量。

2 PSO－FLAC2D程序設(shè)計(jì)

2.1 待反演圍巖參數(shù)取值設(shè)計(jì)

基于PSO－FLAC2D方法的圍巖力學(xué)參數(shù)反演過程是將進(jìn)化算法與數(shù)值計(jì)算融合形成的進(jìn)化有限差分方法，進(jìn)化算法每迭代一次均需調(diào)用一次數(shù)值計(jì)算程序，而對(duì)數(shù)值計(jì)算結(jié)果產(chǎn)生影響的因素主要包括地質(zhì)因素和工程(支護(hù))因素。其中地質(zhì)因素又包括巖體物理力學(xué)參數(shù)和初始地應(yīng)力。具體而言，巖體的物理力學(xué)參數(shù)包括有變形模量E，黏聚力C，內(nèi)摩擦角φ，泊松比μ及剪脹角Φ;地應(yīng)力因素包括巖體容重γ，隧道埋深H及側(cè)壓力系數(shù)λ。

根據(jù)大量文獻(xiàn)研究表明:在圍巖力學(xué)參數(shù)反演過程中變形模量、黏聚力及內(nèi)摩擦角對(duì)數(shù)值計(jì)算結(jié)果影響較為敏感［5－8］，故在數(shù)值建模過程中，前 3個(gè)參數(shù)必須進(jìn)行反演，后2個(gè)參數(shù)不做反演，依據(jù)經(jīng)驗(yàn)取值。具體為:變形模量依照工程區(qū)域土工試驗(yàn)的實(shí)測(cè)數(shù)據(jù)和《鐵路隧道設(shè)計(jì)規(guī)范TBJ10003－2005》的圍巖分級(jí)相關(guān)規(guī)定表綜合確定;黏聚力、內(nèi)摩擦角取值范圍參照《巖石力學(xué)參數(shù)手冊(cè)》中各類圍巖的統(tǒng)計(jì)參數(shù)規(guī)律表確定［9］;泊松比參考FLAC學(xué)習(xí)手冊(cè)中Goodman在1980年進(jìn)行的巖石彈性參數(shù)實(shí)驗(yàn)值來確定泊松比;剪脹角根據(jù)Vermeer和 de Borst(1984)［10］的研究成果，認(rèn)為無論材料是土體、巖石或是混凝土，其剪脹角的近似值分布在0°～20°區(qū)間。

2.2 PSO －FLAC2D程序?qū)崿F(xiàn)

FLAC內(nèi)嵌了FISH語言，使用戶能自定義變量和函數(shù)，擴(kuò)大了FLAC的應(yīng)用及用戶自有的特色。其基本功能包括文件讀取操作，自定義變量與函數(shù)，流程控制語句等?；贔ISH語言編制如下反分析程序，現(xiàn)就主程序 PSOTunnel.dat說明如下。

%＊＊＊＊＊＊＊＊＊＊＊PSOTunnel.dat＊＊＊＊＊＊＊＊＊＊%

new

sys cd E:ackanalysisPSO－Flac

restore model.sav

sys cd E:ackanalysisPSO－Flac

%＊＊＊＊＊＊＊＊＊＊＊＊＊初始化工作＊＊＊＊＊＊＊＊＊＊＊%

call IO.fis;

//調(diào)用IOfis文件，它是定義文件操作函數(shù)的文件call PSO.fis;

//調(diào)用PSO.fis，它是定義PSO優(yōu)化算法各函數(shù)的文件

DefineVar;

//該函數(shù)用于初始化各個(gè)變量，如體積模量、剪切模量等，見IO.fis

LoadPSOVar;

//該函數(shù)用于讀取PSO文件并初始化PSO變量，如迭代步等，見IO.fis

define GetParticle;

//定義GetParticle函數(shù)

if P_isInit=0 then;PInit;

else

LoadParticle;Preset;endif

end;

GetParticle ;

//調(diào)用GetParticle函數(shù)

InitProp

//該函數(shù)用于計(jì)算模型參數(shù)，如楊氏模量、泊松比等，見PSO.fis

call balance－kaiwa.txt;

//FLAC計(jì)算

%＊＊＊＊＊＊＊＊＊初始化工作完成＊＊＊＊＊＊＊＊＊%

PUpdateBest;

//該函數(shù)根據(jù)PSO算法原理判定粒子數(shù)據(jù)更新，見 PSO.fis

WriteParticle;

//該函數(shù)用于將當(dāng)前粒子信息寫入粒子文件中，見 IO.fis

WriteLogFile;

//該函數(shù)運(yùn)行日志將計(jì)算結(jié)果寫入日志，見IO.fis

PContinuable;

//該函數(shù)判定是否達(dá)到準(zhǔn)則函數(shù)要求，并指示是否繼續(xù)計(jì)算，見PSO.fis

WritePSOVar;

//該函數(shù)將當(dāng)前的PSO變量信息寫入PSO文件中，見 IO.fis

call PSOTunnel.dat;

//調(diào)用自身，使得開挖模擬計(jì)算循環(huán)進(jìn)行

%＊＊＊＊＊＊＊＊PSOTunnel.dat結(jié)束＊＊＊＊＊＊＊＊＊%

PSOTunnel.dat文件與 PSO.fis和 IO.fis 3 個(gè)文件共同實(shí)現(xiàn)了PSO－FLAC2D位移反分析算法。其中，PSO.fis和IO.fis文件中主要是關(guān)于9個(gè)子函數(shù)的定義和具體實(shí)現(xiàn)，在此，僅給出如圖1所示的程序結(jié)構(gòu)樹。

圖1 程序結(jié)構(gòu)Fig.1 Structure of PSO －FLAC2D programme

3 反演算例研究

以某隧道工程建設(shè)為背景，參照設(shè)計(jì)方前期勘察對(duì)場(chǎng)區(qū)附近巖土體試驗(yàn)成果的統(tǒng)計(jì)分析，采用《地下鐵道、輕軌交通巖土工程勘察規(guī)范》(GB50307－1999)及《鐵路隧道設(shè)計(jì)規(guī)范》(TBJ10003－2005)等推薦的力學(xué)參數(shù)換算方法，提出本場(chǎng)區(qū)待反演參數(shù)取值區(qū)間，如表1。

表1 待反演參數(shù)取值區(qū)間Table 1 Variation ranges of inversion parameters

由表1可知，剪脹角在本文反演中大致取平均值10°;由于該段巖性較為單一，其容重不做反演，按照實(shí)測(cè)參數(shù)取為27 kN/m3;泊松比參考FLAC學(xué)習(xí)手冊(cè)中Goodman在1980年進(jìn)行的巖石彈性參數(shù)實(shí)驗(yàn)值來確定泊松比為0.30。至此，在明確了圍巖參數(shù)的取值區(qū)間之后，確定正演反分析處理中涉及的基本要素:(1)合理的數(shù)值模擬計(jì)算模型—FLAC2D;(2)優(yōu)良的參數(shù)反分析算法—PSO;(3)有效的準(zhǔn)則函數(shù)—誤差平方和最小化函數(shù)，見式(4)。

式中，u'i和ui分別表示i點(diǎn)的數(shù)值計(jì)算位移和現(xiàn)場(chǎng)監(jiān)測(cè)位移;其中，現(xiàn)場(chǎng)監(jiān)測(cè)位移采用由中鐵西南科學(xué)研究院提供的斷面ZK2+906拱頂沉降監(jiān)測(cè)資料和測(cè)點(diǎn)布置方案，見圖2。

圖2 斷面ZK2+906測(cè)點(diǎn)布置及其沉降曲線Fig.2 Layout of measuring points and crown settlement curve at ZK2+906

粒子群算法(優(yōu)化算法)參數(shù)設(shè)置如下:待優(yōu)化參數(shù)個(gè)數(shù)為Nparams=3，群體規(guī)模Popsize=10，進(jìn)化代數(shù) Maxiter=10，權(quán)重區(qū)間為［0.3，0.9］，c1=c2=2.0?？偤臅r(shí)約3 h左右(平面反分析耗時(shí))。為了進(jìn)一步體現(xiàn)本文提出的平面反演技術(shù)的效率和精度，重復(fù)如上反演操作，對(duì)ZK2+888，ZK2+865及ZK2+836斷面分別進(jìn)行平面反演。粒子群算法反演得到的巖體力學(xué)參數(shù)見表2。

表2 反演得到的巖體力學(xué)參數(shù)Table 2 The feedback rock mass mechanics

依據(jù)表2的反演結(jié)果和《鐵路隧道設(shè)計(jì)規(guī)范》關(guān)于巖體物理力學(xué)參數(shù)的規(guī)定可知，反演結(jié)果為四級(jí)圍巖，這一初步結(jié)論與現(xiàn)場(chǎng)圍巖變更結(jié)果一致。其次，不同斷面的反演結(jié)果不盡相同。這是因?yàn)閺牧W(xué)概念上講，反演的圍巖力學(xué)參數(shù)僅僅可看作一種“等效值”或“近似解”，其并無明確物理意義。同時(shí)，對(duì)于四級(jí)及四級(jí)以上的圍巖而言，在數(shù)值計(jì)算過程中圍巖的黏聚力C值對(duì)計(jì)算結(jié)果十分敏感，表中C值的大小與各個(gè)斷面的變形實(shí)測(cè)值大小規(guī)律匹配，即ZK2+836的拱頂變形數(shù)據(jù)最大，達(dá)到14.5 mm，而其對(duì)應(yīng)的C值也為最小?？傮w而言，平面反演相比三維反演數(shù)據(jù)處理工作量小，反演結(jié)果偏于保守。實(shí)際工程應(yīng)用中可取多個(gè)斷面的反演平均值確定局部區(qū)域的巖體力學(xué)參數(shù)，為工程設(shè)計(jì)、變更決策提供理論借鑒。

4 結(jié)論

(1)將粒子群算法與FLAC2D融合形成進(jìn)化有限差分方法(PSO－FLAC2D)，繼而提出了圍巖物理力學(xué)參數(shù)的位移智能反演方法及實(shí)現(xiàn)步驟。

(2)從程序構(gòu)架上看，該算法不受反演參數(shù)數(shù)目限制且優(yōu)化算法具有系統(tǒng)參數(shù)少、并行搜索、全局最優(yōu)的優(yōu)勢(shì)。

(3)從工程算例結(jié)果看，該方法反演結(jié)果科學(xué)合理且具有較高的計(jì)算效率，平面反演的效率可控制在數(shù)小時(shí)以內(nèi)，有效地降低數(shù)值試驗(yàn)預(yù)設(shè)方案數(shù)量，為巖體力學(xué)參數(shù)反演提出一種新的途徑。

［1］Gioda G，Maier G.Direct search solution of an inverse

problem in elasto－plasticity identification of cohesion friction angle and in－situ stress by pressure tunnel tests［J］.Int J Num Methods in Eng，1980，15:1823 －1834.

［2］Gioda G，Pandolfi A，Cividini A.A comparative evaluation of some back analysis algorithms and their application to in－ situ load tests［C］//Proc 2nd Int Symp on Field Measurement in Geom.Beijing:Science Press，1987:1131－1144.

［3］Arai R.An inverse problem approach to the prediction of Multi－ dimensional consolidation behavior［J］.Soil and Foundations，1984，24(1):95 －108.

［4］曾建潮，介婧，崔志華.微粒群算法［M］.北京:科學(xué)出版社，2004.ZENG Jianchao，JIE Jing，CUI Zhihua.Particle swarm optimization algorithm［M］.Beijing:Science Press，2004.

［5］馮夏庭，張志強(qiáng)，楊成祥，等.位移反分析的進(jìn)化神經(jīng)元網(wǎng)絡(luò)方法研究［J］.巖石力學(xué)與工程學(xué)報(bào)，1999，18(5):529－533.FENG Xiating，ZHANG Zhiqiang，YANG Chengxiang.Study on genetic－neural network Method of displacement back analysis［J］.Chinese Journal of Rock Mechanics and Engineering，1999，18(5):529 －533.

［6］戴榮，李忠奎.三維地應(yīng)力場(chǎng)BP反分析的改進(jìn)［J］.巖石力學(xué)與工程學(xué)報(bào)，2005，24(1):83－88.DAI Rong，LI Zhongkui.Modified BPback analysis of 3d in － situ stresses［J］.Chinese Journal of Rock Mechanics and Engineering，2005，24(1):83 －88.

［7］Hisatake M.Three dimensional back analysis for tunnels［C］//Proc Int Symp on ECRF.Beijing:Science Press，1986:791－797.

［8］ WANG Zhiyin，LI Yunpeng.Back analysis of measured rheologic displacements of underground opening［C］//Proc Int Symp on Underground Eng，New Delhi，1988:181－186.

［9］葉金漢.巖石力學(xué)手冊(cè)［M］.北京:水利電力出版社，1991.YE Jinhan.Rock mechanics Handbook［M］.Beijing:China Water Power Press，1991.

［10］Vermeer PA，Borst R de.Non－associated plasticity for soils［J］.Concrete and Rock，1984，29(3):1 －64.