基于单纯形-模拟退火算法的岩土力学参数反演

3

(1.长江大学 城市建设学院，湖北 荆州 434023； 2.中国一冶集团有限公司，湖北 武汉 430081； 3.东北石油大学 非常规油气研究院，黑龙江 大庆 163318)

近年来，随着计算机技术的飞速发展和各种数值计算方法的广泛应用，运用数值方法研究岩土工程实际问题已成为常态。为了使数值计算结果更加准确、客观、实用，计算参数的选取尤为重要，但因为岩土介质的复杂性，在实际中如何选取计算参数值就成为了一个难题。目前确定岩土参数的途径主要有试验法和反演分析法。由于室内试验或现场原位试验确定的相应岩土力学参数都与实际结果有着较大的偏差，因此得到的数值计算结果存在一定的误差，这在一定程度上也限制了数值方法在实际工程中的应用。而基于实测信息的反演分析方法则为岩土力学参数的获取提供了新的途径，且表现出了独特的优势[1-2]。由于岩体结构在开挖或变形过程中位移是最容易测定的，利用现场位移量测信息反演岩土力学参数的位移反分析法在岩土工程中的应用较为普遍。

位移反分析法是将参数反演问题转化为目标函数寻优问题，但是目标函数是一个很复杂的非线性多峰函数，传统优化方法存在着优化结果依赖于初值选择、容易陷入局部极值等不足，采用全局优化算法成为较理想的选择。一些学者将遗传算法、粒子群算法等智能优化算法应用到岩土参数的反演计算中，取得了较为理想的结果，但是存在着计算量大、收敛速度慢以及计算时间长等不足，而且这些方法也并没有真正将有限元程序嵌入到相应优化算法中[3-7]。一些学者将自己开发的有限元程序嵌入到优化算法中，编写了优化反分析程序，但自编的有限元程序其功能不如大型商业有限元软件优越，所以也在一定程度上限制了反演分析的功能，对复杂本构模型以及多场耦合参数的反演往往无法实现[2，8-9]。还有一些学者将节理岩体视为横观各向同性介质对岩土参数进行反演，忽略了节理面的影响，反演结果与实际存在着较大误差[10-11]。

本文将单纯形算法和模拟退火算法相结合组成混合算法，并基于该算法编写优化反分析程序，将ABAQUS作为一个单独模块嵌入到优化算法程序中，实现MATLAB和ABAQUS的联合反演，加快优化算法的收敛速度，减少反演计算的迭代次数。基于ABAQUS自带的节理材料模型，编写了节理岩体各向异性弹塑性本构模型子程序USDFLD，并应用到实际工程数值模型中，提高了反演参数精度。

1 反演模型

位移反分析是利用现场实测的变形作为观测信息反演岩土力学参数。待反演参数可用向量x表示，即x=[x1,x2,…,xm]T，其中m表示待反演参数总个数。以计算值与测点实测值之间的误差建立目标函数如下：

(1)

在参数优化反演过程中，寻找一组材料参数{x*}，构造目标函数的均值为新的目标函数，使新目标函数得到最小值时的反演参数即视为现场岩体材料参数，此时目标函数值为

(2)

岩土工程反演问题通常较为复杂，为保证解的稳定性和唯一性需要加一些约束条件作为约束优化问题来处理，则有限元优化反演模型可表示为

(3)

式中，Mi(x)为第i个等式约束函数；Nj(x)为第j个不等式约束。

采用精确罚函数法将上述约束优化问题变成无约束优化问题，即：

(4)

2 反演步骤与程序实现

由于岩土工程反演模型较为复杂，无法准确写出其解析表达式，用计算目标函数的导数来计算出目标函数极小值就无法实现，因此本文采用单纯形-模拟退火优化算法来解决这个问题。

2.1 单纯形-模拟退火优化算法

单纯形优化算法是先求出一个基本可行解，然后逐步改善可行解，使目标函数值逐步增大(或减小)，直到目标函数达到极值时，就可得到最优解。模拟退火优化算法是基于Monte-Carlo迭代求解策略的一种随机寻优算法，是源于组合优化与物理中固体物质退火过程之间的相似性。模拟退火优化算法是在某一较高初温下，伴随温度参数的不断下降，结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解中能概率性地跳出并最终趋于全局最优解。

将两种算法结合形成一种混合算法，即单纯形-模拟退火优化算法，该算法的主要计算步骤如下。

(1) 给出待反演参数的初值x0，初解为G(x0,u0)，设定算法控制求解精度。

(4) 利用Metropolis准则判别函数收敛性。

(5) 若优化后目标函数值大于上步的迭代函数值，则将上步的反演参数替代该反演参数，返回步骤(2)。

(6) 设置线性温度降低函数，初温设为某一固定值。

2.2 程序研制

基于本文提出的反演思路和方法，在ABAQUS的基础上，采用MATLAB语言为平台并结合单纯形-模拟退火算法编制了优化反演子程序SMSAInverse.m。优化反分析程序的具体实施步骤如下。

(1) 用MATLAB语言编写反演主程序，给出待反演参数的初值和测点的实测位移值，并用MATLAB语言调用ABAQUS的命令流文件进行计算。

(2) 用MATLAB语言编写目标函数子程序，调用ABAQUS的计算结果dat文件，获取测点的位移计算值，以此计算目标函数。

(3) 编写单纯形-模拟退火算法，在算法中定义相应的边界条件并初始化其参数。

(4) 先用单纯形算法优化出待反演参数，再用strrep命令修正命令流文件中的待反演参数，当目标函数值满足精度要求时，以反演后的参数作为模拟退火算法的初值进行优化迭代，并同样通过修改命令流文件更新反演参数，直到满足迭代次数或计算精度要求为止，停止迭代。

优化反分析程序如图1所示。

3 材料本构模型

目前理论上多将节理岩体视为横观各向同性介质[12]，其应力-应变关系包含5个弹性参数，但由于岩体中存在节理面，节理岩体和岩石的力学特征存在较大的差别，因此节理岩体表现出各向异性特征。由横观各向同性理论[13]可知，局部坐标系中横观各向同性增量的弹性应力-应变关系表达式如下：

图1 岩土力学参数优化反演程序框图Fig.1 Flow chart of inversion analysis for mechanical parameters of geotechnical materials

[Δσ′]=[K′][Δε′e]

(5)

式中,[K′]为局部对称刚度矩阵，写入对应的应力应变对号，则可表示为

(6)

岩石应力应变的关系不仅与这5个弹性参数有关，还与节理面倾角有关[14]，由各向异性理论可知，轴向应力σy与应变(εx,εy,εz)的关系如下：

εx=a12σy,εy=a22σy,εz=a23σy

(7)

式中,a12,a22,a23是由5个弹性参数和节理面倾角θ决定的函数关系式，表达式如下：

(8)

(9)

(10)

当5个弹性参数确定时，即可确定任一节理面倾角的弹性参数，即：

(11)

基于上述理论分析，本文以ABAQUS自带的节理材料模型为基础，结合节理岩体弹性参数的各向异性，建立了节理岩体各向异性弹塑性本构关系，并用FORTRAN语言开发了材料本构模型子程序USDFLD，将其应用到本文数值模拟的模型当中。

4 工程应用

4.1 工程概况

武汉花山大道宝盖山隧道工程采用小净距分离式隧道，隧道桩号为K0+960～K1+542，隧道全长582 m，穿越宝盖山与钵孟峰山谷，隧道线型为直线，具体概况详见武汉花山大道工程报告[15]。本文选取宝盖山隧道左线桩号K1+360为研究对象，隧道围岩属于Ⅲ级围岩，采用上下台阶开挖法，该隧道支护方式如图2所示，采用复合式衬砌设计，并钉入注浆锚杆，混凝土强度等级为C25，系统锚杆采用Φ22/Φ25的中空注浆锚杆，二次模筑衬砌厚度为450 mm。隧道开挖的主要施工步骤如图3所示。

图2 围岩支护结构布置Fig.2 Rock support structure layout

图3 隧道开挖步骤Fig.3 Excavation steps

4.2 数值模型

由宝盖山隧道左线K1+360处地质资料可知，隧洞方向长度约为243 m，隧道高度方向为10.157 m，跨度方向长度为13.75 m。以此建立隧洞断面二维模型，模型的竖直方向高度为99.5 m，水平方向长度为160 m。衬砌厚度为450 mm，单元类型为CPE4，且锚杆采用半径为25mm的B21杆单元，衬砌和锚杆均采用节点公用单元。模型两侧施加水平方向位移约束，底边施加竖直方向固定约束，单元类型为CPE4(四节点平面应变单元)，有限元数值模型及测点分布如图4所示。

图4 有限元网格模型及测点布置Fig.4 Measuring points and finite element mesh model

如图4所示，参与反演的测点为A、B、C、D、E五个测点，其中A、B、C、D为水平方向测点，E为竖直方向测点，测点对应的有限元模型中单元节点分别为991、699、28、30、48，有部分测点与单元节点无法对应，采用插入法取测点位移值，则该5个测点实测位移值如表2所示。

表2 测点位移实测值Tab.2 Measured displacements of selected points mm

笔者用正交设计试验方法分析了岩土力学参数对于围岩变形的敏感性，最后选取了节理岩体凝聚力、内摩擦角和节理面凝聚力作为待反演参数，具体过程在此不再赘述。将USDFLD材料本构子程序应用到隧道数值模型中，且模型中参与反演的材料参数的初值及取值范围如表3所示。

表3 待反演参数设置Tab.3 Initial setting of undetermined parameters for inversion model

隧洞围岩参数中弹性模量为0.244 GPa，泊松比取0.18，密度为2 500 kg/m3，节理岩体凝聚力为3.2 MPa，节理岩体内摩擦角为24.001°，节理面凝聚力为0.3 MPa，节理面内摩擦角为10°；衬砌的弹性模量为9 GPa，泊松比为0.25，密度为2 500 kg/m3；锚杆弹性模量为180 GPa，泊松比为0.22，密度为7 800 kg/m3。在数值模拟过程中，复制开挖土体并改变其材料参数，在开挖原始土体后添加弹性模量较小的土体，以此达到弹性模量衰减的目的。

数值模拟开挖的过程如下：创建二维数值模型，设定材料参数和荷载边界条件，进行地应力平衡；复制①、②两部分土体，且每部分复制两份，首先开挖原始土体①部分并添加模量较小的土体以达到模量衰减的目的，然后施加上部衬砌；同样的方法开挖②部分土体，然后施加下部衬砌，开挖结束。

4.3 反演结果

对上述二维数值模型进行有限元计算，将实际测点对应的模型测点的位移值作为理论计算值，依据该计算值可优化反演出相关待反演参数，并通过比较实测点的计算位移值与实测位移值对本文方法的合理性和精度做出判断。

利用单纯形-模拟退火算法进行优化迭代，首先以待反演参数的反演区间上下限的平均值作为初始值，由单纯形算法进行优化迭代，得出反演参数值。然后以单纯形反演值作为初始值，进行模拟退火算法优化计算，最终得到待反演参数的反演结果，如表4所示。最终优化程序迭代了65次，调用有限元软件ABAQUS计算177次。反演求解迭代次数曲线如图5所示。

表4 参数反演结果Tab.4 Back analysis results for mechanical parameters

图5 目标函数值变化趋势Fig.5 Iteration process curve of back analysis

测点位移的反演结果和与实测值比较的结果如表5所示。可知，实测值与反演值误差较小，30测点的Ux位移值误差较大，但整体位移值较小，其它测点的位移相对误差均控制在5%以内。

表5 反演结果与实测值的比较Tab.5 Comparison of measured results and numerical results for measuring point displacement

由上可知，由单纯形局部优化算法与模拟退火全局优化算法相结合的混合算法，解决了单纯形算法易陷入局部极小值和模拟退火算法搜索效率低的问题，提高了搜索效率和反演参数精度，为确定地下工程围岩力学参数提供了有效的途径。将自行开发的节理本构模型应用到花山大道宝盖山隧道工程当中，能够很好地反映隧洞围岩变形及受力特性，为指导实际工程起到了一定的预见作用。

5 结 论

(1) 将单纯形优化算法和模拟退火算法相结合，构造一种混合算法并结合合适的精确罚函数，建立了基于计算值和测点实测值之间的目标函数，并将有限元软件ABAQUS作为单独模块嵌入到单纯形-模拟退火算法中，实现ABAQUS与MATLAB联合反演，可保证该算法不依赖于初始值就能搜索到全局最优解，且可加快搜索速度。

(2) 建立花山大道宝盖山隧道开挖的二维数值模型，将开发的节理岩体各向异性本构子程序应用到数值模型中，并嵌入到编写的优化反演程序中，对节理岩体凝聚力、节理岩体内摩擦角和节理面凝聚力三个力学参数进行反演，得到了合理的参数值，说明本文提出的方法可用于实际工程，为确定围岩力学参数提供了一套有效的方法。