小样本数据下三维地应力反演分析

陈正林， 何国志， 张劼超， 刘加柱， 王红娟， 候圣均

(1.中电建路桥集团有限公司， 北京 100048； 2.北京科技大学土木与资源工程学院， 北京 100083)

对于隧道和地下工程，岩体的地应力场为受力状态分析和支护加固提供数据支持，在了解工程区域原岩应力状态的前提下，可以提前规划工程设计和施工方案，有效地规避变形、坍塌等常见地质灾害，减少人员伤亡和经济损失。

为了了解一个区域的地应力分布状态，常采用地应力反演的方法，即通过区域内某几个实测应力的数据和数学推演计算来推断区域内三维地应力场。现场测量方法一般采用水压致裂法、应力解除法、声发射法、超声波谱法和发射型同位素法等[1-2]，本文研究现场测量采用空心包体应力解除法。关于获取到实测数据后，采用的地应力反演方法，1983年，郭怀志等[3]和余云燕等[4]提出多元回归分析法，同时考虑由埋深确定的岩体上覆重力和依据实际工程地质条件模拟的应力，在二者的基础上构建了不同工况下应力场。方明礼等[5]大幅改进了地应力函数反分析法，提出用三维有限元反演地应力场方法来解决地应力计算。1999年尚岳全[6]提出直接边界调整法，通过不断调整计算区域边界上的荷载大小去分析研究，使用数值方法来运算获取应力场的分布变化规律[3-6]。21世纪初，随着计算机的飞速发展，人工智能方法[7-8]逐渐被应用在地应力反演领域。

为了克服某些工程区域地质情况复杂，实测应力数据较少问题和减少现场测量的耗费，现减少测点，提出一种边界荷载和人工神经网络的联合反演方法，在小样本下对区域内地应力状态进行计算，综合两者的优点，并对建(个)元高速五老峰隧道部分区域进行地应力反演，在测点较少的情况下，利用人工神经网络的数据处理和非线性优势模拟地应力分布状态，希望可以给其他相关工程提供一定的理论依据和指导。

1 反演理论及实现

1.1 边界荷载法

在研究矿区、隧道、边坡、硐室等区域的构造应力场时，直接边界调整法具有广泛的应用性，可以很好地反演待演区域的初始地应力场和边界条件。直接边界调整法的具体步骤如下。

(1)根据待演区域的地勘报告、地形地貌特征、物探报告等地质信息，结合三维建模软件建立包括实际测点在内的区域三维地质力学模型。

(2)选定一组边界条件进行有限单元计算，监测记录测点的计算应力结果，与现场实测应力进行对比，分析出入原因，结合反演区域地质条件和模型本身，进行边界条件的大小和作用方式的调整。

(3)对三维模型重新施加调整后的边界条件和岩石力学参数，再次进行软件计算，然后再进行对比，循环往复，直至得到一组最佳边界条件使得测点的计算结果和实测结果相近[9]。

1.2 BP神经网络

反向传播神经网络(BP)属于前向型网络结构的一种，全称为“前向多层误差逆向传播学习算法”，下文简称为BP网络。这种人工神经网络通过反复改变网络拓扑结构中的各神经元之间的连接权值和阈值，使网络的输出不断逼近预设的期望输出。相比于其他网络结构，其最大的特点就是BP网络是将网络输出与期望输出之间的误差通过网络结构一层一层反向传输，各层神经元结构通过返回的误差值进行权值何阈值的调整，这种算法因此也被称为反向学习算法。

BP网络的结构和连接方式如图1所示，主要结构由输入、输出和隐含层组成。每一层有大量的节点，相邻层节点之间进行全连接，相隔层节点之间则无连接。网络的学习方式也就是网络对数据的处理方式，在BP网络里，外界的数据按照从左至右的顺序，即从输入层经过隐含层最后在输出层得到网络输出[10-11]。

图1 BP网络结构示意图Fig.1 Schematic diagram of BP network structure

1.3 地应力反演因素分析

一个地区的初始应力场的形成，往往经历了几百甚至上千万年的地质运动和外界因素的影响，且时刻在变化，不过对于某个地质年代来说，地应力场可以认为是稳定的。地应力场在形成过程中，影响因素众多，不同的区域，决定性因素和主导因素也不尽相同，研究人员也并没有定论。但是经过研究人员多年的现场考证和实验理论研究，证明初始应力场主要受以下几个因素影响。

(1)重力因素。岩体的自重应力是引起初始应力场中垂直应力的主要原因，区域任意一点的竖直应力大约等于单位面积上的上覆岩层的重量。

(2)地壳构造运动。地壳的构造运动是三维应力场中水平应力形成的主要原因，中国的大陆板块由于受四周的板块挤压和约束，内部主要产生水平的挤压应力。

(3)岩石力学性质。有学者研究证明，应力场的大小和形成与区域岩石的力学性质有一定关系，如区域岩层的弹性模量、泊松比、刚度等因素[12]。

(4)地形因素。数万年的地球活动造成了现在的山川大河，地表的高低起伏，不同的形地貌，地形对地应力的影响并不确定，不同地形处不同深度的地应力状况也并无规律可言。

(5)其他因素。除了以上对地应力影响较为明显的因素外，温度、渗流、冰川以及某些暂时不明的因素在某些特定情况下也能对区域内应力场产生一定的影响。

1.4 边界荷载法-BP神经网络联合反演

对于以上众多地应力影响因素，可以用式(1)[13]表示：

σ=F(P,R,ΔG,ΔH,T,W,Q)

(1)

式(1)中：σ为区域某点地应力；F为一种非线性函数关系；P为测点的坐标位置，可在现场测量获得；R为区域的岩石力学参数；ΔG为重力因素；ΔH为构造运动因素；T、W、Q为温度、渗流和其他暂时不明因素。

通过式(1)，将地应力表示成重力、构造应力、温度等因素的某种函数关系。在实际地应力反演过程中，测点的坐标位置是不可改变的，温度、渗流和其他的不明因素也不考虑。所以式(1)可以进一步简化，得到

σ=F(R,ΔG,ΔH)

(2)

在地应力场多元线性回归反演理论中，初始地应力场可以简化为各种因素产生的应力场的线性叠加，最终上述地应力函数可简化为

σ=C0+C1σG+C2σH+C3σR+…+Cnσn

(3)

式(3)中：C0～Cn为回归系数；σG为由在重力因素下产生的应力场；σH为在构造应力下产生的应力场；σR为由于岩石力学性质导致的应力场；σn为其他因素产生的应力场。

但是在实际上，一个区域内的地应力场是复杂多变的，尤其是在埋深较大或者地质情况较为复杂的地质区域，地应力场并不是简单地叠加，而是具有明显的非线性特征。因此引入BP网络技术，如何利用BP网络处理非线性问题的特性来解决地应力场的反演是需要解决的问题，但是BP网络需要大量样本训练网络和实际地应力测量数据稀少又构成了显著的矛盾。

首先，地应力的数值反演基本上是在有限元软件上进行，基本思路就是建立待反演区域的三维地质模型，赋予合适的参数和边界条件，使实测点的计算数据和实际测量数据相近，这样就可以得到一个区域的大致应力场。在多元线性回归反演方法中，可以通过回归后的方程推算边界条件的大小和作用方式，但是在非线性反演中，并没有一个确定的方程或者函数关系，因此就无法获得准确的边界条件。

在数据稀少的情况下，通常采用直接边界调整法，通过试算的手法来慢慢调整边界条件，但是这样非常耗费时间和精力，因此结合BP网络的非线性特征，提出以下思路。

首先建立待反演区域的三维地质模型，通过直接边界调整法进行反演计算，在试算的基础上得到大概的岩体参数和边界条件范围(实测点的地应力值被包含在内)。将试算出的岩体参数和边界条件按照正交试验构造原则，构造出若干组组合边界条件。将组合边界和岩石力学参数赋予模型进行“正算”，得到实测点的计算数据。这样就可以得到若干组BP网络的训练学习样本，以测点的地应力值作为输入，岩石力学参数边界条件作为输出。构建BP网络，带入数据进行训练学习，得到地应力值和边界条件的非线性映射关系。将测点的实际测量数据带入训练完成的网络，“反算”出最佳边界条件，代入模型中，进行模拟计算验证。具体流程如图2所示。

图2 地应力非线性反演流程Fig.2 Nonlinear inversion flow of ground stress

2 工程实例分析

2.1 五老峰隧道地应力实测

五老峰隧道位于建水至元阳段，隧道左线全长8 305 m，起点桩号Z3K22+670,终点桩号Z3K30+975;右线全长8 360 m，起点桩号K22+665,终点桩号K31+025。隧道一般埋深50～900 m，最小埋深25 m，最大埋深929 m，属典型的深埋特长分离式隧道。本次现场测量采用新型数字化无线式瞬接续采型原位空心包体应变计[14-15](图3)。由于实验环境复杂，本次只进行一个点的地应力测试，具体计算过程限于篇幅原因，不作具体介绍，测量结果如表1所示。为了便于计算，将主应力进行应力分解，结果如表2所示。

图3 包裹空心包体的岩芯Fig.3 A core enclosing a hollow inclusion

表1 主应力值

表2 实测三维应力分量

2.2 三维地质建模

通过对五老峰隧道地区的深入的地质调查，结合本次实验测点位置和模型建立原则[6]。试验点位于五老峰隧道出口段加宽道位置，桩号K27+982～932，埋深大约900 m，根据隧道图纸，确定模型y轴与隧道轴线平行，x轴垂直于隧道轴线，埋深方向为z轴。隧道区域实际位置处于模型中间位置，远离边界，实际测点位于隧道中段，同样处于模型中央位置。五老峰隧道布置图及模型计算边界范围如图4所示。

通过FLAC3D建立五老峰隧道部分区域的三维地质模型，如图5所示。模型总计165 000个网格、174 216个节点、21 100个面。此次建立三维模型既具有计算精度上的优势，又避免了冗长的计算时间。

根据地勘报告和现场的取样室内试验，得知反演区域的花岗岩部分岩石力学性质。按照室内试验结果和以往模拟经验，在参数选择上，重度取 2 650 N/m3，黏聚力0.5 MPa，内摩擦角35°，本构模型采取弹塑性岩体模型(莫尔库伦模型)。其中弹性模量和泊松比在本节不进行选取，这两个参数需要利用BP网络和FLAC进行推算，需要后期试算确定范围。

病理学主要是研究疾病的病因、发病机制、病理变化及过程，是以形态学为基础，通过病变大体标本和细微结构的观察认识，掌握疾病的本质及发生规模对许多疾病的病理特点的观察、掌握，要从宏观和微观的角度同时的观察大体标本和显微境下的特点，因此，病理学是基础医学和临床医学之间的纽带[2]。

本文在模型边界条件上主要选择自重因素和构造运动这两类主要因素施加，地形地貌在三维地表建模中得以实现，岩石力学性质则在模型赋参数这一环节得以实现。具体到模型上分为重力值、x向位移、y向位移，岩体的弹性模量和泊松比。

图4 五老峰隧道区域三维示意图Fig.4 Three dimensional schematic diagram of Wulaofeng tunnel area

图5 五老峰隧道三维地质模型Fig.5 Three dimensional geological model of Wulaofeng

2.3 边界荷载法人工试算

由于数据的稀少，本文研究只能采取人工试算的方法使得实测点的计算值逼近测量值，又因为直接边界调整法本就繁杂，要想仅通过边界调整法找到理论上最佳的反演边界条件，使得测点的计算值很好拟合测量值会耗费大量的时间和精力，当反演区域地质条件复杂，测点数量较多时，这种仅靠人工试算的方法更加不现实。所以，需要在直接边界调整法的基础上，通过人工试算确定边界条件的大致范围，使得实测数据处在这个计算范围之内，构建训练数据，然后通过BP网络反推最佳边界。

经过若干次边界调整计算，确定了两种岩体参数(岩体的弹性模量和泊松比)，3种边界条件(重力值、x向位移、y向位移)的大致范围。每个因素在波动范围内取4个水平。具体值如表3所示。

表3 边界条件水平

2.4 边界荷载-BP神经网络联合反演

2.4.1 网络训练样本构建

2.3节确定了边界条件和部分模型参数的范围并且给每个因素划分了4个水平，BP网络训练需要大量样本，所以利用这5个因素4个水平(表3)人为构造训练样本，这样一来，边界荷载法和BP神经网络的优势互补，既避免了边界荷载法大量试算的缺点，又通过前期试算给BP网络提供的样本数据。样本构造则按照正交实验法则[16]，具体如表4所示。

上述16组均匀设计组合充分地考虑了5个因素的4个水平，但是为了增加BP网络训练的可靠性和精度，另外的增加了4组随机水平组合(和上述16组不重复)，如表5所示，从而使训练样本数量达到了20组，基本满足BP网络的训练样本要求。

表4 16种边界条件组合

表5 4组随机边界组合

通过上述步骤，得到了20组有完整输入输出指标的BP神经网络训练数据，接下来就可以构造关于三维应力和边界条件的非线性网络结构，利用这些样本数据对网络进行训练。

2.4.2 网络构造及训练

网络的构造和训练则利用MATLAB软件实现，网络模型采用最高效的三层模型，具体参数如表7所示。

分析本文训练好的网络(图6和图7)，训练、验证和测试的拟合优度R都在0.9以上(R越接近1说明拟合越好)，网络训练的误差基本集中分布在0.03附近，极少数样本误差到了1或者2，但都在允许范围内。这样的结果说明本次训练的BP网络达到了期望效果，可以进行下一步的最佳边界条件反演，具有工程指导意义。

表6 BP输入、输出样本

表7 神经网络参数

2.4.3 网络构造及训练

2.4.2节中通过BP网络得到了三维应力分量和边界条件的非线性关系，虽然这种关系无法利用函数关系式表达，但是可以利用保存的网络映射输入与输出之间的关系。本文最终目的是利用训练完成的网络，将测点的真实三维应力分量作为输入带入网络，就能得到三维模型在有限差分计算软件中的最佳边界条件和部分模型岩体参数。

将表2代入训练好的网络，利用MATLAB命令a=sim(net，b)，其中a为输出(边界条件、岩石力学参数)，b为输入(三维应力分量)，net为训练完成的网络，可以得到最佳边界条件。由于BP网络的非确定性，决定进行了10次边界条件反推，具体结果如表8所示。

图6 BP网络拟合回归图Fig.6 Regression diagram of BP network fitting

通过训练好的BP网络进行了10次反推，得到了上述10组边界条件和岩石力学参数。到这一步，通过直接边界调整法和BP网络的结合成功反推出边界条件，但是这些边界条件是否真的满足反演精度要求？与实测值相比又如何？还需要将这些边界条件和岩石力学参数赋予三维地表模型，通过有限元软件计算来检验这些条件的有效性，进一步验证网络的有效性。具体的计算结果如表10所示。

从表9中可以看出，10组边界条件下三维应力计算结果与实测值总体上是相一致的。除去个别误差稍大外，自重应力(σzz)平均模拟误差1.24%，最大2.29%，最小0.09%；x向水平应力(σxx)平均模拟误差1.04%，最大3.03%，最小0.04%；y向平应力(σyy)平均模拟误差0.55%，最大1.22%，最小0.03%。其中第8组边界条件下模拟结果与实测应力很是接近，三个方向上的应力误差都控制在1%以下，y向水平应力甚至与实测应力相同，总体模拟结果精确在可接受范围内。

图7 BP网络误差分布图Fig.7 Error distribution of BP network

表8 边界条件BP网络输出

表9 计算应力与实际应力对比

2.5 五老峰隧道区域应力场分析

利用FLAC3D的切片(cutting tool)功能，将沿五老峰隧道轴线方向的三个主应力进行竖向切片分析，如图8所示。三维模型区域的隧道总长约1 400 m，走向与坐标y轴平行，最大埋深约900 m。

(1)从图8(a)可以看出，五老峰隧道区域的垂直应力总体呈梯度分布，在接近地表区域的应力等值线起伏较大，这说明工程区域的地形地貌对浅层区域的竖直应力分布是有一定影响的，这主要是地形的起伏造成同一水平下上覆岩层的自重力不尽相同，而竖直应力大部分又是由重力造成。沿隧道轴线的竖直应力分布在23～28 MPa范围内，波动较大且略大于重力，原因是因为深部构造运动强烈，与自重力共同作用下形成了竖直应力，但总体上竖向应力是符合应力分布规律的。

图8 五老峰反演区域云图Fig.8 Inversion area cloud map of Wulaofeng

(2)从图8(b)可以看出，五老峰隧道区域的最小主应力也是呈梯度分布，与竖直应力不同的是，地表的高低起伏对水平的应力的分布影响不大，只在地表约100 m的位置造成应力起伏。随着深度的增加，应力等值线越来越趋向水平，这说明五老峰区域的最小主应力主要是受水平构造运动影响，重力只会在山体表层对水平应力造成影响。沿五老峰隧道轴线最小主应力分布在28～29 MPa，应力基本在一条等值线上，方向呈水平，平均应力在数值上大于竖直应力，因为隧道埋深较大，越到深部，水平应力与竖直应力的差距就越大。

3 结论

为了解决在某些工程区域地应力测量难，实测数据缺乏的问题，采用边界荷载法和BP神经网络的联合反演法，既能解决边界荷载法人工试算时间过长，又能通过神经网络产生大量相关数据样本，最后以五老峰隧道为工程实例进行验证，得到如下结论。

(1)岩体的原位三维应力场并非是简单相关影响因素的线性叠加，在测量数据较少和地质条件复杂情况下，多元线性回归法不适用于反演研究，边界荷载-BP神经网络联合方法可以有效地克服并反演原岩地应力状态。

(2)采用人工试算边界条件范围可以有效减少单一边界荷载法的计算时间，BP神经网络利用正交实验法则提供的足量数据可以精确地反算最佳的边界条件，进而模拟三维应力状态。

(3)联合反演的三维应力模拟结果与实测相比，自重应力平均误差1.24%，x向水平应力平均误差1.04%，y向水平应力平均误差0.55%，模拟结果较为精确。