0 引言

纵横波速度比(v_P/v_S)是油气藏储层描述的关键参数之一，对于储层的含气性具有较高的灵敏度^[1]。要获取这一参数，主要方法包括岩石物理测试、横波测井和叠前反演。在实际油气田勘探中，岩石物理测试和横波测井往往成本高且仅有“一孔之见”，难以获得整个研究区的v_P/v_S参数。同时，叠前反演虽然可以获得高精度的v_P/v_S参数，但精度受地震资料品质影响较大，并且需要较多横波测井曲线用于标定^[2]。因此，开展高精度的v_P/v_S(或横波速度)预测对于油气储层的分析、评价至关重要，不仅具有重要的理论意义，而且也具有广阔的应用前景。

利用地震资料预测v_P/v_S，关键在于横波速度预测。目前，横波速度预测方法众多，有许多学者开展了深入研究。其中，基于人工智能的深度学习^[3]、机器学习^[4]等方法已经成为储层横波预测中的研究热点，并取得了一定的应用效果^[5-6]。然而，人工智能在模型性能及泛化能力、特征参数选择、参数优化、计算效率等方面仍存在诸多难点问题亟待解决。为了解决这些问题，本文引入XGBoost算法。该算法以梯度提升决策树(Gradient Boosting Decision Tree，GBDT)为基础，通过梯度提升策略，在每一轮迭代中优化残差，从而逐步提升XGBoost模型的性能。此外，该算法对异常值相对鲁棒，有助于提升模型的泛化能力。在特征选择与处理方面，XGBoost算法内置了特征重要性评估机制，能够自动选择对任务最相关的特征^[7]。这不仅简化了特征优选的过程，还有助于构建更简洁、解释性更强的模型。

前人已经在多个领域验证了XGBoost算法的有效性。例如，闫星宇等^[8]通过分析砂岩储层孔隙度与渗透率的关系，应用XGBoost算法构建了砂岩储层渗透率预测模型；谷宇峰等^[9]应用粒子群优化(PSO)算法对XGBoost算法进行改进，并用于测井渗透率预测，并与PSO优化支持向量回归、PSO优化GBDT等方法进行对比分析，认为PSO优化XGBoost的预测效率和精度最高，稳定性最好；丁阳阳等^[10]应用XGBoost算法进行多属性测井信息融合，识别煤体结构类型，相对于常用图版法，XGBoost能够实现煤体结构的高精度和高效率识别；张家臣等^[11]应用XGBoost算法在渤海湾盆地建立了测井曲线生成模型。此外，Al-Mudhafar^[12]和Merem-bayev等^[13]的研究表明，XGBoost算法在测井解释和岩相预测等方面具有优异表现。同样，XGBoost算法也应用到地震数据处理和解释，如在初至检测^[14]、岩性分类^[15]、砂体识别^[16]、断层检测^[17]等方面也取得了很好应用效果。

基于以上分析，本文提出一种基于XGBoost算法、应用地震属性预测储层v_P/v_S参数的方法。该方法应用独立成分分析(Independent Component Analysis，ICA)优选多种地震属性，并采用贝叶斯算法优化XGBoost模型的超参数。实际数据应用结果表明，该方法可为油气储层预测提供一种新思路。

1 方法原理

以从地震数据中提取的多种属性为输入，利用ICA对这些地震属性进行降维处理。然后，结合地震属性和测井信息，应用XGBoost算法对储层敏感的参数v_P/v_S进行预测。

1.1 地震多属性提取及优选

地震属性是通过数学方法从地震数据中提取的、能够反映地下物性及构造变化的数值化指标。这些属性数据体涵盖了地震波的几何学、运动学、动力学和统计学特征。地震多属性提取的目的在于从地震数据中挖掘与储层特征参数密切相关的信息，并转化为能够为地质解释和油藏描述所用的信息。自20世纪60年代出现地震属性概念后，现已提出200多种。通常采用独立成分分析等方法降低地震多属性的维度并优选属性，以获取与研究目标相关的地震属性，减少冗余信息，提高模型预测能力。

地震属性优选有很多方法，其中ICA是应用广泛的降维和特征提取方法。其基本思想是将混合信号视为多个独立成分的线性组合，通过统计分析和独立性检验分离出独立的成分，并选择最具有代表性的成分作为新的特征向量。这可以大幅降低数据维度，提高数据处理效率，同时保留主要特征，避免信息丢失。ICA可以发现数据中的线性关系和非线性关系，非常灵活和强大。具体实现方法包括最大熵法、最大似然法、快速ICA算法等^[18]。其处理流程包括准备数据集、中心化数据、白化数据、应用ICA算法、选择独立成分和反变换等，详细实现过程见文献[18]。在地震属性优化时，使用ICA可以消除冗余信息，提高数据利用率和模型预测能力，有助于发现数据中的隐含关系，从而提高预测结果的精度。

1.2 XGBoost算法

XGboost算法是GBDT算法的改进，二者都是基于梯度提升的集成学习方法。其核心思想是使用弱学习器(如决策树)进行迭代，不断优化集成模型的预测效果。GBDT仅使用一阶导数，而XGBoost则对损失函数进行二阶泰勒展开，同时引入一阶导数和二阶导数。这意味着任何满足二阶连续可导的函数都可以作为损失函数，从而提高了模型的稳定性和效率。XGBoost还采用了分裂节点时的贪婪算法和分位数近似技巧，加速了树的训练过程。此外，XGBoost对树的结构进行了限制，如最大决策树生成数量、单棵树最大分支深度、学习速率、样本随机抽取比例、构建树时向下采样率、分支树最小损失梯度、叶子节点最小样本数限制、子树权重的L₁和L₂正则化等超参数。通过调整这些超参数，可以优化XGBoost模型的性能，使其在训练数据上达到良好的拟合效果，同时在测试数据上保持强大的泛化能力。该算法还具有缺失值自动处理、稀疏数据处理及并行计算等功能，进一步提高了模型的效率和扩展性^[7]。XGBoost算法的流程如下(图 1)。

(1) 计算损失函数对样本的一阶导数和二阶导数。一阶导数为

$ {g}_{i}=\frac{\partial L({y}_{i}, {\widehat{y}}_{i})}{\partial {\widehat{y}}_{i}} $

(1)

式中：$ L({y}_{i}, {\widehat{y}}_{i}) $是损失函数；y_i是第i个样本的真值；$ {\widehat{y}}_{i} $是模型对第i个样本的预测值。二阶导数为

$ {g}_{i}^{\text{'}}=\frac{{\partial }^{2}L({y}_{i}, {\widehat{y}}_{i})}{\partial {\widehat{y}}_{i}^{2}} $

(2)

(2) 初始化模型

$ {f}_{0}\left({\boldsymbol{x}}_{i}\right)=\underset{\gamma }{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\sum\limits_{i=1}^{n}L\left({y}_{i}, \gamma \right) $

(3)

式中：f₀为初始模型(也称基分类器)；x_i为第i个样本的参数；γ为使损失函数达到最小的常数；n为样本个数。

(3) 迭代构建树模型。对于第t次迭代，首先计算当前模型f_t-1的预测值$ {{\widehat{y}}_{i}}^{(t-1)}={f}_{t-1}\left({\boldsymbol{x}}_{i}\right) $。其次计算样本残差$ {r}_{i}^{\left(t\right)}={y}_{i}-{\widehat{y}}_{i}^{(t-1)} $。第三，将上一步得到的残差$ {{r}_{\mathrm{i}}}^{\left(t\right)} $作为样本新的真实值，拟合新的回归树模型

$ {h}_{t}\left(\boldsymbol{x}\right)=\sum\limits_{j=1}^{{m}_{t}}{w}_{j}I(\boldsymbol{x}\in \boldsymbol{R}) $

(4)

式中：R为样本容量；w_j表示第t次迭代时的第j个叶子节点的分数；I为样本$ \boldsymbol{x} $对应叶子节点指示函数，即$ \boldsymbol{x}\in {\boldsymbol{R}}_{} $则I=1，否则I=0；m_t表示叶子节点数。第四，对于每个叶子节点，计算其分数w_j (j=1, 2, …, m_t)，w_j的计算公式为

$ {w}_{j}=-\frac{{G}_{j}}{{H}_{j}+\lambda } $

(5)

式中：G_j是第j个叶子节点上所有样本的梯度和；H_j是第j个叶子节点上所有样本的二阶导数和；λ是正则化参数，用于控制叶子节点权重的大小，防止过拟合。第五，更新模型

$ {f}_{t}\left(\boldsymbol{x}\right)={f}_{t-1}\left(\boldsymbol{x}\right)+\eta {h}_{t}\left(\boldsymbol{x}\right) $

(6)

式中η为学习率。

(4) 输出最终模型

$ f\left(\boldsymbol{x}\right)=\sum\limits_{t=1}^{T}\eta {h}_{t}\left(\boldsymbol{x}\right) $

(7)

式中T为最大迭代次数。

相对于GBDT算法，XGBoost算法具有许多优势。首先，它支持并行计算和优化，能够在较短的时间内处理大量数据，因此适用于大规模数据集和高维度数据。其次，该算法通过自动学习缺失值的分裂方向、为缺失值设定默认值、利用缺失值的统计信息和正则化等方式提高对缺失值和异常值的鲁棒性，能够处理不完整的数据集和含噪数据。此外，XGBoost算法还支持自定义损失函数和正则化项，可以根据不同的需求定制。无论是在处理分类问题还是回归问题，XGBoost算法都能通过灵活的自定义损失函数和正则化项，以及高效的优化方法，获得高精度的结果。

1.3 XGBoost超参数优化

XGBoost是一种强大的集成学习算法，但其性能受超参数选择的影响。在实践中，不恰当的超参数选择可能导致模型的性能不佳，进而影响预测的稳定性和精度。通过优化这些超参数，可以提高模型的准确性和泛化能力，防止模型过拟合或欠拟合，增强模型的鲁棒性，使其能够更好地适应地震数据集。

1.3.1 网格搜索

网格搜索是一种用于超参数调优的方法，其目的是确定模型的最佳超参数组合。在进行网格搜索时，首先定义一个超参数空间，该空间包含了所有可能的超参数组合。随后，利用交叉验证评估各种超参数组合的性能，并最终选择表现最佳的超参数组合作为模型的最终设置。

1.3.2 随机搜索

随机搜索是一种常见的超参数调优方法，与网格搜索相似，都是在超参数空间中寻找最佳的超参数组合。然而，与网格搜索不同的是，随机搜索通过在超参数空间中随机采样一组超参数组合来进行搜索，而不是遍历预定义网格中的所有可能组合。这种方法的优点在于，它可以在有限的时间内尝试更多的超参数组合，从而增加找到最佳组合的机会。此外，由于随机搜索可以在更大的超参数空间中进行探索，因此可以避免陷入局部最优解。

1.3.3 贝叶斯优化

贝叶斯优化是一种用于优化难以解析的黑盒函数的方法。其核心思想是通过代理模型来模拟目标函数，以选择下一个评估的样本点，从而最小化所需样本数量并在探索和利用间取得平衡。具体流程包括定义目标函数、超参数空间以及使用高斯过程作为代理模型。初始观测后，根据代理模型的预测进行迭代优化，选择下一个最佳超参数设置。当达到预定迭代次数或目标函数收敛时，优化停止并返回最佳超参数设置。相较于传统的网格搜索或随机搜索，贝叶斯优化能更高效地找到最佳超参数设置，并自适应地调整探索和利用的权衡^[19]，一般通过构建目标函数的概率模型和定义采集函数自适应地调整这种权衡。因此，本文选用贝叶斯优化算法对XGBoost的超参数进行优化。其流程如图 2所示。

1.4 评价指标

本文采用相关系数(r)、均方误差(Root Mean Squared Error，RMSE)、平均绝对百分比误差(Mean Absolute Percentage Error，MAPE)三种指标评价预测结果。r用来衡量两个变量之间相关性的指标，取值范围在-1到1之间。当r等于1时，表示两个变量完全正相关；当r等于-1时，表示两个变量完全负相关；当r等于0时，表示两个变量之间没有线性相关性。一般来说，相关系数的绝对值越接近1，说明两个变量之间的线性关系越强。r的定义为

$ r=\frac{\sum\limits_{i=1}^{n}({y}_{i}-\stackrel{-}{y})({\widehat{y}}_{i}-\overline{\widehat{y}})}{\sqrt{\sum\limits_{i=1}^{n}({y}_{i}{-\stackrel{-}{y})}^{2}\sum\limits_{i=1}^{n}({\widehat{y}}_{i}{-\overline{\widehat{y}})}^{2}}} $

(8)

式中：$ \stackrel{-}{y} $为真实值的均值；$ \overline{{\widehat{y}}_{}} $为预测值的均值。

RMSE是一种常用的衡量预测值和真实值之间差距的指标。通常情况下，RMSE越小，表示预测越准确。RMSE定义为

$ \mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{1}{n}\sum\limits_{i=1}^{n}({y}_{i}-{\widehat{y}}_{i}{)}^{2}} $

(9)

MAPE反映了预测值与真实值之间的相对误差。MAPE越小，表示预测精度越高。MAPE的定义为

$ \mathrm{M}\mathrm{A}\mathrm{P}\mathrm{E}=\frac{1}{n}\sum\limits_{i=1}^{n}\left|\frac{{y}_{i}-{\widehat{y}}_{i}}{{y}_{i}}\right|\times 100\mathrm{\%} $

(10)

2 Marmousi 2模型预测实验

Marmousi 2模型为完全弹性模型，横向延伸17 km，深度为3.5 km，包括199层^[20]。采用弹性波动方程有限差分法进行正演，炮间距为25 m，道间距为12.5 m，时间采样率为1 ms，经过系列处理，获得的纵波叠前时间偏移剖面。

为了验证XGBoost算法在预测横波速度的可行性，本文从Marmousi 2模型的叠后地震数据中提取振幅、频率、相位等地震属性，共20种(表 1)。考虑到地震属性对横波速度的敏感性不同，首先采用ICA对地震属性进行优化选择。通过ICA降维和优化处理，前四个分量的累计方差百分比为91.5%。在预测横波速度之前，首先建立地震属性与理论模型的横波速度识别标签。然后，采用XGBoost算法构建横波速度预测模型。随机选取80%的道作为训练集，20%的道为测试数据集。将降维处理后地震属性的前四个分量输入到XGBoost模型。由图 3的迭代误差曲线可知，基于XGBoost算法的预测模型在迭代到一定次数后，RMSE趋于稳定，说明模型在训练集和测试集上的性能已经达到了一个相对稳定的状态，得到了最优训练效果，此时进一步的训练不会显著提升模型的性能。

从预测的横波速度中抽取部分单道进行定量分析。由图 4的部分道横波预测结果可知：在单道横波曲线中，大部分与原始模型重合，但由于模型纵、横向变化较大，部分单道或单道的部分层段与模型存在较大差异。部分道预测结果与真实值的r、RMSE、MAPE统计如表 2所示，可以看出，相关系数最小值为0.88，说明预测值与实际值之间有较强的相关性；而RMSE和MAPE的最大值分别为6.55$ \times $10^-7、4%，说明预测值与实际值之间的误差小，预测结果可靠性高。但由于不同道对应的地层存在差异，在地层特征变化较大时，相应预测结果的准确性会下降，在第2581道误差较大，相关性也最低。

3 应用实例 3.1 工区概况

鄂尔多斯大牛地气田是中国境内最大的甲烷型天然气田之一。天然气广泛分布于浅层至深层，其中目标层马家沟组储层的平均厚度为50~80 m，最大厚度可达130 m。这一显著的储层厚度为气田的开发奠定了坚实的储量基础^[21]。而且，大牛地气田还有大量的未探明储量，因此具有巨大的勘探和开发潜力^[22]。

本文以鄂尔多斯盆地大牛地气田马家沟组碳酸盐岩储层为研究目标。收集了工区内马家沟组5 ¹⁺²亚段测井数据，包括纵波阻抗、横波阻抗、v_P/v_S等。图 5为v_P/v_S与纵波阻抗交会图，可以看出：含气白云岩和含水白云岩在v_P/v_S和纵波阻抗上也存在显著差异，因此可以选择这些v_P/v_S作为岩性识别、流体检测的参数。由于在研究工区内仅部分井有横波速度数据，无法获取全工区的横波速度数据，因此本文利用基于地震多属性的XGBoost算法直接预测v_P/v_S。

3.2 叠前反演

研究区马家沟组马5段为主要含气储层，其Xline1019测线叠前时间偏移剖面如图 6所示。为了验证XGBoost算法在鄂尔多斯盆地大牛地地区实际资料中应用的可靠性，本文与叠前反演获得的v_P/v_S数据进行对比、分析。

3.2.1 叠前道集优化

研究区的叠前道集经过叠前反褶积、道集拉平、振幅补偿等优化处理(图 7)。由图 7可知：研究工区叠前道集的目标层位的信噪比较高(绿线所示)，同相轴较连续，振幅变化较一致。对比图 7的G3气井和W3水井的井旁道集，AVO特征明显不同：气井井旁振幅随炮检距增大而减小；水井井旁振幅随炮检距增大而增大。储层AVO响应特征明显，能够区分含气储层和含水储层。

3.2.2 叠前同时反演

将叠前道集划分为0°~6°、6°~12°和12°~18°分别进行叠加，并从地震资料中提取地震子波，对工区内每口测井资料进行精细井震标定和追踪太原组、本溪组和马5⁴亚段的地震层位，结合工区已有纵波速度、横波速度和密度测井资料进行约束，建立初始的低频模型。然后采用约束稀疏脉冲算法，通过迭代提高模型的精度，并对低频模型进行修正，直到满足精度要求为止。

3.3 基于XGBoost算法预测v_P/v_S及应用 3.3.1 地震属性及超参数优化

本文共提取了实际地震数据的16种地震属性(表 3)，沿Tc2b向下移动6 ms处的部分属性切片如图 8所示。利用ICA算法分析属性之间的相关性(图 9)；随后提取出最具有代表性的成分，进行数据降维，去除冗余信息，提高数据的表达能力。ICA分析结果如图 10所示。

图 9是地震属性热图，其颜色表示相关大小，颜色越白或越黑，则表明两个地震属性的相关性越强。通过贝叶斯算法优化后，输入XGBoost模型的每个地震属性应该是最佳的，并能很好地映射到v_P/v_S值。图 10是对地震属性进行ICA降维后前7个分量的方差百分比统计，其中前4个分量的累计方差百分比达94.86%。因此，采用ICA前4个分量作为参数输入到XGBoost模型预测v_P/v_S。

本文将处理后的地震属性与前文计算得到的测井、地震反演的横波速度计算纵横波速度比，并组成数据集，该数据集将作为测试数据，用于XGBoost算法对纵横波速度比的预测。在具体实验中，首先采用贝叶斯优化XGBoost模型的最优参数组合(表 4), 再应用XGBoost模型预测大牛地马5段含气储层的v_P/v_S。

3.3.2 结果分析

研究区共有7口井，其中W1~W4井产水，G1~G3井试气结果分别为9.77×10⁴、1.58×10⁴、7.03×10⁴ m³/d。图 11为G1井的测井曲线，自上而下在2821.0 ~ 2824.0、2830.0 ~ 2833.0和2835.5~2838.5 m进行综合测试，无阻流量为9.77×10⁴ m³/d，获得了高产的工业气流。气层具有低自然伽马(GR)、中低孔隙度(5%~12%)、中低v_P/v_S(1.78~1.84)特征。

从叠前反演和基于XGBoost模型预测的v_P/v_S数据体中，提取两条典型的连井剖面。图 12a为过G1-G2-W3井的基于叠前反演获得的v_P/v_S，图 12b是基于XGBoost模型预测的v_P/v_S，其中图中的红色部分代表有利储层。显然，G1和G2井在马5¹⁺²亚段呈现明显的异常响应，而W3井含水层的v_P/v_S异常特征则相对不明显。图 13a揭示了过W1-G3井的基于叠前反演获得的v_P/v_S，图 13b是基于XGBoost模型预测的v_P/v_S。G3井在目标层也呈现出明显的v_P/v_S异常，但W1水井未出现明显异常。这一现象说明含气层与水层在v_P/v_S剖面上具有显著差异，从而证实了v_P/v_S能够较好的预测含气储层。两条过井剖面验证了应用v_P/v_S参数预测大牛地碳酸盐岩含气储层是可行的。对比图 12a与图 12b、图 13a与图 13b可见，两种方法获得的v_P/v_S剖面整体特征一致，说明本文方法能够达到叠前反演的效果。值得一提是，基于XGBoost模型的预测方法不需要叠前数据，只需要利用工区具有v_P/v_S测井的井旁地震道的地震属性来建立训练模型，然后推广应用到整个研究区，从而获得能够反应储层流体性质的v_P/v_S参数，为直接从叠后地震属性预测v_P/v_S提供了一种有效手段。

马5¹⁺²亚段沉积时期相对海平面持续下降，海水深度不断减小，加上干旱的气候，共同为白云岩的形成提供了有利条件^[23]，也是形成优质储层最有利的层段。图 14展示了沿Tc2b向下6 ms提取的预测v_P/v_S参数，其中红色代表强异常，对应有利储层区，W1、W2、W3、W4水井分布在强异常区外。表 5对比了工区盲测井的相同深度(即Tc2b向下移动6 ms处)的基于XGBoost算法预测的储层v_P/v_S参数与通过测井计算的v_P/v_S参数，可见，二者基本一致，表明基于XGBoost算法预测v_P/v_S是准确和可行的。

4 结论

本文提出了一种基于XGBoost算法的地震多属性v_P/v_S预测方法。在分析大牛地测井横波速度敏感参数的基础上，建立地震属性与储层横波速度之间的关系，采用XGBoost算法直接从叠后地震属性预测横波速度或v_P/v_S参数，获得的主要结论如下。

(1) XGBoost算法是一种基于决策树的集成学习算法，可以通过不断迭代提高模型的预测精度。采用贝叶斯算法对XGBoost算法的超参数进行优化，以增强模型的预测精度和泛化能力。

(2) 应用XGBoost算法对Marmousi 2模型进行横波速度预测，获得了精确的结果，从而证明了该方法的可靠性。

(3) 针对鄂尔多斯盆地大牛地气田的地震数据，利用XGBoost算法预测的v_P/v_S能够区分含气与含水储层，并有效地反映含气储层分布范围。与测井数据计算的v_P/v_S误差较小，表明XGBoost作为一种强大的机器学习算法精度较高，在地震数据解释中具有较大的应用潜力。