基于粒子群—有限差分耦合算法的隧道岩体力学参数反演①

孙 斌

(中铁二十二局集团有限公司，北京 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的应用及用户自有的特色。其基本功能包括文件读取操作，自定义变量与函数，流程控制语句等。基于FISH语言编制如下反分析程序，现就主程序 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。总耗时约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.