0 引言

对于致密砂岩储层，因受成岩作用和裂缝分布等多种地质因素影响，较难准确获取渗透率，因此在测井解释中渗透率预测方法始终为研究难点。渗透率预测模型主要分为物理模型^[1-3]和拟合模型^[4-6]。常用的物理模型有指数型^[1-3]、Kozeny-Carman^{[2-3, 7]}、Timur^{[2, 8-9]}和Herron^{[2, 8-9]}等。虽然上述物理模型基于测井理论，能得到可靠的渗透率预测值，但推广性较差，原因在于：难以准确获取模型的许多非测井曲线参数，如束缚水饱和度和岩石矿物含量等，导致在资料不充足的情况下无法应用；致密砂岩储层孔喉系统复杂，且易发育裂缝，使渗透率测井响应存在多解性，因此利用同一模型难以良好解释研究区所有层段。由于渗透率本质上是多种地质特征参数的综合反映，而这些参数大多与测井响应相关，因此从数学角度上讲可由测井曲线拟合渗透率^[4-6]。拟合模型一般仅用于计算测井曲线，且无需过多考虑地质因素的影响，极大减少了渗透率预测的计算量，因此在测井解释中得到快速发展^[10-12]。多元线性拟合和逐步迭代为经典的拟合方法。与多元线性拟合相比，逐步迭代考虑了自变量共线性问题，因此应用更广泛。刘建建等^[10]利用逐步回归建立了长6储层渗透率模型。刘敏^[11]以测井资料为基础，利用逐步回归估算了长庆T气田致密砂岩储层渗透率。由嘉雨^[12]利用逐步回归计算了榆树林油田葡萄花油层储层参数。虽然逐步迭代能快速预测渗透率，但也存在问题，如难以确定各类测井曲线与渗透率之间的复杂关系，因此无法用解析式表示目的层渗透率最优拟合模型。

为更智能地建立拟合模型，人们尝试利用机器学习技术预测渗透率。当前，拟合性能较好的机器学习模型包括支持向量回归(support vector regression，SVR)^[13-16]和梯度提升决策树(gradient boosting decision tree，GBDT)等^[17-20]。邓浩阳^[13]利用SVR估算高孔、低渗碳酸盐岩储层渗透率。李佳^[14]利用SVR和其他机器学习模型预测多孔介质渗透率，且SVR预测效果较好。Majid等^[15]、Zhang等^[16]利用以径向基函数为主的SVR预测储层渗透率，取得较好效果^[15-16]。Subasi等^[17]使用改进GBDT模型高效预测了常规油气储层渗透率。韩启迪等^[18]、谢云欣^[19]、Zhang等^[20]认为GBDT在测井解释方面具有发展前景。虽然上述拟合模型的预测效果较好，但仍存在一些应用问题：①由于建模时需多种经验参数(SVR需要惩罚系数等，GBDT需要回归决策树(classification and regression tree，CART)个数和学习速率等)，难以优化预测结果；②由于模型没有加入正则化项，在学习过程中易出现过拟合现象；③受自身算法的限制，模型在训练大量样本时非常低效。针对上述问题，人们提出了一种性能更好的模型——XGBoost^[21-23]，高效地实现了GBDT算法并进行了算法和工程上的许多改进。该模型加入了正则化项，并将拟合模型精细化，因此在理论上较SVR和GBDT等拥有更强的预测能力^[21-22]。闫星宇等^[23]利用XGBoost较好地预测了致密砂岩含气储层参数，但其采用试错法调试参数，效率低，且无法保证所用参数为最优组合。为此，本文从优化参数角度出发，利用粒子群优化(particle swarm optimization，PSO)算法改进XGBoost^[24-26]，进而提出渗透率预测模型PSO-XGBoost，并阐述该模型的方法原理及相关实验。

1 方法原理

因为XGBoost利用测井数据预测渗透率，所以训练样本和预测样本数据集应由测井曲线和实测渗透率数据构成，可表示为

$ \boldsymbol{A}=\left\{\boldsymbol{X}_{m n}, \boldsymbol{Y}_{m}\right\} $

(1)

式中：X_mn为测井数据集，m为样本量，n为测井曲线个数；Y_m为渗透率向量。在确定样本数据集后，预测模型可表示为^[21-23]

$ F\left(x_{i}\right)=f_{0}\left(\boldsymbol{x}_{i}\right)+\sum\limits_{k=1}^{K_{X}} \sum\limits_{j=1}^{T_{k}} w_{j, k} \cdot \eta $

(2)

式中：f₀(x_i)为初始模型(也称基分类器)，x_i为第i个样本；k=1, 2, …, K_X为迭代次数；j=1, 2, …, T_k为CART叶节点，T_k为第k次迭代的CART叶节点个数；w_{j, k}为第k次迭代的第j叶节点对应的所有样本的替代值；η为学习速率。

f₀(x_i)一般由损失函数确定，即

$ f_{0}\left(\boldsymbol{x}_{i}\right)=\underset{\alpha}{\operatorname{argmin}} \sum\limits_{i=1}^{m} L\left(y_{i}, \alpha\right) $

(3)

式中L(y_i, α)为损失函数，确定目标值与预测值之间的差值，y_i为第i个样本目标值，α为使损失函数达到最小的常数。w_{j, k}由

$ w_{j, k}=-\frac{\sum \partial L\left(y_{i}, \hat{y}_{i, k-1}\right) / \partial \hat{y}_{i, k-1}}{\sum \partial^{2} L\left(y_{i}, \hat{y}_{i, k-1}\right) / \partial \hat{y}_{i, k-1}^{2}+\lambda} $

(4)

确定。式中：${{{\hat y}_{i, k - 1}}} $为第i个样本在第k-1次迭代的预测值；λ为正则化系数。

XGBoost采用不同损失函数可得到不同的训练效果。考虑到平方损失函数对异常点或噪声点敏感，本文采用绝对损失函数建模^[21-22]。

由于XGBoost在建模中要设定许多经验参数，如CART叶节点个数、学习速率和正则化系数等，因此需优化才能保证建模质量。PSO能解决多目标最优化问题，因此文中用该技术优化XGBoost。如设定优化种群数量为q，则XGBoost参数种群可表示为

$ \mathit{\boldsymbol{ \boldsymbol{\varGamma} }} = \left\{ {{\mathit{\boldsymbol{\sigma }}_{i'}}\mid {\mathit{\boldsymbol{\sigma }}_{i'}} = \left( {{\sigma _{1i'}}, {\sigma _{2i'}}, \cdots , {\sigma _{zi'}}} \right)\quad i' = 1, 2, \cdots , q} \right\} $

(5)

式中$ {{\mathit{\boldsymbol{\sigma }}_{i'}}}$为种子向量，包含z个优化参数。

确定种群后，PSO的迭代公式可表示为^[24-26]

$ \left\{\begin{array}{l} \mathit{\boldsymbol{\varDelta}}_{i'}^{j'}=\left\{\boldsymbol{\sigma}_{i'}^{j'} \mid \boldsymbol{\sigma}_{i'}^{j'}=\min \left(\left\|\boldsymbol{O}_{i'} \boldsymbol{-} \boldsymbol{O}_{\mathrm{std}}\right\|_{2}\right) \quad j'=1, 2, \cdots, t\right\} \\ \mathit{\boldsymbol{\varDelta}}_{g}^{j'}=\left\{\boldsymbol{\sigma}_{\mathit{\Gamma}}^{j'} \mid \boldsymbol{\sigma}_{\mathit{\Gamma}}^{j'}=\min \left(\left\|\boldsymbol{O}_{\mathit{\Gamma}}^{\prime}-\boldsymbol{O}_{\mathrm{std}}\right\|_{2}\right)\right\} \\ \boldsymbol{W}_{i'}^{j'+1}=\omega \boldsymbol{W}_{i'}^{j'}+c_{1} r_{1}\left(\mathit{\boldsymbol{\varDelta}}_{i'}^{j'}-\boldsymbol{\sigma}_{i'}^{j'}\right)+c_{2} r_{2}\left(\mathit{\boldsymbol{\varDelta}}_{g}^{j'}-\boldsymbol{\sigma}_{i'}^{j'}\right) \\ \boldsymbol{\sigma}_{i'}^{j'+1}=\boldsymbol{\sigma}_{i'}^{j'}+\boldsymbol{W}_{i'}^{j'+1} \end{array}\right. \\ \text { s. t. }\left\{\begin{array}{l} \boldsymbol{\sigma}_{i'}^{j'} \in\left[\boldsymbol{\sigma}_{\max _{-} 1}, \boldsymbol{\sigma}_{\max _{-} 2}\right] \\ \boldsymbol{W}_{i'}^{j'} \in\left[-\boldsymbol{W}_{\max }, \boldsymbol{W}_{\max }\right] \end{array}\right. $

(6)

式中：$\mathit{\boldsymbol{\varDelta}}_{i'}^{j'} $为第i′个种子在前j′次迭代的最小计算误差；$ \mathit{\boldsymbol{\varDelta}}_{g}^{j'}$ 为在前j′次迭代的具有最小计算误差的种子；σ_i′^j′为在第j′次迭代的第i′个种子；σ_Γ^j′为在第j′次迭代的任意一个种子；O_i′^j′为第i′个种子在第j′次迭代的计算结果；O_Γ^j′为任意一个种子在第j′次迭代的计算结果；O_std为目标结果；t为迭代次数；W_i′^j′为第i′个种子在第j′次迭代的步长；ω为重力系数；c₁与c₂为加速系数；r₁和r₂为弹性系数；σ_{max_1}与σ_{max_2}为种子参数值界限；W_max为迭代步长最大界限。

为在PSO迭代前期锁定全局最优解方向，ω一般取较大值，而在迭代后期取较小值进行局部搜索。本文采用线性下降法(linear decreasing inertia weight，LDIW)优化ω^[25-26]。

综合上述，图 1总结了PSO-XGBoost计算流程^[21-26]。

2 实验验证 2.1 数据来源及实验设计

实验目的层为姬塬油田西部长4+5段致密砂岩储层。姬塬油田处于鄂尔多斯盆地中西部，位于天环坳陷和陕北斜坡两个一级构造单元之间(图 2a)。长4+5段是三叠系延长组主力层系之一，储层形成受控于沉积相带展布，多为水下分流河道砂体，部分为河口坝砂体^[27-30]。分析现有资料可知，储层孔隙以原生粒间孔(图 2c)、粒间溶孔(图 2d)和晶间孔(图 2e、图 2f)为主，在某些井区发育微裂缝(图 2g)，表明储层具有一定的储集空间和渗流通道^[27-28]。然而，在成岩作用阶段，受石英次生增大和自生黏土矿物等影响，孔喉结构以细孔—细微喉道和微孔—细微喉道组合为主，在压汞曲线上多体现为Ⅲ型和Ⅳ型(图 2h)，说明储层渗流性较差。因此，在评价储层时，能否获取准确渗透率数值将成为一关键研究点^[29-30]。

由岩心和测井等资料分析得知，目的层可划分为4个小层(长4+5₁¹、长4+5₁²、长4+5₂¹、长4+5₂²)。采用指数模型分析小层渗透率时，发现孔渗相关性较差(图 3)，其拟合优度R²一般小于0.5。对于其他常用物理模型，如Kozeny-Carman模型等，因基础资料不充足无法应用。以上情况说明常规物理模型并不适用于预测目的层渗透率，故文中采用PSO-XGBoost进行预测。为突出验证效果，在实验中对比了逐步迭代、SVR和GBDT三种模型的验证结果。验证井共有8口，分布在研究区东、西部(图 2b)，东、西部储层因受不同沉积和成岩作用影响而具有不同的品质。为综合验证PSO-XGBoost的预测能力，设计了三个实验。三个实验分别为：利用西区井资料(H1、H2、Y1)预测东区井资料(HA)；利用东区井资料(L1、H3、H4)预测西区井资料(HB)；利用全区井资料预测东、西区井资料(HA、HB)。前两个实验可了解不同储层品质的井区资料的相互预测效果，第三个实验考查在训练更多样本的情况下模型的预测能力。

2.2 实验1

实验测井数据由H1、H2和Y1井资料组成，共534个样本。每个样本由6条测井曲线构成，包括声波时差(AC)、补偿中子测井(CNL)、补偿密度测井(DEN)、自然伽马(GR)、自然电位(SP)和原状地层电阻率(RT)。样本对应的渗透率值由气测实验得到，变化范围为0.01~7.60mD。表 1展示了部分训练样本数据。由于对比模型SVR和GBDT在计算中也涉及经验参数，因此为使验证对象都处于相同优化状态，也用PSO进行优化。对于逐步迭代，由于算法只利用自变量和因变量数据进行分析，不涉及处理其他经验参数，因此无需优化。由PSO-XGBoost计算流程(图 1)可知，因增加PSO优化流程，训练样本分为两部分，一部分用于建模，另一部分用于PSO迭代检测。为此，在实验中将训练样本以比例0.7/0.3分配给建模和优化检测，即374个样本用于建模，160个样本用于优化检测。

虽然逐步迭代公式无需优化，但为防止所建模型过拟合，也采用上述训练样本分配方案进行建模和检测。应用逐步迭代之前，先要确定测井曲线与渗透率之间的关系式。由于该模型只是从数学角度出发对数据进行处理，且为加强通用性，一般不加入交叉项和幂指数项，所以设定的渗透率与所用的测井曲线的关系式一般为线性^[8-10]。渗透率和RT常呈数量级变化，但目的层的RT变化不大，因此建模时将渗透率数据K以对数形式分析。由374个训练样本拟合得到的逐步迭代模型为

$ \begin{array}{r} \lg K=0.017 \mathrm{AC}+0.028 \mathrm{CNL}-1.91 \mathrm{DEN}- \\ 0.021 \mathrm{GR}+0.013 \mathrm{SP}+0.003 \mathrm{RT}+1.746 \end{array} $

(7)

式(7)的调整拟合优度R²(0.68)高于常规指数模型(图 3)。因此，与常规指数模型相比，以多种测井曲线资料为基础，利用逐步迭代建立的拟合模型能更好地反映致密砂岩储层渗透率的变化趋势。采用MSE (mean squared error)评价计算误差^[7-12]，因此误差单位应为mD²。在下文中，为表达简便，MSE只讨论其数值。统计逐步迭代模型的160个检测样本的验证结果(图 4a)表明，72.5%的结果与目标值处于同一数量级，说明该模型在预测同一井区的数据时效果较好。

预测样本为HA井数据，共50个，其构成比例与训练样本一致。利用逐步迭代模型进行预测的MSE较大(0.8744，表 2)，说明逐步迭代建立的拟合模型在不同储层品质的井区资料的相互预测中效果不理想。

对于SVR，首先确定参数框架，之后再进行PSO优化建模，其中3个参数需要优化^{[13-16, 24-26]}。需要说明的是，重力系数和加速系数虽为经验值，但一般不需要改动，即大幅度改动这些参数也不明显影响优化结果，侧面说明利用PSO技术进行优化可减少模型调参工作量。SVR参数优化结果与初始设定值差别较大，说明参数优化的必要性(表 3)。优化过程(图 5a)表明，在100次迭代中，MSE由1.400降至0.4919，说明了PSO的优化作用。PSO-SVR的验证结果统计(图 4b)表明，84.38%的计算结果与目标值处于同一数量级，其可靠程度较逐步迭代模型有所提升，表明PSO-SVR能更好地处理来自同一井区的数据。利用PSO-SVR的MSE较逐步迭代小(0.7212，表 2)，但仍不理想，说明该模型在利用其他井区资料预测渗透率时效果不好。

对于GBDT，表 3分别展示了初始参数设定、优化参数设定和参数优化结果^[17-20]。优化过程(图 5a)表明，MSE由0.8000降至0.1124，表明PSO在参数优化方面的显著作用。PSO-GBDT验证结果统计(图 4c)表明，92.5%的结果与目标值数处于同一数量级，说明在处理检测样本时该模型的预测能力较逐步迭代、SVR强，且MSE(0.1526，表 2)远小于上述两种模型，说明该模型能有效利用其他井区资料预测目标井数据。

对于XGBoost，表 3分别展示了初始参数、PSO优化参数和参数优化结果^[21-23]。优化过程(图 5a)表明，MSE由0.7000降至0.1059，充分说明XGBoost的参数优化得益于PSO技术。PSO-XGBoost验证结果统计(图 4d)表明，94.37%的结果与目标值数处于同一数量级，表明在处理来自同一井区的资料时，该模型的预测能力强于前三种模型，且MSE(0.1273，表 2)最小，证明该模型能有效利用其他井区资料高精度地完成预测任务。

图 6为实验1预测结果柱状图。由图可知，PSO-GBDT和PSO-XGBoost落在计算误差列中的数据点最多，但后者的数据点整体更偏向0线，表明PSO-XGBoost预测结果的精度最高。

2.3 实验2

本实验训练数据由L1、H3和H4井资料构成，共505个样本，样本构成要素同实验1。渗透率值由气测实验得到。在实验中将训练样本以比例0.7/0.3分配给建模和优化检测。

对于逐步迭代，建模得到的拟合公式为

$ \begin{aligned} \mathrm{lg} K=&-0.017 \mathrm{AC}-12.36 \mathrm{DEN}-\\ & 0.09 \mathrm{GR}-0.023 \mathrm{SP}+0.018 \mathrm{RT}+36.803 \end{aligned} $

(8)

式中的拟合变量缺少CNL，说明CNL曲线在实验中对于渗透率变化无影响。对比式(7)和式(8)发现，逐步迭代在处理不同井区资料时得到的拟合表达式差别较大，表明模型的推广性差，即由本井区资料所建的拟合模型不适用于预测其他井区资料。预测数据由HB井资料构成，共50个样本，其构成要素与训练样本一致。利用逐步迭代模型进行预测的MSE(0.3479，表 2)明显小于实验1，表明该模型的预测能力受建模数据影响较大，因此模型的预测稳定性非常差。

对于SVR，参数及其优化设定同实验1。优化过程(图 5b)表明，MSE由约1.2000降至0.3871，表明模型的预测能力受PSO优化而得到提升。PSO-SVR的MSE(0.1756，表 2)明显小于实验1，表明该模型的预测能力受建模资料品质影响较大，因此计算稳定性不佳而难以推广应用。

对于GBDT，参数及其优化设定同实验1。优化过程(图 5b)表明，MSE由约0.8000降至0.1106，说明PSO起到了良好的优化作用。PSO-GBDT的MSE(0.1413，表 2)与实验1相近，表明该模型能利用其他井区资料有效预测渗透率，而且预测稳定性较好，能够推广应用。

对于XGBoost，参数及其优化设定同实验1。优化过程(图 5b)表明，MSE由约0.7000降至0.0988，说明PSO技术在建模中起到了良好的优化作用。PSO-XGBoost的MSE(0.0667，表 2)最小，再次证明PSO-XGBoost能够利用其他井区资料有效预测渗透率，而且预测能力最强。

图 7为实验2预测结果柱状图。由图可见：虽然PSO-SVR的计算误差与PSO-XGBoost的相近，且都靠近0点，但在计算误差列的顶端和末段都处于0.5~1范围内，且中间有些点的误差已经大于1；PSO-XGBoost在这些区域的点都基本落在0~0.5范围内，表明PSO-XGBoost的计算结果更精确。

2.4 实验3

在本实验中，将由东、西井区资料合并构成的训练数据进行预测。训练样本有1039个，对于所有模型，将训练样本以比例0.7/0.3分配给建模和优化检测。

对于逐步迭代，得到的拟合公式为

$ \begin{aligned} \mathrm{lg} K=& 0.032 \mathrm{CNL}-7.06 \mathrm{DEN}-0.014 \mathrm{GR}-\\ & 0.012 \mathrm{SP}-0.004 \mathrm{RT}+19.198 \end{aligned} $

(9)

模型经算法判断将AC曲线去除，形成的拟合表达式与实验1、实验2的不一致，再次说明在处理不同的建模数据时逐步迭代容易形成不同的拟合公式，难以推广应用。对于目标井的MSE(分别为0.3360和0.1581，表 2)较实验1、实验2的小，表明在训练更多学习样本的情况下，模型的预测能力得到明显提升。

对于SVR，优化过程(图 5c)表明PSO具优化作用，PSO-SVR的MSE(分别为0.1534和0.1146，表 2)小于实验1、实验2，证明该模型的预测能力在处理更多训练样本的情况下得到加强。同时，MSE与实验1的相差较大，与实验2的相差较小，也再次表明模型的预测能力受训练数据品质影响较大，反映了模型的预测稳定性差。

对于GBDT，PSO-GBDT的MSE(分别为0.1352和0.0990，表 2)小于实验1、实验2，但相差不大，说明模型不仅在训练更多学习样本后形成更强的预测能力，而且模型的预测稳定性较好，能够推广应用。

对于XGBoost，PSO-XGBoost的MSE(分别为0.1259和0.0378，表 2)小于实验1、实验2，且降幅较小，说明XGBoost的预测能力不仅在处理更多训练样本后得到提升，而且预测稳定性好。与PSO-GBDT相比，XGBoost的MSE更小，说明该模型更具备推广价值。

由图 6和图 7的计算误差列可看出，在处理更多建模样本的情况下，各验证模型的MSE较实验1、实验2有所减小，很好地说明了扩大训练样本容量是一种提升各验证模型预测能力的有效手段。同时，PSO-XGBoost的MSE最小，再一次说明该模型预测渗透率的能力最强，最具推广应用价值。

2.5 计算时间分析

表 2展示了各验证模型在不同实验中的计算耗时。从整体上看，PSO-SVR耗时最少，PSO-XGBoost耗时最多。PSO-XGBoost的计算速度最慢是因为在实验中没有采用XGBoost并行计算技术。考虑到在处理大数据时，XGBoost计算耗时会随着回归树量的增加而成倍增长，为此设计了能够并行计算的数据存储方式及其处理方法^[21-23]。在处理大于10万个样本或样本包含上百个变量的数据体时，才能明显体现XGBoost的并行计算效率。虽然PSO-XGBoost在处理所用的训练样本(小数据体)时耗时最长，但与其他模型的耗时差距并不大，在可接受范围之内，并且所得结果的精度更高，说明PSO-XGBoost在处理一般小型数据体时效率较高。

3 结论

本文从优化参数角度出发，利用PSO算法改进XGBoost^[24-26]，得到渗透率预测模型PSO-XGBoost，阐述了该模型的方法原理并进行了预测致密砂岩储层渗透率实验，得到以下认识：

(1) 与物理模型相比，拟合模型涉及的储层特征参数较少，因此在建模资料不充足的情况下，适用性更好。

(2) 致密砂岩储层因孔喉系统复杂且易发育裂缝，使渗透率在测井响应上与孔隙度相关性差，导致指数模型适用性差。

(3) 与指数模型相比，逐步迭代能更好地反映渗透率与测井曲线之间的线性响应关系，但由于预测性能受建模数据品质影响很大，导致计算稳定性不佳，使其难以推广使用。

(4) SVR参数能够由PSO技术优化，且形成的PSO-SVR在处理小型数据体时能快速给出预测结果。但模型的预测性能随训练样本品质的变化而有较大的波动，导致模型预测稳定性较差，因此PSO-SVR难以推广使用。

(5) 使用PSO技术可优化GBDT参数，形成的PSO-GBDT不仅能快速给出可靠的预测结果，且模型的预测性能随训练样本品质变化的波动不大，表明具备推广应用性。

(6) 使用PSO技术可优化XGBoost参数，且形成的PSO-XGBoost预测效率高，稳定性好，具备推广使用性。与PSO-GBDT相比，PSO-XGBoost的预测精度更高，表明更具推广使用价值。

(7) 逐步迭代、PSO-SVR、PSO-GBDT、PSO-XGBoost的预测能力都可在训练更多学习样本后得到明显提升。