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

孫 斌

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

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

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

1 粒子群優化算法簡介

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

不失一般性，以最小化系統適應度函數f為優化目標，標準粒子群算法數學原理可描述如下［4］。

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

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

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

2 PSO－FLAC2D程序設計

2.1 待反演圍巖參數取值設計

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

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

2.2 PSO －FLAC2D程序實現

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

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

new

sys cd E:ackanalysisPSO－Flac

restore model.sav

sys cd E:ackanalysisPSO－Flac

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

call IO.fis;

//調用IOfis文件，它是定義文件操作函數的文件call PSO.fis;

//調用PSO.fis，它是定義PSO優化算法各函數的文件

DefineVar;

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

LoadPSOVar;

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

define GetParticle;

//定義GetParticle函數

if P_isInit=0 then;PInit;

else

LoadParticle;Preset;endif

end;

GetParticle ;

//調用GetParticle函數

InitProp

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

call balance－kaiwa.txt;

//FLAC計算

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

PUpdateBest;

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

WriteParticle;

//該函數用于將當前粒子信息寫入粒子文件中，見 IO.fis

WriteLogFile;

//該函數運行日志將計算結果寫入日志，見IO.fis

PContinuable;

//該函數判定是否達到準則函數要求，并指示是否繼續計算，見PSO.fis

WritePSOVar;

//該函數將當前的PSO變量信息寫入PSO文件中，見 IO.fis

call PSOTunnel.dat;

//調用自身，使得開挖模擬計算循環進行

%＊＊＊＊＊＊＊＊PSOTunnel.dat結束＊＊＊＊＊＊＊＊＊%

PSOTunnel.dat文件與 PSO.fis和 IO.fis 3 個文件共同實現了PSO－FLAC2D位移反分析算法。其中，PSO.fis和IO.fis文件中主要是關于9個子函數的定義和具體實現，在此，僅給出如圖1所示的程序結構樹。

圖1 程序結構Fig.1 Structure of PSO －FLAC2D programme

3 反演算例研究

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

表1 待反演參數取值區間Table 1 Variation ranges of inversion parameters

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

式中，u'i和ui分別表示i點的數值計算位移和現場監測位移;其中，現場監測位移采用由中鐵西南科學研究院提供的斷面ZK2+906拱頂沉降監測資料和測點布置方案，見圖2。

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

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

表2 反演得到的巖體力學參數Table 2 The feedback rock mass mechanics

依據表2的反演結果和《鐵路隧道設計規范》關于巖體物理力學參數的規定可知，反演結果為四級圍巖，這一初步結論與現場圍巖變更結果一致。其次，不同斷面的反演結果不盡相同。這是因為從力學概念上講，反演的圍巖力學參數僅僅可看作一種“等效值”或“近似解”，其并無明確物理意義。同時，對于四級及四級以上的圍巖而言，在數值計算過程中圍巖的黏聚力C值對計算結果十分敏感，表中C值的大小與各個斷面的變形實測值大小規律匹配，即ZK2+836的拱頂變形數據最大，達到14.5 mm，而其對應的C值也為最小。總體而言，平面反演相比三維反演數據處理工作量小，反演結果偏于保守。實際工程應用中可取多個斷面的反演平均值確定局部區域的巖體力學參數，為工程設計、變更決策提供理論借鑒。

4 結論

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

(2)從程序構架上看，該算法不受反演參數數目限制且優化算法具有系統參數少、并行搜索、全局最優的優勢。

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

［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］.北京:科學出版社，2004.ZENG Jianchao，JIE Jing，CUI Zhihua.Particle swarm optimization algorithm［M］.Beijing:Science Press，2004.

［5］馮夏庭，張志強，楊成祥，等.位移反分析的進化神經元網絡方法研究［J］.巖石力學與工程學報，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］戴榮，李忠奎.三維地應力場BP反分析的改進［J］.巖石力學與工程學報，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］葉金漢.巖石力學手冊［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.