基于DE-FEM 隧道反分析可视化平台研究及应用

姜谙男，李 鹏，唐述林，王军祥

（1. 大连海事大学 道路与桥梁研究所，辽宁 大连 116026；2. 中铁21 局集团有限公司，兰州 730000）

1 引 言

由于地质体环境的复杂性和不确定性，隧道预设计施工方案不可避免地存在盲目性。利用施工观测信息，对隧道围岩性质和状态进行反馈分析，是工程界提倡的新奥法施工、信息化施工的精髓。隧道反馈分析包括定性经验分析和定量计算分析，后者主要是参数反分析，即从力学模型层次上进行定量识别分析，具有重要的学术价值和工程意义[1]。

岩土反分析的基本思想最先由Kavangh、Gioda和Maier 等提出[1]。日本学者樱井春辅进行了岩体弹性模量及初始地应力的线弹性有限元位移反分析研究[2]。此后国内外许多学者在此方面进行了大量的研究工作，发展了弹塑性、黏弹性、黏塑性等非线性反分析，以及有关量测误差处理、优化、校验等技术[3－5]。冯夏庭[6]将智能科学引入岩石力学与工程反分析取得良好效果。近年来国内也出现了用于隧道工程的反分析软件，例如同济曙光、BMP 反分析软件等[7]。

尽管反分析已经取得较多学术理论成果，但也应看到，实际工程中反分析的应用仍然不尽如人 意[8]。制约反分析应用的主要原因是：（1）岩土工程的非线性特点使传统优化方法容易限于局部最小，而且线弹性简化条件与实际相差太大；（2）传统反分析数据抽象，对专业知识和理论有较高要求，现场人员不易使用。因此，寻求鲁棒性强的优化算法并开发界面友好的反分析可视化平台对于促进反分析工程应用，提高隧道安全保障具有重要的意义。

2 反分析原理及DE-FEM 步骤

2.1 反分析原理

反分析是指利用现场观测物理量反算出围岩初始地应力、力学参数等。反分析分为解析法和数值法，数值法又分正反分析、逆反分析和图谱法3种。其中正反分析法可以利用现有正演计算程序，具有很强的适应性。其原理步骤如下：选取一组参数计算得到围岩变形，如果计算变形与监测变形数据相差较大，则调整输入参数重新计算，直到计算与监测数据相差很小，这时对应参数即为识别的参数。由此可见，参数反分析本质上是优化问题。根据隧道围岩参数的特定物理意义设定上下限，如果区域内有m 个观测值，那么有约束的优化问题为

2.2 基于DE-FEM 的反分析过程

2.2.1 差异进化的优化步骤

目前以遗传算法、粒子群算法和模拟退火算法等为代表的进化优化算法，以其良好的鲁棒性、全局寻优能力和不需要对目标函数求导数等优点显示出强大生命力。其中差异进化算法（Difference Evolutionary，DE）由美国加州大学伯克利大学的Kenneth Price 和Rainer Storm 提出的。该算法虽源于遗传算法，由于不进行编码和解码操作，使用上大为简化。同时它对初始值无要求，收敛速度快，对各种非线性函数适应性强[9]。

式（1）的Yi是以围岩参数xi（x1, x2,…, xn）为自变量的复杂非线性隐式函数，引进差异进化算法实现对参数xi全局优化。假设需要识别n 个岩石力学参数，则首先产生维数为n 的Np个向量，Np称为种群规模，每个向量也就是一组潜在的解，称为个体。对每一个个体向量按式（1）进行计算，以计算结果作为评价值。根据评价值大小按照DE 规则迭代，逐渐接近最优解。将有限元算法嵌入到差异进化算法中，形成DE-FEM 耦合的反分析计算流程见图1。如图所示，DE-FEM 的计算过程如下：

（1）产生初始种群。选择需要反分析的参数（自变量）个数为n，首先随机产生满足自变量上、下界约束的NP个n 维向量，公式如下：

图1 岩土反分析DE-FEM 计算流程 Fig.1 The DE-FEM calculation process of back analysis

（2）变异操作。缩放种群中任意2 个向量个体之间的差值并叠加到种群中的第3 个向量个体上，形成新的变异变量。第G +1 代第i 个变异向量的第j 个分量为

式中：下标r1，r2，r3为1～NP中的随机整数且互不相同；F 为缩放因子，用来调节向量差异的步长幅值，在0～2 内取值。式（3）是差异进化算法的基本变异模式，被称作DE/rand/1 模式。差异进化还有其他模式，如DE/best/1、DE/best/2、DE/rand/2等[9]。

（3）交叉操作。将目标向量xi(G)与变异向量vi(G +1)按照如下规则杂交，生成新的试样向量ui(G +1)，其第j 个分量表示为

式中：rj∈[0，1]，为与向量第j 个分量对应的随机数；CR∈[0，1]，为杂交概率常数；rni为在1，2，…，n 中随机挑选一个整数，以确保变异向量Vi(G +1)中，至少有一个分量被试样向量ui(G +1)采用。

（4）选择。采用贪婪搜索方法进行选择操作。将试样向量ui(G +1)和目标向量xi(G)比较，调用有限元程序计算，采用式（1）的目标函数作为评价函数，如果ui(G+1)对应较小的目标函数值，则选择ui(G +1)；反之，如果xi(G)对应较小的目标函数值，则保留xi(G)。

（5）循环迭代。重复（2）～（4）步的计算，直到i 从1～Np、j 从1～n 循环完成，即完成种群的一个迭代。循环直到达到最大迭代步数或者适应值小于设定值，则中止迭代，输出当前对应的向量即为所识别的参数。

2.2.2 应力回映弹塑性算法

考虑到线弹性不能反映岩土体的材料非线性特点，在式（1）围岩位移Yi的计算中采用弹塑性有限元方法求解。塑性有限元应力和应变之间的非线性关系导致平衡方程是应变（也是节点位移）的非线性方程，变形历史取决于材料的非线性本构关系，因此，采用增量分析方法跟踪位移、应变和外部作用引起的应力的发展[10]。采用Newton Raphson 方法求解非线性方程组，采用应力回映算法实现弹性预测-塑性修正过程。弹性预测-塑性修正的步骤如下（见图2）：

图2 应力回映算法 Fig.2 The algorithm of stress return mapping

（1）弹性预测

当塑性因子 0γΔ = 时，

式中： γ ( t)为塑性因子； α ( t)为内部变量；A(t)为强化参数；Φ (t )为屈服函数；为密度；ψ ( t)为应变能函数。由式（5）～（7）得出的结果是弹性的应力状态，但还不是最终本构方程的解，如果满足式（8），即是在弹性范围内或是屈服面上，则所得到的解是方程的最后解：

（2）塑性修正

如果不满足式（8），则利用下面的式子进行计算，即返回映射方程。

式中：H 为定义硬化变量变化的综合硬化模量；N为流动矢量函数；φ 为标量屈服函数。此时也要满足 0γΔ ＞ 的要求，最终得出解。

3 基于VTK 可视化平台开发

基于DE-FEM 算法，采用微软VC++研制了隧道反分析可视化软件系统（Engineering Emulation Optimization System，EEOS）。该平台分为4 个模块，即前处理模块、FEM 计算模块、DE 算法模块和后处理模块，各模块按照功能要求采用面向对象的编程技术封装设计了相应的类，见图3。前处理模块包括三角单元剖分类Delaunaymesh，栅格法四边形单元剖分类Quadrilateralmesh 和Ansys 网格输入类，实现基本模型图形的网格剖分和复杂模型网格的导入。FEM 计算模块包括区域类（Domain）、单元（Element）类、节点（Node）类、边界（Boundary）类、材料（Material）类、时步（Time step）类和数学求解（Math）类等，实现弹性、摩尔-库仑弹塑性的有限元计算。DE 算法模块包括种群（Population）类、进化（Evolution）类、适应函数（Fitness）类和数据输入输出（Inoutdata）类。其中适应函数类之中调用FEM 计算模块计算位移，差异迭代操作产生的新参数变量输入FEM 计算模块的材料类，以此实现DE-FEM 交互耦合。

图3 反分析可视化系统EEOS 的程序模块 Fig.3 The program module of back analysis visualizing system EEOS

能否将抽象的结果数据以可视化形式呈现给现场用户，在很大程度上决定了反分析能否在工程中得到推广。科学计算可视化已成为学术和工程界的热点问题[11－13]。直接进行可视化底层代码编写需要消耗大量时间和精力，质量也不容易保证。VTK (Visualization Toolkit)是一个开放源码、跨平台、支持平行处理的图形应用函数库，以此编写可视化程序其质量和效率有显著提高，在国外的大型研究机构如Sandia, Los Alamosn 及Livermore 国家实验室等都得到应用。

图4 EEOS 反分析可视化系统界面 Fig.4 The interface of back analysis visualizing system EEOS

本文采用VTK 来进行处理模块编写。后处理模块包括等值线（Contour）类、数据输出（Dataout）类和三维建模（3D model）类。利用VTK 的vtlPlane，vtkCutter，vtkProbeFilter 类实现云图，利用VTK 的vtkPlane，vtkContourFilter，vtkProbeFilter 类实现等值线图。采用轮廓线建模，即在二维数据切片(DataSlice)中逐一提取闭合的等值线，然后将相邻切片的等值线相连接，形成曲面网格逼近等值表面，这样便可实现三维模型。该平台界面友好易用，采用了分割窗口的界面技术，界面窗口左侧为树状和按钮式控制栏，右侧为显示窗体。底部为状态栏分析过程状态信息，如图4 所示。

4 工程应用

4.1 工程概况

陕西省延安市的某隧道全长2 165 m，范围内地表之处三面高，中部低，地层主要为第四系上更新统风积砂质黄土、中更新统黏质黄土和下伏的侏罗系砂岩夹页岩。地表第四系覆盖层较厚，地势陡峻，地形起伏较大。隧道洞身通过主要地层为风化砂岩夹页岩，岩质较差，产状近于水平，节理较发育，成份以石英、长石、云母为主。隧道设计为城门洞型，高为8 m，宽为7.5 m。隧道施工过程进行了围岩的监控量测，工程地质剖面和监测断面收敛布置见图5，其中典型监测断面KD+500、KD+550和KD+600 对应的埋深分别为76、110、146 m，各断面测线的监测值见表1。

图5 隧道的地质剖面和监测断面 Fig.5 Geologic section and monitoring sections of tunnel

表1 主要断面监测数据 Table 1 The monitoring data of main sections

4.2 参数反分析

调用本文EEOS 平台分别对KD+500、KD+550和KD+600 断面相关参数进行反分析，采用摩尔-库仑的本构模型，对弹性模量E、黏聚力c 和内摩擦角φ 进行识别。围岩区域之外的岩层对隧道位移和稳定性影响不大，各断面采用均质模型。初始垂直应力按照自重应力场考虑，模型区域左右边界分别距离隧道中心为22 m，上边界为自由地表，下边界距离隧道为8 m。泊松比v 按照先验信息取为0.25。根据勘测资料，待反演参数围岩弹性模量E、黏聚力c 和内摩擦角φ 的取值范围分别为：100 MPa ≤E≤800 MPa、0.1 MPa≤c≤1.6 MPa、15°≤ φ ≤45°。3 个断面反分析参数的计算位移结果与监测位移对比见表2。计算位移与监测位移较为接近，最大相对误差为4.76%，反分析精度令人满意。

表2 各断面监测信息与反演结果 Table 2 The monitoring and back analysis results of each section

4.3 DE 参数对反分析影响

在EEOS 差异进化算法进行参数反分析过程中，CR 和F 取值对差异进化算法收敛有不同的影响，一般说来，CR 和F 取值都宜大于0.5。以KD+600断面为例，变化差异进化算法的参数，缩放因子F为0.7，杂交概率常数CR 分别为0.9、0.8 和0.7，研究优化算法的收敛性能，所对应的收敛曲线见图6。这3 种情况下CR=0.9 所对应的收敛性能最好，在第9 代适应值就可收敛；CR=0.7 收敛性最差，在第15 代才收敛。总体上这些参数都能保证算法收敛，较快地搜索到力学参数值，体现了差异进化算法具有的良好的寻优能力。

图6 不同参数下DE 收敛曲线（KD+600） Fig.6 The DE convergence curves with different parameters(KD+600)

进一步将上述问题的差异进化与遗传算法搜索过程(GA)进行比较（种群规模和迭代次数相同），收敛曲线的对比见图7。由图可见，差异进化的收敛速度和参数识别精度均优于遗传算法。基于EEOS 反演的参数计算KD+600 断面，水平位移等值线云图见图8，可视化效果良好，工程人员可直观地获得该断面水平位移分布规律。

图7 GA 与DE 收敛曲线对比（KD+600） Fig.7 The comparisons between the convergence curves of GA and DE

图8 KD+600 断面水平位移云图 Fig.8 The horizontal displacement contour map of KD+600 section

BMP 采用弹性的假设条件，同济曙光软件采用遗传算法进行参数反分析。本文在吸收上述反分析软件优点的基础上，引入优化能力更好的差异进化算法、收敛性更好的弹塑性应力回映算法，而且采用VTK 类库实现计算结果可视化，使得反分析程序具有更好的易用性和可实现性。该软件上手简单，现场技术人员容易掌握使用，采用一般PC 机求解一般规模隧道工程（1 000 个单元）反分析问题，需要5～6 h，可以满足工程现场实时反馈的要求。

5 结 论

（1）差异进化算法是具有全局优化特性的计算智能算法，解决传统优化求导数困难以及容易限于局部最小化的问题。与线弹性模型相比，本文应力回映的弹塑性有限元方法能够反映岩土材料非线性，收敛快速。基于DE-FEM 隧道围岩反分析能获得较高的反分析精度。

（2）科学计算可视化使结果数据直观表达，是计算方法在工程现场推广应用的重要决定因素。本文采用主流的VTK 可视化开发技术，实现了高质量反分析可视化平台，解决了一般反分析的输入、输出数据繁杂、抽象的问题。VTK 可视化开发方法加快了软件的设计和实现，降低了系统开发难度，减少了开发时间。

（3）初步工程应用表明，本文方法和系统是成功的。可视化平台界面友好、结果表达直观，方便现场技术人员使用。本文方法和思路为类似的工程软件开发提供了参考，对促进隧道反分析技术推广应用、推进国产软件发展都具有积极的作用。

[1] 杨林德. 岩土工程问题的反演理论与工程实践[M]. 北京: 科学出版社, 1996.

[2] SAKURAI S, TAKEUCHI K. Back analysis of measured displacements of tunnels[J]. Rock Mechanics and Rock Engineering, 1983, 16: 173－180.

[3] PIERPAOLO ORESTE. Back analysis techniques for the improvement of the understanding of rock in underground constructions[J]. Tunneling and Underground Space Technology, 2005, (20): 7－21.

[4] 文建华, 吴代华, 陈军明, 等. 地下洞室黏弹性位移反分析模式分层运算[J]. 岩土力学, 2010, 31(3): 967－970. WEN Jian-hua, WU Dai-hua, CHEN Jun-ming, et al. Parameter inversion of viscoelastic cavern displacements based on hierarchical pattern search[J]. Rock and Soil Mechanics, 2010, 31(3): 967－970.

[5] 张志增, 李仲奎. 横观各向同性岩体中圆形巷道反分析的惟一性[J]. 岩土力学, 2011, 32(7): 2066－2072. ZHANG Zhi-zeng, LI Zhong-kui, Uniqueness of displacement back analysis of a circular tunnel in transversely isotropic rock mass[J]. Rock and Soil Mechanics, 2011, 32(7): 2066－2072.

[6] 冯夏庭. 智能岩石力学导论[M]. 北京: 科学出版社, 2000.

[7] 李世辉, 隧道支护设计新论[M]. 北京: 科学出版社, 1999.

[8] 张鲁渝, 欧阳小秀, 郑颖人. 国内岩土边坡稳定分析软件面临的问题及几点思考[J]. 岩石力学与工程学报, 2003, 22(1): 166－169. ZHANG Lu-yu, OUYANG Xiao-xiu, ZHENG Ying-ren, Problems and thoughts of development of slope stability analysis software in China[J]. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(1): 166－169.

[9] JARMO ILONEN, AMARAINEN JOWI-KRISTIAN, LAMPIEN J. Differential evolution training algorithm for feed-forward neural networks[J]. Neural Processing Letters, 2003, 17(1): 93－105.

[10] NETO EA DE SOUZA, PERI'C D, OWEN DRJ. Computational methods for plasticity: theory and applications[M]. [S. l.]: John Wilen & Sons, Inc., 2008.

[11] PUSCH R. Practical visualization of rock structure[J]. Engineering Geology, 1998, 49: 231－236.

[12] 柴贺军, 黄地龙, 黄润秋, 等. 岩体结构三维可视化模型研究进展[J]. 地球科学进展, 2001, 11(4): 55－59. CHAI He-jun, HUANG Di-long, HUANG Run-qiu, et al. New progress in the study of rock structure 3-D visualization model[J]. Advance in Earth Sciences, 2001, 11(4): 55－59.

[13] WILLIAM J S. The visualization toolkit users guide, version 4.0[M]. New York: Kitware Inc, 2001.