基于XGBoost的初始地应力场修正反演法

刘敬文，白金朋，王 建，徐 韧，徐伟男

(1.河海大学水利水电学院，江苏 南京 210098； 2.中国电建集团北京勘测设计研究院有限公司，北京 100024； 3. 无锡市惠山区水利局，江苏 无锡 214174)

初始地应力场是地下厂房等工程设计与施工所需考虑的重要因素之一[1]，然而，受场地、经费等现实条件的限制，不可能大规模开展原位地应力测量，只能利用有限的实测地应力样本来反演厂房区的地应力分布，这是目前确定初始地应力场的主要途径。

常用的反演方法主要有多元线性回归分析法[2-3]和智能优化算法[4-6]等，其中，多元线性回归分析法采用最小二乘法求解回归系数且解唯一，故被广泛采用[7]。但是，多元线性回归分析法是基于弹性力学的叠加原理，更适用线弹性材料的情况。对于复杂地质条件下的岩体应力场，部分岩体结构不满足弹性假设，力学非线性较为显著，需要开展修正以解决线性回归可能存在的拟合精度问题。此前已有部分学者研究过修正、优化反演的方法，郭运华等[8]在最小二乘回归的基础上，根据屈服条件修正偏应力大小，对不满足强度条件的单元地应力进行修正，以满足弹性假设条件。郭明伟等[9]采用优化位移边界条件拟合工程区域初始地应力场，通过调整位移边界模式得到其最优组合并依此拟合工程区域的初始地应力场。闫相祯等[10]利用阻尼最小二乘法建立了应力场反演的优化约束模型，提出了一种用于反演地应力场边界力的优化分析方法。然而，上述方法大都从数值模拟的角度出发，通过改变模型边界的约束参量实现修正反演的目的。为了进一步提高反演计算精度，本文在有限元多元线性回归分析的基础上，采用极端梯度提升(extreme gradient boosting,XGBoost)算法，以不同构造应力下有限元的计算结果为特征函数，以应力实测值与线性回归结果的差值为目标函数，对多元线性回归的反演结果进行修正，其核心思想是基于实测应力场与线性回归应力场的差值构造一个残差修正场，这样做一方面有效减小了使用线弹性模型回归计算带来的误差，同时避免了直接使用XGBoost纯数学算法拟合应力场，更加准确拟合实际应力场。

XGBoost算法是一种基于Boosting的集成学习算法，该算法具有计算复杂度低、运行速度快、精度高等优点。由于迭代初始值采用单位载荷法，XGBoost有效地避免了地层几何构造和边界条件的限制，同时，对不满足弹性假设的一些单元的地应力进行了修正，使其更符合实际情况。因此，使用XGBoost修正反演结果，一方面可以有效修正使用线弹性模型回归计算带来的应力偏差，提高反演结果的精度，另一方面利用该集成算法所得的黑箱模型可以有效降低误差、增强数据可靠性，从而为提高地应力场反演计算的精度和数据的可靠性提供了一种新的解决方案。

1 修正反演分析原理和方法

1.1 多元线性回归法

(1)

根据最小二乘法原理，只需求得使残差平方和为最小值时的n个待定回归系数L=(L1,L2,…,Ln)T，计算域内任一点P的回归初始地应力就可由该点各工况有限元法计算值迭加而得

(2)

式中：j=1,2,…,n对应n个初始地应力分量。

1.2 XGBoost算法修正反演原理与模型构建

XGBoost算法是梯度提升机器(gradient boosting machine，GBM)算法的扩展，是一种同时具备非线性模型和树模型特点的优化模型，能够同时完成回归和分类任务，且已在疾病诊断、交通流量预测、商品销售量预测等领域的应用取得成功[11]，目前在水利及岩土领域应用较少。XGBoost算法由多棵决策树(cart)组成，通过决策树集成实现机器学习，所有决策树的回归值相累加即为模型回归值。与一般算法不同的是，XGBoost算法通过对损失函数二阶泰勒展开以快速逼近目标函数，并加入了正则项，调节参数避免过拟合，训练速度快，拟合及预测精度高[12]。

在构建地应力回归的修正反演模型时，XGBoost算法通过在回归函数中不断加入关于各因子的新函数来逼近实测值，以解决地应力回归时非线性拟合不足的问题，即：

(3)

XGBoost算法的目标函数定义如下：

(4)

式中：C为损失函数，用于评估回归值与真实值之间的损失；Ω(fk)为正则化函数，表示第k棵决策树的复杂度，用于控制模型计算效率，避免过拟合。正则化函数定义如下：

(5)

式中：γ、λ为正则项惩罚系数；T为该树的叶子节点个数；ω为叶子节点对应的权重。

可以对式(4)改写并按二阶泰勒级数展开，进一步简化后得到：

式中：gi为损失函数的一阶梯度统计；hi为损失函数的二阶梯度统计。通过迭代找到f(x)使损失函数最小则模型训练完成。

综上所述，本次XGBoost算法修正反演的建模流程如图1所示，具体步骤如下：

图1 基于XGBoost算法修正反演的计算流程

步骤1根据地质条件，对厂房及周边区域进行有限元建模。开展基于实测地应力的反演分析，对模型进行各构造模式对应的加载和约束，获取各构造模式下的应力场，并利用多元线性回归分析反演得到整个厂房及周边区域的初始地应力分布。

步骤2提取有限元计算的应力结果和多元线性回归值的残差，对这些数据进行处理，整理成合适的影响因子，共同组成XGBoost输入样本的数据集，并划分出训练集和测试集。对模型进行参数调优，选取最优参数结果构建基于XGBoost的修正反演模型。

步骤3将测试集中的数据代入步骤2中构建的模型，并进行交叉验证，分别比较多元线性回归值、模型修正反演值与实测值的均方根误差RMSE和相对误差均值RSE，对XGBoost修正反演模型做出评价。最终利用该模型对厂房及周边区域内的初始地应力回归值进行修正，得到修正后的整体初始地应力场。

2 工程实例

2.1 工程概况及地质条件

某抽水蓄能电站位于辽宁省大连市境内，其上水库为东、北、西三面环山的天然库盆地形，主坝修筑于库盆南侧的主沟内，副坝位于库盆东北侧低矮分水岭处，均为沥青混凝土面板堆石坝；输水发电系统深埋于上、下水库间山体内，围岩以微风化～新鲜变质长石砂岩为主。

在进行地应力回归分析时，计算模型的范围为：上、下游边界以地下厂房为中心，向南北方向各延伸400 m；左、右边界以地下厂房为中心，向东西方向各延伸500 m；上、下边界上至地表，下至高程-170 m。

初始地应力场反演区域和三维数值模型如图2所示。模型主要考虑fp25、fp28、fp31、fp32等较大的断层，有限元模型共有节点169 349个，单元967 079个。本次计算将断层简化为匀质连续体，采用弹性模型模拟岩体及断层的力学行为，涉及的力学参数包括弹性模量、泊松比和密度等。本次计算采用的岩体和断层力学参数见表1。

图2 地应力反演三维有限元模型

某抽水蓄能电站在地下厂房斜井段(钻孔ZK213)、高压岔管段(钻孔ZK201)、下平段(钻孔ZK204)、厂房区(钻孔ZK202、ZK206)、地表孔ZK149共选定6个钻孔(每个钻孔有5个测点)内进行了水压致裂法地应力测试。钻孔位置见图3，地应力

表1 计算岩体和断层力学参数

试验实测结果汇总表见表2，表中SH为水平方向最大压应力；Sh为水平方向最小压应力；SV为竖向应力；负号表示为压力。

图3 地应力测孔位置示意图

表2 实测地应力

2.2 基于多元线性回归的一次反演

根据实测地应力的结果，采用多元线性回归进行初始地应力场的反演，通过模拟若干组构造模式，并对其进行线性组合，从而获得最终地应力。由于水压致裂法无法测得铅锤面内的剪切应力，一般采取以下5种构造模式：①自重；②x方向均匀分布、三角形分布的挤压构造力；③y方向均匀分布、三角形分布的挤压构造力；④xy向的剪切应力；⑤yx向的剪切应力。由这5种构造模式进行多元线性回归计算，以确定自重系数L1，x方向水平构造系数L2、L3，y方向水平构造系数L4、L5，xy向剪切应力系数L6、yx向剪应力系数L7。将32个测点的192组应力数据代入公式进行计算，求得回归系数结果如下：L1=1.189、L2=1.498、L3=0.553、L4=1.322、L5=-0.526、L6=1.652、L7=-1.470。复相关系数r=0.974，回归效果较好。偏相关系数V1=0.965，V2=0.915，V3=0.364，V4=0.940，V5=0.440，V6=0.701，V7=0.752。其中V1、V2、V4较大，说明所施加的自重应力场、x方向均布荷载、y方向均布荷载作用较显著。

求得回归系数后，代入式(2)求得在实测位置的回归计算值，部分回归结果与实测值的比较见表3(表中数据压应力为负，拉应力为正)。虽然本次线性回归反演的复相关系数较高，但标准差却高达0.97 MPa，说明反演拟合精度存在一定的不足，需要进行修正。

表3 部分测点的应力分量实测值与回归值对比 单位：MPa

2.3 基于XGBoost的修正反演

2.3.1数据预处理

将2.2中的每一种构造模式代入有限元，作为边界条件进行正演计算，得到有限元计算结果，并提取L1～L7七个系数分别对应工况下的应力因子，作为输入变量；将2.2中多元线性回归结果中应力分量的残差作为输出变量。

输入样本总数为180组，将输入样本集的前150组划分为训练集，后30组划分为测试集。

2.3.2模型的参数调优

对各关键参数进行初值定义，由于模型参数max_depth(最大树深)对突变数据较为敏感，对模型收敛的稳定性影响较大，因此针对max_depth进行敏感性分析，以权衡参数的最优取值。

将max_depth的优化范围设为(2，8)，比较和权衡不同树深取值下模型的收敛速度与精度，调参结果显示，当最大树深max_depth=4时，模型的收敛速度与精度都较高，收敛稳定性较好，故本次最大树深参数选定为4，其余参数为默认值。

模型设置最大迭代次数为250次，最多经过190次迭代即可收敛。图4为模型训练迭代过程。

图4 模型训练迭代过程

2.3.3模型的结果与评价

求得测试集的修正反演值后，与实测值、回归值的对比如表4所示。

由表4可知：多元线性回归结果与实测值之间的均方根误差为1.13 MPa，XGBoost修正反演结果与实测值之间的均方根误差为0.48 MPa，拟合精度提高约60%；从相对误差的平均值看，XGBoost修正反演精度也得到不同程度提高。可见XGBoost修正反演的结果更精确且更具有代表性，对多元线性回归的反演结果进行了二次优化，同时，也在一定程度上减少了非协调回归因子的影响，在多元线性回归的基础上提升了反演精度以及数据的可靠性。

为探究本文方法在工程实际中样本量相对较少时的反演精度，以及评估模型的质量、避免过拟合与欠拟合现象的发生，采用两种方法加以验证：

a.保持原先30组测试集不变，将训练集数量由150组减少到60组，其余参数不变，得到修正反演结果与实测值之间的均方根误差为0.54 MPa，相比于多元线性回归结果的均方根误差1.13 MPa，拟合精度提高约50%。由此可见，当样本数量较少时，该修正反演方法的精度依旧较高，效果良好。

表4 钻孔测试集主应力的反演分析值对比

b.保持训练集与测试集总数不变，重新划分一组数据集，进行交叉验证。取样本集的后150组为训练集，前30组为测试集，其余参数不变，得到修正反演结果与实测值之间的均方根误差为0.65 MPa，相比于多元线性回归结果的均方根误差1.56 MPa，拟合精度提高约58%，说明模型稳定性高、泛化能力好，能够很好地满足初始地应力修正反演的工作需求。

本文模型基于Python3.7进行，计算机内存12.0 GB、4核CPU，计算967 079个单元的初始应力共计用时948.5 s，反演效率高，满足工程需要。

图5和图6分别为主厂房轴线、垂直主厂房轴线剖面的第一、第三主应力等值线云图。由图5和图6可知，主厂房开挖区域的主应力从上到下基本递增，等值线在浅层受地形变化影响较大，说明修正反演计算的应力场充分反映了地形、地貌的影响。

图5 主厂房轴线剖面第一主应力云图

图6 主厂房轴线剖面第三主应力云图

3 结 语

在地应力反演的传统多元线性回归方法基础上，引入XGBoost算法并基于实测应力场与线性回归应力场的差值构造残差修正场，对多元线性回归的反演结果进行了二次优化，既提高了反演精度，同时也避免了直接使用纯数学的XGBoost算法来拟合应力场物理概念不明确、容易产生过拟合的问题，更加符合力学原理。

工程实例计算表明，采用本文方法，拟合精度提高约60%，且计算稳定性高、泛化能力好，能够很好地满足初始地应力修正反演的需求。