基于应力函数法的初始地应力场反演研究综述

2020-09-18王呈璋

四川建筑 2020年4期

王呈璋

(山东理工大学交通与车辆工程学院,山东淄博 255000)

1 初始地应力场反演

随着调水工程输水隧洞、深埋长距离交通隧道、深部矿产能源等大型地下工程的不断开发建设[1-5]，越来越多的工程面临较为复杂的地质构造环境，而工程区域初始地应力场的分布特征受区域地质构造的直接影响。如何能够准确地获得工程所需的初始地应力场量值与分布规律，一直都是地下工程研究的基本问题之一。在地下工程的设计与施工中，地应力都被视为最基础的荷载因素，工程设计是否有效可靠，施工过程是否合理安全，均受其直接的影响[6-8]。

可靠的地应力实测数据可以为工程建设提供有效的参考[9]，但是由于工程区域复杂环境的限制，经费的要求以及试验技术的局限性，无法开展大量的地应力实测工作，而有限离散的地应力实测数据往往无法满足实际工程的需要[10-11]。因此，在地应力实测资料的基础上，结合地质构造背景、地形地貌特征、岩性分布规律等地质资料，以数学力学方法为数值计算的手段，反演出有效的工程区域初始地应力场，是初始地应力场反演的基本思路[12]。在众多的初始地应力场反演方法中，比较具有代表性的有应力函数法、侧压系数法、边界荷载法、边界位移法、神经网络法、遗传算法、变差函数法、径向基函数法等[13-17]，其中应力函数法是较为经典的初始地应力场反演方法。应力函数法在工程实践过程中，被不断的改进发展，最初的应力函数法较为单一，仅仅为根据弹性理论设定的多项式，通过对多项式的系数求解，进而求解出应力函数的具体形式。随着地应力反演计算相关理论的发展，以及计算条件和计算手段的进步，有关学者通过构建更为复杂多样的应力函数，如三维正交多项式函数，应力趋势函数等，利用计算机高效的计算效率进行初始地应力场的反演。近些年来，伴随着大型通用计算软件的发展，应力函数的形式也变得多样化，应力函数从最初的反演目标，逐渐转变成为反演手段，将应力函数与有限元联合进行初始地应力场反演，成为目前较为主流的反演方法之一。

基于相关学者在初始地应力场反演过程中对于应力函数法的研究成果，本文从基本思想、实现过程、工程应用特点和优缺点等方面，着重介绍了该方法在不同阶段的重要特征。

2 应力函数法的基本理论

初始地应力场的应力函数法回归反演分析，即根据有限且离散的地应力实测数据，结合工程区域的地质计算模型，以数学力学方法为主要手段，构建三维应力场函数或应力分量函数的表达式，利用数值计算方法，推导出工程区域的初始地应力场[18]。在初始地应力场反演分析的方法中，应力函数法的基本概念提出的时间较早，经过大量学者的丰富研究，现在已经具有多样的形式，如多元线性回归分析法、应力函数趋势分析法、应力函数数值分析法、边界荷载分析法等，以上的方法均采用不同的应力函数形式作为数值计算方法的核心。

2.1 多元线性回归分析法

数值计算方法多种多样，其中较为典型的应力函数法最早是由郭怀志等[19]提出，他认为岩体在弹性工作状态下，初始应力场可视为因变量，若干自变量共同作用形成应力场这一复杂的因变量，这些自变量分别为弹性模量、泊松比、容重、地温等物理因素，以及地质构造特征、地形地貌特征、构造作用下岩体的三维坐标等构造因素，各种因素在地质体内产生的应力视为基本初始应力，而初始应力场即为各基本初始应力的线性叠加，其表达公式为

(1)

基于这种回归反演分析理论，初始地应力场的反演问题即转变为对回归系数bi的求解。具体实现过程为，首先将各自变量因素作用于地质计算模型，得到对应状态下的基本初始应力计算值，根据线性叠加原理求得初始地应力场的计算值，结合地应力实测数据，将二者的差值利用最小二乘法原理使残差的平方和最小，即可得到回归系数bi的最优解。

2.2 应力函数趋势分析法

张有天[20]等提出了一种根据趋势对地应力场进行分析的应力函数法，即地应力场的趋势分析法。该方法的核心思想是假设高阶的多项式应力函数，结合有限的地应力实测数据以及边界约束条件，根据趋势对地应力场进行反演分析。该方法首要的任务是假定合理的四次多项式函数φ1、φ2、φ3，三个应力函数应该满足Beltrami-Michell方程，应力分量由应力函数表示的公式如下：

σx=∂2φ3/∂y2+∂2φ2/∂z2

σy=∂2φ1/∂z2+∂2φ3/∂x2

σz=∂2φ2/∂x2+∂2φ1/∂y2

τyz=-∂2φ1/∂z∂y

τzx=-∂2φ2/∂z∂x

τxy=-∂2φ3/∂x∂y

(2)

由式(2)，六个应力分量σx、σy、σz、τxy、τyz、τzx均可表示为关于空间坐标(x、y、z)的二次多项式。对于地应力实测数据，每一个测点都包含六个应力分量，当实测点为n个时，可以建立6n个独立的方程。作为边界约束条件的地表点，由于其法向应力与剪切应力均为零，因此可以选取若干地表点建立方程组。当六个应力分量的多项式系数个数与所有方程的个数相同时，即可求得六个应力分量的具体表达式。

2.3 应力函数数值分析法

肖明等[21]在对锦屏二级水电站坝址处山体进行初始地应力场反演回归分析时，提出了三维应力函数拟合方法。该方法主要分为两个部分，首先根据工程区域范围内地勘资料所呈现的地质体地表剥蚀风化产生的地貌特征、时空变化下的构造演化特征和岩土物理力学性质，采用有限元反演出初始地应力场。然后根据测点位置的反演计算值和实测值，运用三维有限元对计算后的地质体模型进行插值运算，利用三维正交多项式构建出能够合理反映三维应力场的应力函数。假定一组多项式序列p0(x,y,z)，p1(x,y,z)，...，pm(x,y,z)满足以下的正交条件：

(3)

式中：n为拟合区域离散点的个数；当i=j时，p为常数；当i≠j时，p为零。

由正交多项式的性质构造的正交多项式系列如式(4)所示：

p0=1

p1=(z-a1)p0

p2=(z-a2)p1+(x-b2)p0

⋮

pi+1=(z-ai+1)pi+(x-bi+1)pi-1+(y-ci+1)pi-2+di+1pi-3

(4)

式中：ai+1、bi+1、ci+1、di+1为待定系数。

由数学归纳法求得空间坐标(x、y、z)表示的多项式系数和多项式pi的具体表达形式，从而拟合出三维应力函数的具体表达式：

Qi=αi0p0+αi1p1+…+αimpm

(5)

式中：αi0，αi1，…，αim为待定的应力系数。

通过建立以地应力实测数据与拟合边界的应力值为极值条件的目标函数，利用最小二乘法的原理对目标函数求极值，即可得到式(5)中的应力系数，进而得到三维应力场的应力函数表达式。

2.4 边界荷载反分析法

边界荷载反分析法的核心思想是地质体模型的边界荷载与初始地应力状态是相对应的关系，如果已知地质体模型的真实边界荷载，则可得到初始地应力场。因此，边界荷载法的实现过程就是通过不断地调整地质体模型的边界荷载，使得地应力测点位置处的应力回归计算值与实测值之间的误差达到最小。

刘允芳[22]等对于该方法的具体实现过程为，根据地质勘测资料，选择合理的工程区域范围，构建地质力学计算模型。对时空变化下的构造演化过程进行分析，确定可能对工程区域地应力场的形成产生影响的主要地质构造因素，一般选择自重应力、两个水平方向的构造应力和水平面内的剪切应力，将以上四种因素作为待定因素。对于每一种待定因素，将其作用于地质计算模型，求得在已知测点位置处的地应力值。将测点位置处所有计算求得的地应力值与现场实测的地应力值之间建立回归方程：

(6)

利用最小二乘法原理使得地应力测点的实测值与回归计算值之间的残差达到最小值，得到对应于四个待定因素的回归系数，即可得到应力函数的回归方程。

3 反演方法的工程实践特征分析

3.1 初始地应力场的区域特征

由于地应力场具有时间性和空间性二重属性，所以就整体而言，地应力场是非线性且非稳定的应力场[23]。在复杂漫长的地质历史年代，地应力场受多重因素的作用，如地壳的运动、地温的不均匀变化、地幔的随机对流、地表的风化剥蚀和地形地貌的变化、岩体自重的影响、岩体内孔隙水压力以及节理结构面的影响等[24-25]，如果在地应力场反演工作中全面考虑这些因素的影响，理论上是可以以较高的精度反演出与实际情况相符合的地应力场，但是在实际操作中可行性较低。一方面，地应力场的非线性特征很难用单一的非线性函数或者非线性模型进行模拟表征。另一方面，地应力场又是整体上的非稳定场，无法准确的表示出地质历史年代整个应力场的变化过程。但是，在研究具体的工程问题时，需要考虑的地质变化过程一般较短，而在所考虑的时间尺度和空间尺度的范围内，地质体是相对稳定的状态。所以，在对具体的工程问题进行初始地应力场反演时，完全可以对各种因素进行简化处理，将工程地质体视为相对稳定的非线性模型。

3.2 反演方法的工程应用特征及优缺点

应力函数法就是在这样的背景下提出来的，不同形式的应力函数均可以在一定程度上模拟地质体的非线性特征。应力函数趋势分析法是最为直接的方法，该方法用四次的函数来模拟三维的应力场，基于弹性力学的假定，根据平衡方程、物理方程、相容方程来确定应力分量的具体表达式，对于地质体模型中的任意一点，均可由其坐标值得到对应的应力值，从而得到整个过程区域的初始地应力场。该方法所需的计算工作量较小，简便易行，具有较强的可操作性，表示的是连续空间的地应力场。缺点在于将地质体视为简单的弹性模型，并未考虑岩体的塑性变形。而且地质体内部存在的断层、褶皱等地质构造特征均无法充分的考虑，所以该方法的代表性不强，反演精度较低，适合在岩性比较均匀，断层发育较少、褶皱等构造地质构造作用不够强烈的地壳深部工程使用，如王拉才[26]等在抚顺龙凤煤矿、孔广亚[27]等在某黄金矿的的地应力场反演中均采用该方法，取得了良好的效果。

应力函数数值分析法从地应力实测点出发，虽然离散且有限的应力实测数据无法满足三维有限元计算的要求，但是完全可以利用正交多项式的特性，根据实测数据构建正交多项式序列，并且结合插值运算拟合出一个应力函数。该方法在锦屏二级水电站地应力场反演中首次运用，并且分别考虑了存在断层影响和不存在断层影响两种情况。反演结果表明采用此方法反演回归的初始地应力场在两种情况下的结果基本一致，仅仅在靠近断层附近的主应力值有所差别，说明改方法具有一定的精度，并且便于操作。但是，该方法的局限性在于建立三维地质体模型的过程中对地质特征进行了较多的简化处理，对于地形地貌特征以及构造演化特征较为复杂的地质体，一个拟合的应力函数无法准确描述整个初始应力场的真实状态。赵辰[28]等以该方法为基础，保留其核心的思想，即构建正交多项式系列，通过正交多项式序列拟合出函数的具体形式，从而提出了考虑地层剥蚀作用的侧压力系数法这一地应力场反演方法。该方法拟合的函数不再是应力函数，而是侧压力系数函数。对于复杂的地质构造环境而言，初始地应力场的变化可能比较复杂，但是在一个相对稳定的区域内，侧压力系数这一代表水平应力与自重应力比值的量值保持在稳定的范围内，因此完全可以通过建立侧压力系数拟合函数来进行地应力场的反演。利用有限元软件通过逐次开挖来模拟地表风化剥蚀的过程对初始地应力场的影响，以应力实测点的侧压力系数数据为基础，采用正交多项式序列构造出三维空间的侧压力系数拟合函数，由水平应力、自重应力和侧压力系数之间的数学关系，即可得到整个地质模型的初始地应力场。该方法的显著优点就是可以充分考虑地形地貌特征、地质构造条件对初始地应力场的影响，反演方法具有一般性，适用于各种复杂地质体，反演结果具有较高的精度和合理性。

多元线性回归分析初始地应力场的方法，由于考虑多种影响地应力场的因素，并且每一种因素都对应着唯一的应力状态，将每一种应力状态视为一个应力函数，因此通过线性叠加得到的对应于三维空间初始地应力场的应力函数，具有唯一性[29]。边界荷载反分析法不仅具备多元线性回归分析法的优点，而且其独特的优点就是仅考虑对初始地应力场产生显著影响的四个因素，包括自重应力因素、水平方向两个构造应力因素以及水平面内的一个剪切构造应力因素，反演过程只需求得各因素作用于有限元地质模型的最优系数，即可获得各因素共同作用下的唯一初始地应力场。该方法最早应用于水利水电行业，由于其采用有限元计算软件建立逼真的工程地质模型，能够反映地形地貌特征、断层和褶皱等地质构造特征，并且采用科学可靠的数理统计理论做计算支撑，因此从提出来之后便被行业广泛认可并普遍使用，如抽水蓄能电站[30]、引水隧洞工程[31]、水电站高边坡[32]等大型水利水电工程的初始地应力场反演均采用该方法。虽然在其他工程领域，利用多元线性回归分析法和边界荷载反分析法对工程区域初始地应力场反演的研究工程相对较少，但是由于一些交通隧道工程正在面临超长深埋等较为突出的地质问题，利用上述方法进行地应力场反演的优势也正在凸显。汪波[33]、代聪[34]和尤哲敏[35]等分别在在苍岭特长公路隧道、蓝家岩深埋特长公路隧道和大坪山深埋公路隧道的初始地应力场反演分析中运用多元回归分析原理，获得了较为合理可靠的隧道轴线处初始地应力场分布规律。因此，多元线性回归分析法，以及以此为基础发展而来的边界荷载法，由于能够更为真实的反应工程区域的地质构造特征，更为合理可靠的反演出初始地应力场的分布规律，无论是规模较大的水利水电工程，还是深埋超长的交通隧道工程，都将其视为较为科学的方法并成功应用于工程实践。

4 结论与展望

4.1 结论

本文以应力函数法为基础，详细介绍了该方法的提出背景，系统地总结了该方法在具体的工程实践中发展而来的多元线性回归分析法、应力函数趋势分析法、应力函数数值分析法和边界荷载反分析法，具体阐述了每一种方法的基本原理和数值计算方法，对比分析了各方法的工程应用特点以及优缺点，具体对比结果如表1所示。主要结论如下：