考虑非均质性和物理信息的三维原位应力预测在深煤层中的应用 李冰、 康云伟、 郝鹏程、 朱伟平、 何鹏博、 甄怀斌、 徐东、 白坤森、 刘毅、 郭子曦 + 1位作者

《Processes》：Prediction of Three-Dimensional In Situ Stress in Deep Coal Rocks Considering Heterogeneity and Physical Information Bing Li, Yunwei Kang, Pengcheng Hao, Weiping Zhu, Pengbo He, Huaibin Zhen, Dong Xu, Kunsen Bai, Yi Liu and Zixi Guo + 1 author

　　摘要 深层煤层气储层的特征是复杂的地质条件、强烈的非均质性以及原位应力的显著变化，这些因素给精确的三维原位应力预测带来了挑战。为了解决传统物理模型对岩石力学参数的强烈依赖性问题，以及纯数据驱动方法在样本量较少时缺乏物理约束和泛化能力的问题，本文提出了一种结合物理信息和数据聚类的LightGBM预测模型。研究对象选择了DJ区块中的1289个压裂簇。首先，使用K-means算法将储层划分为三类，以减少非均质性的影响。然后，为每一类构建一个LightGBM模型，并将基于黄模型的物理约束和应力-重力平衡关系纳入损失函数，以确保预测结果符合力学规律。以第一类压裂簇为例，该模型在测试集上的平均绝对百分比误差（MAPE）为2.78%，R2为0.89%。对比实验表明，该模型的预测精度优于BP神经网络、随机森林和Transformers。消融实验验证了聚类模块和物理信息约束的独立贡献和协同效应。转移实验表明，该模型适用于具有相似地质条件的区块。本研究提供了一种有效的方法，可以在保持准确性的同时，实现对深层煤层气储层原位应力的预测。

1. 引言
随着浅层煤层气和煤炭资源的枯竭，以及能源需求和开采强度的增加，深度在1500-3000米之间的深层煤层气逐渐引起了关注，煤层气开发已进入深层次阶段[1]。中国拥有丰富的深层煤层气资源，具有巨大的潜力，代表了非常规天然气勘探和开发的新前沿[2]。三维原位应力是压裂设计和操作中的关键参数，也是地下工程规划和建模的重要基础[3]。其准确表征直接影响到水平井轨迹设计、水力裂缝传播路径以及裂缝之间的相互作用，从而调节储层渗透率和流体流动通道[4]，最终影响单井的生产和回收率。准确预测原位应力可以为储层刺激方案提供科学依据[5]，并有效提高油气资源的开发效率。

深层煤层气储层的特征是埋藏深度深、结构复杂以及应力非均质性显著[6]，其力学行为受到水平主应力差、应力状态转变和煤层裂隙的共同控制[7]。准确确定三维原位应力的大小、方向和分布模式可以揭示其对煤层储层孔隙渗透特性的影响[8]，并指导水力压裂参数的优化，有助于形成有效的裂缝网络，从而突破低渗透率储层开发的瓶颈[9]。此外，这项研究还有助于评估井筒稳定性并优化生产技术[4]，支持深层煤层气开发中的风险控制和效率提升，促进具有挑战性资源的有效利用。深层煤层气储层的压裂刺激不仅涉及复杂的力学行为，还面临着温度-应力-损伤多物理过程的耦合挑战。李[10]通过循环热液氮冲击实验揭示了热力学耦合下的岩石力学退化和裂缝演化的机制；薛[11]进一步利用CT图像和灰度共生矩阵方法定量描述了冷冲击下的岩石微观损伤演化特征，发现损伤积累显著影响岩石的力学响应。这些研究表明，在深层煤层气压裂过程中，原位应力场并不是静态的，而是与储层损伤演化和温度扰动耦合的，这对原位应力预测模型的准确性和鲁棒性提出了更高的要求。

传统的原位应力获取方法主要分为两类：直接测量和经验模型。尽管直接测量方法如水力压裂[12,13,14]、岩心应力释放[15,16,17]、声发射（Kaiser效应）[18]和成像测井[19]相对成熟，但它们存在成本高、数据稀疏以及难以实现连续剖面预测的局限性。在理论模型方面，研究人员提出了一系列基于岩石力学理论的经验公式，包括Mattews & Kelly公式[20]、Anderson公式[21]、Newberry公式[22]和组合弹簧模型[23,24]。这些模型通常只考虑应力分布和岩石强度等主要因素，通过泊松比等参数估计原位应力。然而，它们没有充分考虑储层各向异性，并且过度依赖有限的岩石力学参数，导致实际应用中的准确性受到显著限制。随着从纯物理模型向数据驱动方法的转变，研究人员相继开发了基于数学统计和机器学习的原位应力预测方法。早期研究通过多元线性回归[25,26]建立了测井参数与原位应力之间的经验关系，或通过历史数据的曲线拟合获得原位应力[27]。近年来，各种机器学习方法被引入到原位应力预测领域，包括XGBoost、CatBoost和AdaBoost[28]；CNN-BiLSTM-Attention[29]；人工神经网络（MLR-ANN）[30]；双向长短期记忆网络（Bi-LSTM）[31]；粒子群优化（PSO）与BP神经网络（BPNN）结合使用[32]；极端学习机（ELM）和支持向量机（SVM）[32]；以及随机森林（RF）、ANFIS和功能网络（FN）[33]。此外，Abdiev[34]通过在吉尔吉斯斯坦一个构造复杂的矿区进行原位应力测量，建立了弹性波速度与应力参数之间的回归关系，清楚地揭示了构造断层对原位应力场的显著控制作用——这一发现为后续机器学习模型中的特征选择和地质信息嵌入提供了关键的先验知识支持。在煤矿结构薄弱区域的地球物理综合诊断中，Mussin[35]基于电阻率层析成像、折射地震和氡测量的多参数协同分析，建立了结构薄弱区域的定量识别系统，为原位应力场的间接推断提供了有效的地球物理支持。尽管数据驱动方法在预测精度和计算效率方面显示出显著优势，但它们仍面临关键挑战，包括高样本采集成本、缺乏物理约束以及神经网络弱泛化能力，这使得难以满足高度非均质深层储层原位应力预测的实际需求。特别是在数据稀少或外推情况下，纯数据驱动模型容易产生非物理预测（如负生产率或无法解释的应力值）。

为了克服纯数据驱动方法的“黑箱”局限性并赋予模型物理可解释性，研究人员引入了物理信息神经网络（PINNs）。PINNs将偏微分方程残差、初始条件和边界约束嵌入损失函数中，使神经网络在拟合数据的同时遵循已知物理定律。这一框架已成功应用于各种地质和岩石力学问题。在热力系统建模中，Ghasemi[36]使用自适应PINN解决了深层同轴钻孔热交换器的耦合热传递偏微分方程系统问题。该模型无需依赖观测数据（纯物理驱动），便实现了温度场的准确预测和可回收热能的准确估计，平均误差为1.24%，充分展示了PINNs在无数据条件下的问题解决能力。在生产预测中，Wang[37]进一步结合Caputo分数阶导数、自动微分和稀疏回归，从实际生产数据中识别出控制生产下降的非均匀分数阶微分方程。与纯数据驱动模型相比，这种PINN框架在外推预测方面表现出显著优势，证明了物理信息约束在提高模型外推性能中的关键作用。在岩石力学参数预测中，Yang[38]将基于声波传播时间和密度的动态弹性参数的经验公式作为物理约束嵌入神经网络，构建了一个物理-数据双重驱动的NSGA-PINN模型。该模型在弹性模量、泊松比、抗拉强度和断裂韧性方面的预测精度较高，显著优于纯数据驱动模型，并在不同区块中表现出优异的泛化能力。

为了解决这些关键问题，本文从深层煤层气田的压裂作业数据中获取三维原位应力样本，通过数据聚类和特征提取有效识别深层煤岩气储层的非均质性，并使用可解释的机器学习方法建立了三维原位应力预测模型。将基于经典岩石力学理论的物理约束引入模型，以确保预测结果的合理性和可靠性，从而实现三维原位应力的准确预测。在本研究中，三维原位应力标签是从校准的压裂作业数据中得出的，包括压裂闭合压力、破裂压力和从密度测井计算出的上覆岩层应力。

2. 材料与方法
2.1. 研究区域概述
DJ区块位于鄂尔多斯盆地东部边缘。该区块的结构剖面如图1所示。该区块的煤层气地质条件和储层特征表现出显著的深层煤层气积累特征。主要煤层由亮煤和半亮煤组成，属于高阶煤。在结构上，顶板和下底板主要由石灰岩和泥岩构成，提供了良好的密封能力和有利于深层游离气保存的条件。该区块的煤层气特征是吸附气和游离气的共存，游离气含量超过20%。煤层埋藏深度范围为2000至2520米，煤层厚度范围为2.0至9.8米，平均渗透率约为0.04米D。总体而言，该区块具有良好的储层特性和商业开发价值。

2.2. 数据来源与预处理
2.2.1. 最小水平主应力
在压裂泵停止后的闭合阶段，对数-对数图上的压力导数通常呈现斜率为3/2的直线段。当压力下降行为偏离这一理论斜率时，表明裂缝已完全闭合。拐点对应于闭合时间tc，此时点的压力下降ΔP(tc)可用于计算闭合压力Pc[39]：
(1) 在压裂分析中，闭合压力Pc数值上等于地层的最小水平主应力σh。因此，通过压力下降曲线分析可以获得最小水平主应力的数据样本。

2.2.2. 最大水平主应力
随着注入流体不断进入井筒，井底压力（BHP）持续上升。如果在恒定或增加的泵送速率条件下，BHP达到峰值后突然下降，表明某个穿孔簇发生了裂缝启动。这个峰值是该簇的破裂压力Pb。根据Hoek–Brown失效准则[34]，岩石的破裂压力与最小水平主应力正相关；因此，具有最低最小水平主应力的簇将首先启动。根据这一机制，可以从压裂作业曲线中准确匹配破裂压力。

在获得破裂压力Pb和最小水平主应力σh之后，可以使用改进的Hubbert–Willis公式[40]反计算最大水平主应力σH：
(2) 其中σf表示岩石的抗拉强度（MPa），Pp表示地层孔隙压力（MPa）。
通过上述操作获得的三维原位应力数据按压裂簇进行汇总。对于每个穿孔簇，提取相应的测井曲线数据，并将其与该簇的三维原位应力参数一一配对。这样就形成了一个结构化数据矩阵，压裂簇作为行，测井数据和三维原位应力作为列，作为后续聚类分析和三维原位应力预测模型训练的基础。

2.2.3.垂直应力
垂直应力主要由上覆岩层的重力引起，可以通过整合连续的密度测井数据来获得：
(3)
其中 σv 是垂直应力（MPa），ρ(z) 是深度 z 处的岩石密度（kg/m3），g 是重力加速度，D 是目标深度（m）。

2.2.4 数据预处理
共收集了 DJ 块中 1289 个具有完整数据的分裂簇作为研究对象，涵盖了 196 口深部煤层气井（包括 156 口垂直井和 40 口水平井）。每个分裂簇都与相应深度的测井数据相匹配。数据质量评估显示，井测井数据的整体缺失率约为 2.3%。
为了解决这些缺失值，采用了分段三次赫米特插值（PCHIP）方法，该方法利用测井曲线的垂直连续性和非线性趋势。与标准线性插值相比，这种方法更好地保持了数据的局部单调性，这对于维持地质信号的物理合理性至关重要。
对于异常值处理，应用了四分位数范围（IQR）方法，识别出大约 1.7% 的数据点为统计异常值。这些异常值被替换为定义深度窗口内相应测井参数的中值，从而减轻了它们对模型训练的不当影响，同时保持了数据分布的中心性。
最后，应用了最小-最大归一化方法，将七个测井输入特征的每个值独立缩放到 [0, 1] 范围内。这一步消除了测井参数之间单位和大小的不同影响，确保所有特征对后续的 K-均值聚类分析和三维原位应力预测模型的训练都有相同的贡献。

2.3 基于 K-均值算法的数据聚类
DJ 块在地质条件上表现出显著的横向和纵向变化， reservoir 的异质性很强，导致不同井段的分裂簇在物理性质、原位应力和岩石力学参数上有显著差异。为了全面描述这种异质性并提高后续原位应力预测模型的准确性，本研究使用单独的分裂簇作为基本研究单元，并应用 K-均值算法对深部煤层气分裂簇进行聚类分析，识别出具有相似地质和工程特征的分簇组。获得的聚类结果将用于构建分类和差异化的预测模型，从而有效提高关键参数（如三维原位应力）的预测准确性。
簇的数量 k 是 K-均值算法的核心参数，直接决定了聚类结果的合理性。本研究使用轮廓系数方法来确定最优簇数 k。该方法通过测量每个样本与其同一簇内其他样本的相似性（紧凑性）及其与最近相邻簇中样本的差异性（分离性）来评估聚类效果。轮廓系数的计算公式如下 [41]：
(4)
其中 si 是第 i 个样本（即第 i 个分裂簇）的轮廓系数；ai 是第 i 个样本与同一簇中所有其他样本的平均距离，反映了簇的紧凑性；bi 是第 i 个样本与最近相邻簇中所有样本的平均距离，反映了簇之间的分离性。整个数据集的平均轮廓系数是所有样本轮廓系数的平均值，范围从 ?1 到 1，数值越接近 1 表示聚类效果越好。
具体步骤如下：
(1) 选择最优簇数（k）对于有意义的分层至关重要。根据三个关键考虑因素，k 的候选范围限制在 2-5 之间：
(1) 可解释性和实用性：聚类的主要目标是识别出明确、可解释的地质力学单元（例如，由岩性驱动的组合），以支持分层建模。较大的 k（例如 >5）会导致统计上不稳定的碎片化簇，缺乏清晰的地质含义，从而降低维护多个子模型的操作简便性。
(2) 样本量限制：有 1289 个分裂簇，较大的 k 会大幅减少每个簇的平均样本数量，增加过拟合的风险，并降低每个簇对模型训练的代表性。
(3) 地质先验：对深部煤层气储层的初步探索性数据分析和区域地质知识表明存在 2-4 种主要岩石类型或应力状态。在评估了 k = 2 到 5 的平均轮廓系数后，确定 k = 3 为最优值，因为它产生了分离最好且最紧凑的簇。这三个类别对应于地质上可解释的组合：类别 I（特征是高伽马射线和低密度，典型于富含粘土/较软的层段），类别 II（中等相），以及类别 III（特征是低伽马射线和高密度，表明岩石清洁、致密、脆性）。这种三分法与储层的机械地层相一致，为后续的按类别预测建模提供了物理基础。
(2) 对于每个候选 k，执行 K-均值聚类。在聚类过程中，将初始化次数设置为 10 以减少随机初始质心的影响并提高聚类稳定性。
(3) 计算每个 k 的平均轮廓系数。
(4) 选择平均轮廓系数最大的 k 作为最优簇数。
为了表征控制原位应力分布的 reservoir 异质性，应用了 K-均值聚类算法。输入特征包括七个传统的井测井数据：自然伽马射线（GR）、自发势（SP）、声学旅行时间（AC）、补偿中子（CNL）、密度（DEN）、电阻率（RT）和温度（TEMP）。这种选择基于它们与原位应力关键决定因素的直接物理和经验关系。具体来说，AC 和 DEN 分别是动态弹性模量和上覆层应力的主要指标，这是控制水平应力大小的基本因素。GR 和 SP 作为岩性和粘土含量的代理指标，这些因素强烈影响岩石强度和应力敏感性。CNL 和 RT 提供了关于孔隙率和流体饱和度的信息，这与孔隙压力和有效应力状态有关。TEMP 与深度相关，反映了能够引起热应力的地热梯度。使用 K-均值算法进行了聚类分析，并使用轮廓系数定量评估了聚类效果。结果如图 2 所示。图 2. DJ 块中深部煤层气储层的聚类分析结果。(a) 不同簇数的平均轮廓系数。(b) 每个类别的分裂簇数量。如图 2a 所示，当 k = 3 时，轮廓系数达到最大值 0.72，表明聚类效果最好。因此，确定最优簇数为 3。图 2b 显示了每个类别的分裂簇数量：类别 I 包含 428 个簇，类别 II 包含 532 个簇，类别 III 包含 329 个簇，总计 1289 个有效样本。后续的原位应力预测将分别针对这三个类别进行，以减轻深部煤层气储层的异质性，并支持细化和差异化的开发策略。

2.4 基于物理信息的 LightGBM 模型
对于每个分裂簇类别，使用 LightGBM 构建独立的三维原位应力预测模型，以完全适应不同类型煤岩储层的机械响应特性。在模型训练过程中，将三维原位应力的物理约束纳入损失函数的构建中，以确保预测结果符合客观规律，从而提高模型的泛化能力。模型框架如图 3 所示。图 3. 基于 LightGBM 和物理约束的地质应力预测模型流程图。

2.4.1 LightGBM 模型
LightGBM（Light Gradient Boosting Machine）是一种基于梯度提升决策树的高效集成学习算法 [42]。它依次构建决策树以最小化损失函数。其核心特性包括基于直方图的决策树算法和具有深度约束的逐叶生长策略，这降低了计算复杂性，同时有效控制了过拟合，使其特别适合小样本场景。在本文的三维原位应力预测任务中，LightGBM 的输入是七个关键测井参数——GR、SP、AC、CNL、DEN、RT 和 TEMP——输出是三维原位应力。模型可以表示为：
(5)
其中 σtij 是第 i 个样本在第 t 次迭代时第 j 个输出指标的预测值，j ∈ {H, h, v}；η 是学习率；ft 是第 t 次决策树；Xi 是第 i 个样本的输入向量；i = 1, 2, …, N。
在每次迭代期间，LightGBM 旨在最小化训练损失。包含物理信息的损失函数将在后续部分构建。

2.4.2 三维原位应力的物理约束
深部煤岩储层表现出显著的异质性、层理结构和构造复杂性，导致原位应力场高度各向异性和空间变化。如果仅使用纯数据驱动的方法，模型容易过拟合或输出违反基本力学规律的非物理结果。因此，为了确保预测符合基本的地质力学规律，通过向损失函数添加额外的基于物理的惩罚项来引入物理约束。这些惩罚项源自黄模型和应力-重力平衡，可以对模型的预测进行微分。尽管 LightGBM 是决策树的集成，但提升算法的目标是最小化总损失函数，现在包括这些基于物理的惩罚项。因此，在训练过程中，梯度提升过程不仅受数据不匹配的指导，还受到偏离预设物理关系的指导，有效地引导集成朝向既准确又物理上合理的解决方案 [43]。
为了使预测的水平应力在物理上现实，我们结合了黄模型（也称为黄荣遵模型）[23]。这个经典的地质力学模型在中国深部煤层和页岩储层中得到广泛应用，因为它考虑了重力上覆层应力和构造效应。该模型表示水平应力的公式如下：
(6)
其中 σH、σh 和 σv 分别是最大水平应力、最小水平应力和垂直主应力。参数 ν 表示岩石的静态泊松比，α 是有效应力系数，Pp 表示地层孔隙压力。ωH 和 ωh 项是最大和最小水平应力方向的构造应力系数，定量表示区域构造力的贡献。物理上，该分量对应于在单轴应变条件下由重力压实产生的水平应力，而 和 表示由构造活动引起的额外水平应力。
水平应力的物理约束可以转换为损失项 Lh：
(7)
其中 σ*j,i 是第 i 个样本的水平主应力预测值，j ∈ {H, h}。
垂直应力主要由上覆层重量引起，必须满足源自连续介质力学的根本应力-重力平衡条件。忽略构造加速度并考虑一维垂直平衡，该条件表示为：
(8)
其中 ρ(z) 是深度 z 处的岩石密度，g 是重力加速度。为了确保模型的预测符合这一物理规律，我们制定了相应的损失项 Lv：
(9)
其中 σ*v,i 是第 i 个样本的垂直应力预测值，N 是用于模型训练的样本数量。
评估 Lv 的一个关键步骤是在 LightGBM 框架内计算深度导数，模型输出在此不可解析地微分。我们通过采用数值微分方案来解决这个问题。具体来说，对于深度为 zi 的样本，导数使用前向有限差分近似：
(10)
其中 Δz 是一个小的深度增量，设置为测井数据采样间隔，σ*v(zi) 和 σ*v(zi + Δz) 是当前 LightGBM 模型对深度 zi 和 zi + Δz 处特征向量的垂直应力预测。在梯度提升优化过程中，最小化这个损失项引导决策树集成产生满足物理规律离散形式的集体输出，从而将应力-重力平衡作为软约束嵌入学习过程中。

2.4.3### 结合物理信息的损失函数

用于训练三维原位应力模型的联合损失函数由几个部分组成。首先是三维原位应力的数据拟合损失项 \(Ldata\):
\[ \text{公式11} \]
其中 \( \sigma_H,i \) 是第 \(i\) 个样本的最大水平主应力的真实值，\( \sigma_h,i \) 是第 \(i\) 个样本的最小水平主应力的真实值，\( \sigma_v,i \) 是第 \(i\) 个样本的垂直应力的真实值。将水平主应力物理约束损失项 \(L_h \) 和垂直应力物理约束损失项 \(L_v \) 结合起来，模型训练的联合损失函数如下：
\[ \text{公式12} \]
其中 \( \lambda_1 \)、\( \lambda_2 \) 和 \( \lambda_3 \) 是用于平衡三个损失项贡献的权重系数。为了超越基于经验的做法和静态权重分配，并提高模型在不同地质条件下的鲁棒性，采用了自适应权重机制。该机制根据每个损失项的指数移动平均（EMA）动态调整权重，确保没有单个项在梯度下降过程中占据主导地位。具体来说，在训练迭代中，设 \( L(t)data \)、\( L(t)h \) 和 \( L(t)v \) 分别表示这三个损失项的当前值。它们的 EMA 分别表示为：
\[ \text{公式13} \]
其中 \( \beta \) 是一个平滑因子（设置为 0.9）。
当前迭代的自适应权重则计算为这些 EMA 的标准化倒数：
\[ \text{公式（此处省略具体计算公式）} \]
这里，\( \epsilon \) 是一个小常数，用于防止除以零。初始的 EMA 值使用小批量训练数据的损失值来设置。
因此，迭代中的总损失动态计算为：
\[ \text{公式14} \]
这种自适应方案自动平衡了不同损失成分的规模。如果某个损失项变得过大，其权重就会降低，防止它压倒其他项的优化信号。这导致更稳定的训练，并提高了模型泛化到具有不同特征的数据的能力。

### 模型评估指标

为了评估模型预测三维原位应力的准确性，引入了平均绝对百分比误差（MAPE）和决定系数 \(R^2\)。它们的表达式如下：
\[ \text{公式15} \]
\[ \text{公式16} \]
其中 \( n \) 是测试集中的样本数量；\( y_i \) 是第 \(i\) 个测试样本的评估指标的真实值；\( y_i^* \) 是第 \(i\) 个测试样本的评估指标的预测值；\( y_{\text{mean}} \) 是测试集上评估指标的平均值。

### 结果与讨论

将基于物理信息的 LightGBM 模型应用于 DJ 块。为了确保评估的公正性，遵循了标准的训练-测试分割协议。例如，在第 I 类的 428 个裂缝簇中，有 120 个样本被作为完全独立的测试集保留，而剩余的 308 个样本构成训练集。在这个训练集内，使用了 10 折交叉验证进行超参数调优，并提供了中间验证指标。调优后，最终模型在整个 308 个样本的训练集上进行了训练。本研究中报告的最终性能指标（例如在比较实验和消融实验中）是通过对独立的 120 个样本测试集进行评估得出的。在整个训练过程中，使用了自适应损失加权机制（第 2.4.3 节），从而无需手动预设固定权重。然后通过超参数优化、比较实验、消融实验和适用性分析全面评估了模型性能。

#### 3.1 LightGBM 模型的超参数优化

#### 3.1.1 每棵树的叶子数量

在 LightGBM 模型中，每棵树的叶子数量（\(Node\) 是一个关键的超参数，控制单棵决策树的复杂度。它的值直接影响模型的拟合能力和泛化性能。叶子数量过多容易导致过拟合，而叶子数量过少可能导致欠拟合。鉴于训练样本量为 308，\(Node\) 的候选值设置为 \(\{8, 16, 24\}\)。以最大水平主应力为例，十折交叉验证的可视化结果显示在图 4 中。如图 4 所示，当 \(Node\) 从 8 增加到 16 时，预测值逐渐接近真实值。进一步增加叶子数量会导致预测值偏离真实值。为了合理确定最优的 \(Node\) 值，计算了在不同 \(Node\) 值下预测三维原位应力的平均 MAPE 和 \(R^2\)，结果如表 1 所示。表 1 显示，当 \(Node\) 设置为 16 时，LightGBM 模型实现了最低的 MAPE 和最高的 \(R^2\)，表明预测准确性最高。因此，每棵树的最优叶子数量确定为 16。

#### 3.1.2 树的最大深度

树的最大深度（\(Depth\) 与叶子数量一起控制模型复杂度。深度过小会限制模型捕捉特征之间的高阶交互的能力，可能导致欠拟合；而深度过大则使单棵树过于复杂，增加过拟合的风险。\(Depth\) 的候选值设置为 \(\{3, 4, 5, 6\}\)。以最大水平主应力为例，十折交叉验证的可视化结果显示在图 5 中。如图 5 所示，当 \(Depth\) 从 3 增加到 5 时，预测值逐渐接近真实值，表明拟合能力有所提高。然而，当 \(Depth\) 达到 6 时，预测值开始偏离真实值，表明模型开始对训练数据中的噪声产生过拟合。为了合理确定最优的 \(Depth\) 值，计算了在不同 \(Depth\) 值下预测三维原位应力的平均 MAPE 和 \(R^2\)，结果如表 2 所示。表 2 显示，当 \(Depth\) 设置为 5 时，LightGBM 模型实现了最低的 MAPE 和最高的 \(R^2\)，表明预测准确性最高。因此，最优的树的最大深度确定为 5。

#### 3.1.3 每个叶子的最小样本数量

每个叶子的最小样本数量（\(Num\) 是另一个关键的正则化超参数，控制单棵树的复杂度。该参数强制每个叶节点至少包含指定数量的样本，从而限制树的生长深度和结构精细度。基于训练样本量，\(Num\) 的候选值设置为 \(\{5, 10, 15\}\)。以最大水平主应力为例，十折交叉验证的可视化结果显示在图 6 中。如图 6 所示，当 \(Num\) 从 5 增加到 10 时，预测值逐渐接近真实值，拟合效果有所改善。进一步增加到 15 时，一些预测值偏离真实值，表明过度的约束限制了模型的表达能力。为了合理确定最优的 \(Num\) 值，计算了在不同 \(Num\) 值下预测三维原位应力的平均 MAPE 和 \(R^2\)，结果如表 3 所示。表 3 显示，当 \(Num\) 设置为 10 时，LightGBM 模型实现了最低的 MAPE 和最高的 \(R^2\)，表明在拟合能力和泛化性能之间取得了最佳平衡。因此，每个叶子的最优最小样本数量确定为 10。

#### 3.1.4 学习率

学习率 \(\eta\) 是一个关键的超参数，控制每棵决策树的贡献。它的值直接影响模型的收敛速度和泛化性能。过大的 \(\eta\) 可能导致训练不稳定并错过最优解，而过小的 \(\eta\) 会导致收敛缓慢并增加迭代次数。考虑到训练样本量和优化后的叶子数量，\(\eta\) 的候选值设置为 \(\{0.01, 0.05, 0.1, 0.2\}\)。以最大水平主应力为例，十折交叉验证的可视化结果显示在图 7 中。如图 7 所示，当 \(\eta\) 从 0.01 增加到 0.05 时，预测值逐渐接近真实值，拟合效果得到改善。进一步将 \(\eta\) 增加到 0.1 时，预测值开始偏离真实值。当 \(\eta\) 增加到 0.2 时，预测值的离散度增加，拟合稳定性下降。为了合理确定最优的 \(\eta\) 值，计算了在不同 \(\eta\) 值下预测三维原位应力的平均 MAPE 和 \(R^2\)，结果如表 4 所示。表 4 显示，当 \(\eta\) 设置为 0.05 时，LightGBM 模型实现了最低的 MAPE 和最高的 \(R^2\)，表明预测准确性最高且收敛稳定。因此，最优的学习率确定为 0.05。至此，LightGBM 模型的关键超参数优化完成。所有后续分析均在此超参数组合下进行。

#### 3.2 与不同算法的比较

为了全面评估所提出模型的预测性能，选择了三种典型算法——BP 神经网络（BPNN）、随机森林（RF）和 Transformer——进行比较实验。为了确保比较的公平性，每种算法使用相同的数据聚类和基于物理信息的损失函数构建，并在相同的训练和测试集上进行训练和测试。训练周期数设置为 200，并使用测试集上的 MAPE 和 \(R^2\) 作为评估指标。模型在训练阶段的超参数设置如表 5 所示。为了确保严格和公平的比较，并明确展示我们方法的内在优势而非调整效果，Transformer 基线模型进行了广泛的优化以减少过拟合。这包括架构简化（减少到只有一个带有四个注意力头的编码器层，\(d_model = 32\)）、增强的正则化（应用dropout = 0.3 和 \(L2\) 权重衰减 = 0.001）以及精细的训练策略（使用 AdamW 优化器，学习率为 0.001，梯度裁剪和提前停止）。这些超参数是通过在验证集上的网格搜索确定的。因此，Transformer 的性能代表了在本研究小样本限制下的最佳可实现配置。即使经过这种全面优化，如表 6 所示，Transformer 模型仍记录了最高的 MAPE（15.21%）和最低的 \(R^2（0.47%），与基于树的模型相比明显表现较差。这一结果表明，标准的 Transformer 架构在数据稀缺的地球力学回归任务中容易过拟合，从而验证了选择 LightGBM 作为更高效样本基础的有效性，并突出了我们提出的基于物理信息和聚类信息的增强机制所带来的真正性能提升。以最大水平主应力为例，十折交叉验证的可视化结果显示在图 8 中，不同算法的 MAPE 和 \(R^2\) 如表 6 所示。表 5 显示了基于不同算法的模型超参数；表 6 显示了不同算法的实验结果比较；图 8 显示了最大水平主应力的实验结果比较。如图 8 和表 6 所示，所提出的基于物理信息的 LightGBM 模型在测试集上实现了最佳的预测准确性，其 MAPE 和 \(R^2\) 都优于其他比较算法。其中，作为基于树的集成模型的 RF 表现次优。BPNN 由于受到小样本条件的限制，泛化能力下降。然而，即使进行了上述详细的架构调整和正则化优化，Transformer 模型的预测性能仍然最差（MAPE：15.21%，\(R^2：0.47\)）。这表明标准的 Transformer 架构在本研究的小样本回归任务中存在根本性的挑战，因其本质上需要大量数据。这一结果反过来验证了我们选择LightGBM（一个更适合小样本场景的框架）作为基础，并通过物理约束和聚类先验来增强其有效性的策略。上述比较验证了所提出的模型在预测三维原位应力任务中的优越性。

3.3. 消融实验
为了系统地评估所提模型中每个模块的贡献，通过分别移除K-means聚类模块和物理信息约束模块来设计消融实验。为了确保公平比较并隔离每个模块的贡献，消融研究中使用了三种模型变体——基线LightGBM、LightGBM + K-means、LightGBM + Physics以及所提出的完整模型——它们都使用第3.1节确定的相同优化过的LightGBM超参数进行训练（即Node = 16、Depth = 5、Num = 10和学习率η = 0.05）。这保证了任何观察到的性能差异都可以归因于聚类和物理信息模块的包含或排除。比较的模型如下：(1) 基线模型（仅LightGBM）；(2) 仅添加了K-means聚类的模型（LightGBM + K-means）；(3) 仅添加了物理信息的模型（LightGBM + Physics）；以及(4) 完整模型（所提出的模型）。为了确保公平性，所有模型都采用了相同的训练/测试集划分、十折交叉验证和超参数策略。以最大水平主应力为例，图9展示了测试集上的可视化结果，直观反映了每个模块对预测稳定性的影响。表7总结了每个模型在测试集上的MAPE和R2指标，揭示了每个模块对模型性能的独立贡献和协同效应。

如图9和表7所示，基线模型（LightGBM）的泛化能力较差，MAPE为14.34%，R2仅为0.38%。仅添加K-means聚类模块后，预测值与真实值之间的整体一致性显著提高，MAPE降至9.13%，R2升至0.51%。然而，一些预测值仍然与真实值存在异常偏差，表明虽然聚类改善了局部拟合，但它无法约束预测结果以满足物理定律。仅添加物理信息模块后，预测值几乎从未与真实值出现异常偏差，但由于难以适应目标储层的异质性，整体预测准确性受到限制，MAPE为7.02%，R2为0.67%。结合了聚类和物理信息的所提出模型进一步将MAPE降低到2.78%，R2提高到0.89%，在所有变体中表现最佳。

上述结果表明，K-means聚类和物理信息约束模块具有明显的独立贡献和协同效应：聚类模块通过划分局部结构提高了模型对异质特征的适应能力，而物理信息模块确保了预测结果的物理合理性。它们的结合有效地提高了三维原位应力预测的准确性。为了进一步明确这两种物理机制之间的相互作用，我们进行了额外的消融实验，评估了不同的物理约束组合。测试了四种不同的配置：(1) 无物理约束的模型；(2) 仅具有黄氏模型水平应力约束的模型；(3) 仅具有应力-重力平衡垂直应力约束的模型；以及(4) 结合了两种约束的所提出模型。这些配置的性能指标总结如下：无约束模型的MAPE为7.81%，R2为0.63；仅黄氏模型的配置MAPE为8.25%，R2为0.60；仅应力-重力平衡配置的MAPE和R2也为8.25%和0.60；而完全的双约束模型实现了最佳性能（MAPE = 2.78%，R2 = 0.89）。这项比较分析表明，这两种约束既不冗余也不冲突，而是表现出强烈的协同效应：黄氏模型确保了预测水平应力的地质力学合理性，而应力-重力平衡则保证了垂直应力剖面的物理一致性。它们的联合应用为三维原位应力场提供了全面的物理基础，从而在预测准确性和稳定性上实现了显著且非加性的提升。这一结果明确了每种物理机制的具体贡献路径，并验证了在所提出的框架中集成应用的必要性。

此外，关于约束组合的消融研究显示，水平（黄氏模型）和垂直（应力-重力平衡）物理约束之间存在协同且非冗余的效应，这对于准确和物理一致地预测三维原位应力场至关重要。

3.4. 适用性分析
为了分析所提模型的适用性，我们在三个不同的煤岩气区块（SM-JX、LX-SF和HC）上进行了迁移应用实验。首先计算了每个区块中煤岩气井样本的欧几里得距离，以确定它们与DJ区块中各种样本类别的归属关系，然后调用相应的子模型来预测三维原位应力。使用MAPE和R2作为评估指标来定量评估模型的可迁移性。每个区块的储层特征和迁移应用结果如表8所示。表8显示，该模型在训练区块DJ上的表现最好。SM-JX区块的储层特征与DJ区块相当，模型的迁移性能良好。LX-SF区块的渗透率显著较高，导致迁移性能下降。HC区块由于埋藏深度较浅和气体含量较低，其储层特征与其他区块有很大差异，迁移性能最差，尽管仍然保持了可接受的预测能力。上述结果表明，所提模型对具有相似地质条件的区块具有良好的迁移适用性。对于储层特征差异较大的区块，虽然迁移性能有所下降，但模型仍具有一定程度的预测准确性。如果引入少量目标区块的样本进行微调，预计预测准确性会进一步提高。

为了进一步阐明HC区块性能下降的根本原因，对HC区块和DJ区块之间的关键特征分布进行了比较分析。如图8所示，HC区块的埋藏深度（300-1200米）明显小于DJ区块（2000-2520米）。因此，HC区块的密度（DEN）、声波时间（AC）和温度（TEMP）分布与DJ区块相比有显著差异。此外，HC区块具有较高的渗透率以及不同的孔隙压力和有效应力状态，这导致了地质应力主导特征的变化。这表明HC区块中的浅层煤层与DJ区块中的深层煤层在机械性质和发生条件上有所不同。训练集（DJ区块，深层）与目标域（HC区块，浅层）之间的特征分布不匹配导致模型的外推能力下降。为了解决这个问题，未来的工作可以探索域适应技术来增强跨区域的泛化能力。潜在的方法包括：(1) 使用最大平均差异（MMD）对源域和目标域的特征分布进行对齐；(2) 使用HC区块中少量标记样本进行迁移学习微调；(3) 域对抗神经网络（DANN）来提高模型在不同地质环境中的泛化能力。

3.5. 抗拉强度的不确定性分析
最大水平主应力（σH）的计算依赖于岩石的抗拉强度（σt），这是一个无法直接从常规井测井数据中获得的参数。它通常通过岩芯实验或经验公式进行估算，其固有误差为±5%至±20%，这种误差会通过Hubbert-Willis公式传递到σH的计算中。为了量化这种误差，将σt的相对误差设置为±5%、±10%和±20%，同时保持其他参数不变，并计算了σH的相对误差。结果如表9所示。表9显示，抗拉强度不确定性对最大水平主应力预测的影响通常是可控的。在σt的典型误差范围内，σH的相对误差小于3%，不会超过模型的预测误差。这表明该参数的不确定性不会显著影响模型的可靠性。

4. 结论
(1) 基于K-means算法，将DJ区块中的深厚煤岩储层分为三类，获得了0.72的轮廓系数。这有效地减少了储层异质性对原位应力预测的干扰，为后续的精细建模奠定了基础。
(2) 基于物理信息的LightGBM模型在测试集上的MAPE为2.78%，R2为0.89，预测准确性明显优于BP神经网络、随机森林和Transformer，验证了该方法在小样本条件下的优越性。
(3) 消融实验结果表明，K-means聚类和物理信息约束模块具有明显的协同效应：聚类模块提高了模型对异质特征的适应能力，而物理信息模块确保了预测结果的机械合理性。它们的组合将基线模型的MAPE从14.34%降低到了2.78%。
(4) 迁移应用实验展示了所提模型的条件可迁移性。对于与训练DJ区块具有相似地质特征的SM-JX区块，模型实现了高预测准确性（MAPE：5.33%，R2：0.83），显示出其在相似地质环境中的强大泛化能力。然而，对于埋藏深度和其他储层属性存在显著差异的HC区块，预测误差增加（MAPE：10.45%）。这一结果强调了模型性能依赖于目标域和训练域之间的地质相似性。尽管模型在差异较大的环境中仍保持基本的预测能力，但其最佳应用依赖于对地质相似性的预先评估或通过局部校准数据进行微调。