0 引言

近年来，油气勘探的目标储层日趋复杂，勘探难度逐步增加，复杂(岩性)油气藏渐渐成为研究与开发的重点。煤系地层在中国分布范围广阔^[1]，包括沁水盆地、鄂尔多斯盆地、四川盆地、东海盆地等，是油气勘探的重要领域^[2-4]。煤系烃源岩生烃潜力大。根据源岩有机碳含量从高到低顺序将其划分为煤岩、碳质泥岩和暗色泥岩等类型^[5-6]。煤成气理论拓展了中国油气勘探的领域，且促进了天然气工业的大发展。根据地质资料可知，煤系地层中的煤层一般具有层数多、厚度薄、横向分布广、垂向分布散等特点^[7]，为储层预测带来了极大挑战。煤层相较于常规的砂、泥岩，在测井响应上具有低密度、高声波时差和高中子值等特征^[8]。目前针对煤系地层岩石物理性质的研究相对较少，无法有效地阐述地层弹性参数与物性参数之间的关系。

地层物性参数是描述地下介质的重要信息，通过地震岩石物理建模可有效地将物性信息过渡到岩石的弹性特征，使得地层物理性质在地震响应中有所体现。岩石的性质主要受矿物组分、体积占比、排列方式、孔隙度、孔隙几何形状、流体分布及周围环境等因素影响^[9]。当地震波长足够大时，地下岩石可视为均匀物体，就能利用等效介质理论和孔隙弹性理论构建地震岩石物理模型，量化物性参数与弹性参数之间的关系。在岩石物理建模方面，人们针对矿物岩性、孔隙结构、流体性质等做了大量研究^[10-13]，分析了物性参数对弹性性质的影响。经典的Xu-White模型将石英和黏土作为碎屑岩的背景矿物，并综合考虑了孔隙和流体的影响^[14]，得到了广泛应用。在Xu-White模型基础上，改变背景矿物并将孔隙类型替换为裂缝、粒间孔和溶洞，可构建适用于碳酸盐岩地层的岩石物理模型^[15]。随着页岩气的大力开发，众多学者聚焦于页岩的物理性质研究^[16-18]，充分考虑干酪根等有机物的形态与孔隙特点，计算岩石的等效模量。在实际应用中，常缺少横波测井资料，利用岩石物理模型计算横波速度等弹性参数已经是一种常用手段^[19-21]。此外，从岩石物理角度，建立物性参数与地震响应之间的联系，可为储层预测奠定基础^[22-25]。因此，针对煤系地层的岩石物理模型构建及其性质研究对油气储层的有效开发具有重要意义。

本文充分考虑地层含煤对地震弹性参数的影响，基于VRH平均模型^[26]、K-T骨架模型^[27-28]、SCA模型^[29-30]和Gassmann方程^[31]，构建了一种煤系地层的地震岩石物理模型。在煤系地层等效介质基础上，进一步计算多种岩石弹性参数，分析多种微观物性参数变化对不同弹性参数的影响，得到区分不同岩性的敏感参数。最后，在测井资料反演约束下构建目标泛函，形成一种煤系地层岩石弹性参数非线性计算方法。利用该方法求得的横波速度的(预测)误差较小且呈正态分布，在煤层处与实测横波速度吻合度较高，表明煤系地层地震岩石物理模型构建的有效性。

1 岩石物理建模

根据煤系地层地质特点和相关背景，构建了考虑煤层影响的油气储层地震岩石物理模型，具体过程分为四个部分，依次为背景岩石基质构建、孔隙填充、煤层耦合、流体替换(图 1)。

1.1 基质构建

首先利用VRH平均模型将石英、黏土等背景矿物按照各成分体积含量进行混合，VRH模型^[26]是通过求取等应力平均与等应变平均的均值来表征矿物基质的等效岩石模量，具体表达式为

$ K_{\mathrm{m}} =\frac{\sum\limits_{i=1}^{N} f_{i} K_{i}+\left(\sum\limits_{i=1}^{N} \frac{f_{i}}{K_{i}}\right)^{-1}}{2} $

(1)

$ \mu_{\mathrm{m}} =\frac{\sum\limits_{i=1}^{N} f_{i} \mu_{i}+\left(\sum\limits_{i=1}^{N} \frac{f_{i}}{\mu_{i}}\right)^{-1}}{2} $

(2)

式中：K_m、μ_m分别为混合后的矿物基质体积模量和剪切模量；f_i、K_i、μ_i分别为第i种矿物组分的体积含量、体积模量以及剪切模量；N为不同矿物类型的数目。

1.2 孔隙填充

K-T模型^[27-28]基于散射理论，可得到矿物基质包含不同形状孔隙时的骨架等效模量

$ \left(K_{\mathrm{KT}}-K_{\mathrm{m}}\right) \frac{K_{\mathrm{m}}+\frac{4 \mu_{\mathrm{m}}}{3}}{K_{\mathrm{KT}}+\frac{4 \mu_{\mathrm{m}}}{3}}=\sum\limits_{i=1}^{N} x_{i}\left(K_{i}-K_{\mathrm{m}}\right) P_{\mathrm{m} i} $

(3)

$ \left(\mu_{\mathrm{KT}}-\mu_{\mathrm{m}}\right) \frac{\mu_{\mathrm{m}}+\zeta_{\mathrm{m}}}{\mu_{\mathrm{KT}}+\zeta_{\mathrm{m}}}=\sum\limits_{i=1}^{N} x_{i}\left(\mu_{i}-\mu_{\mathrm{m}}\right) Q_{\mathrm{m} i} $

(4)

$ \zeta_{\mathrm{m}}=\frac{\mu_{\mathrm{m}}}{6} \frac{8 \mu_{\mathrm{m}}+9 K_{\mathrm{m}}}{2 \mu_{\mathrm{m}}+K_{\mathrm{m}}} $

(5)

式中：K_KT、μ_KT分别为加入孔隙后的岩石骨架等效体积模量和剪切模量；x_i为不同孔隙所占的体积分量，若为干燥孔隙，则K_i＝μ_i＝0，若为流体饱和，则μ_i＝0；P_mi、Q_mi反映的是不同形状的孔隙对矿物基质的作用。

K-T模型仅适用于较低孔隙度的岩石模量估算，常与微分等效介质模型结合使用。Keys等^[32]对此模型进行改进，假定骨架的泊松比为常数，推导出岩石骨架等效体积模量和剪切模量表达式

$ K(\phi)=K_{\mathrm{m}}(1-\phi)^{p} $

(6)

$ \mu(\phi)=\mu_{\mathrm{m}}(1-\phi)^{q} $

(7)

式中：K(ϕ)、μ(ϕ)分别代表加入孔隙后的骨架体积模量和剪切模量；ϕ表示孔隙度；p、q为几何因子，与背景介质和孔隙形状相关。

1.3 煤层耦合

SCA模型^[29-30]是计算多相混合物介质模量的有效模型，向骨架中以包含物形式加入煤层(图 2)，通过调整背景介质参数计算耦合煤后的岩石剪切模量和体积模量，具体表达如下

$ \sum\limits_{i=1}^{N} f_{i}\left(K_{i}-K_{\mathrm{d}}\right) P_{i}^{*}=0 $

(8)

$ \sum\limits_{i=1}^{N} f_{i}\left(\mu_{i}-\mu_{\mathrm{d}}\right) Q_{i}^{*}=0 $

(9)

式中：f_i在此处特指背景骨架或煤的体积含量；K_d、μ_d是加入煤层后的干岩石模量；P_i^*和Q_i^*是与(包含)煤的形状、性质以及骨架背景有关的几何因数，当煤的形状不同时，P_i^*和Q_i^*的表达式也会相应改变。上述两式均为耦合的，需利用迭代方式求解。

在建模过程中，煤层的耦合使用的是SCA模型，煤层占比及孔隙度均较低。假设在理想情况下，煤层以包含物的形式添加时，不会替换掉孔隙，此时得到的背景介质模量依旧为干岩石骨架模量。分析煤层和孔隙度的不同添加顺序对岩石骨架弹性模量的影响(图 3)，可见在不同煤层占比情况下，先加孔隙后加煤层(图 3a)与先加煤层后加孔隙(图 3b)得到的体积模量和剪切模量较接近，说明煤层和孔隙的建模顺序对数值计算结果影响较小。但岩石物理建模只是对地下岩石构成的一种等效近似，并不能表征实际岩石形成过程，煤层和孔隙的具体建模顺序可依据目标工区实际情况做调整。

1.4 流体填充

Gassmann^[31]将地下介质简化为均匀各向同性，在一定假设前提下，建立了饱和岩石等效模量与流体类型、矿物基质、岩石骨架之间的关系，推导出含煤层骨架填充流体后的体积模量和剪切模量

$ K_{\mathrm{s}}=K_{\mathrm{d}}+\alpha^{2} K_{\mathrm{f}} /\left[\phi+(\alpha-\phi) K_{\mathrm{f}} / K_{\mathrm{m}}\right] $

(10)

$ \mu_{\mathrm{s}}=\mu_{\mathrm{d}} $

(11)

式中：K_s、μ_s分别为流体替换后的岩石等效体积模量和剪切模量；α=1-K_d/K_m；K_f为流体等效体积模量，这里可利用Reuss平均求取。

针对上面构建的煤系地层地震岩石物理模型，推导得到岩石等效体积模量和剪切模量，再进一步计算多种岩石弹性参数。本文利用实际测井数据验证模型的合理性。

图 4中黑线是根据所建模型得到的纵、横波速度范围，彩色圆点为实测纵、横波速度样点。样点泥质含量高其速度偏低；泥质含量低则速度较高。可见二者吻合性较好，证明了模型构建的可靠性。

2 主控因素分析

岩石物理方法在物性参数与弹性参数之间起到了桥梁的作用，分析岩石弹性参数对物性参数变化的响应十分重要。因此，针对煤系地层特点，本文分析了煤层占比、泥质含量、孔隙度、含水饱和度等物性参数的渐变对不同弹性参数的影响。在所建模型中，主要矿物为石英、黏土及煤三相矿物混合，包含的孔隙主要为刚性孔和柔性孔两种类型，分别设定其孔隙纵横比为0.220和0.025。

图 5为多种弹性参数随孔隙度和含水饱和度的变化趋势，其中孔隙度变化范围是0~0.2(箭头指向增大方向)，含水饱和度变化范围是0~1.0。可见横波速度随含水饱和度变化幅度较小，主要是受流体密度变化的影响；剪切模量的变化不受含水饱和度影响，是因为孔隙中流体饱和程度不会改变岩石骨架的性质；流体体积模量不受孔隙度影响，但随含水饱和度变化明显；除纵横波速度比外，其余弹性参数均随孔隙度增大而减小，且均受含水饱和度的影响较小。

图 6为不同弹性参数随泥质含量和煤层占比(V_coal)的变化趋势，其中煤层占比变化范围是0~0.1(箭头指向增大方向)，泥质含量变化范围是0~0.5。可见泥质含量固定(即颜色相同)时，随着煤层占比的增大，拉梅常数、纵波阻抗、纵波速度、弹性阻抗及横波速度等参数均呈明显减小趋势；泊松阻抗^[33](泊松比与密度的乘积)受煤层占比影响较小，变化幅度微弱。图中同一条虚线上的点具有相同的煤层占比，从0开始沿箭头方向以0.025的幅度增大。当煤层占比固定时，随着泥质含量的增加(即由蓝变黄)，除了泊松阻抗和拉梅常数增大外，其余弹性参数均减小。

在煤系地层中，储层预测受煤层的影响较大，煤层的分离与砂岩的识别一直是煤系地层油气勘探的难点与重点，因此需优选出较敏感的弹性参数区分煤层和储层。通过分析，主要研究纵波速度、横波速度、剪切模量、小角度(6°)弹性阻抗、泊松阻抗、纵横波速度比等参数对煤层占比和泥质含量的敏感度。

图 7是各弹性参数随煤层占比和泥质含量的变化情况。斜率越大，说明该弹性参数对煤层占比越敏感；颜色跨度越大，证明该弹性参数受泥质含量影响越大，变化越明显。进一步分析可知，6种弹性参数中，剪切模量对煤层敏感性最高，泊松阻抗受煤层占比影响十分微弱，但对泥质含量变化敏感，据此可分别作为区分煤层和砂岩的敏感因子。

在实际应用中，将各弹性参数的敏感性系数^[34]定义为，表达式为

$ C=\left|\frac{P_{\text {mean }}-P_{\text {mean }}^{0}}{P_{\mathrm{std}}}\right| $

(12)

式中：C为敏感系数，数值越高说明对应的弹性参数越敏感；P_mean是目标岩性的弹性参数均值；P_mean⁰是其他岩性的弹性参数均值；P_std是目标岩性的弹性参数标准差。

为进一步验证上述结论的正确性，利用东海盆地M区实际测井资料进行弹性参数敏感性计算测试。当目标岩性为砂岩时，背景岩性为泥岩和煤层，当目标岩性为煤层时，背景岩性为砂岩和泥岩。

图 8是泊松阻抗、拉梅参数、纵横波速度比、剪切模量、横波速度、弹性阻抗、密度等不同参数对煤层和砂岩的敏感性分析。可以看出：泊松阻抗对砂岩更敏感，可作为储层预测的敏感参数；剪切模量对煤层更敏感，可用于识别煤层。

3 横波速度估算

针对煤系地层，基于岩石物理模型，可定量表征煤系地层弹性参数与物性参数之间的关系，其中纵、横波速度与体积模量、剪切模量间的关系为

$ V_{\mathrm{P}}=\sqrt{\left(K+\frac{4}{3} \mu\right) / \rho} $

(13)

$ V_{\mathrm{S}}=\sqrt{\frac{\mu}{\rho}} $

(14)

式中ρ为密度。

实际应用中，横波速度实测资料较少，常规岩石物理建模时未考虑煤层的影响，估算速度时无法得到煤层对应的低值异常，难以适用于煤系地层。以M工区实际测井X为例，该井测得的地质背景曲线如图 9所示，可见煤层处具有较明显的低速度、低密度特征。利用所建模型直接计算纵、横波速度，图 10是建模时考虑煤层与否得到的测井X纵、横波速度对比。从绿色虚线框可看到，在煤层处考虑煤层后建模计算的速度曲线(黑线)与实际测井数据(蓝线)吻合较好，常规建模计算的纵、横波速度曲线(红线)明显高于实测纵、横波速度(蓝线)，证明建模时已有效地将煤层加入到地震岩石物理模型中。

由于岩石弹性参数受多种因素影响，很难建立其与物性参数的显式表达，本文利用所建模型构建了一种煤系地层弹性模量非线性预测方法。以实测纵波速度为约束，密度、孔隙度、含水饱和度、泥质含量等测井数据依据所建模型的岩石物理关系参与计算，对测井解释物性参数与纵波速度不相匹配的位置做局部修正。最终目标是利用煤系地层地震岩石物理模型及局部修正的储层物性参数，计算井曲线中常缺失的横波速度等参数。利用模拟退火算法逐点修正模型参数，建立目标泛函

$ \begin{aligned} \boldsymbol{m}^{*}=& \operatorname{argmin} \| V_{\mathrm{P}}^{\text {mea }}-\\ & V_{\mathrm{P}}^{\text {cal }}\left(\boldsymbol{K}, \boldsymbol{\mu}, V_{\text {sh }}, V_{\text {coal }}, \rho, \phi, S_{\mathrm{w}}, \boldsymbol{a}\right) \|_{2}^{2} \end{aligned} $

(15)

式中：m^*=[V_sh，V_coal，ϕ，a]为物性参数，其中V_sh、V_coal分别为泥质含量和煤层占比，a为不同孔隙类型的孔隙纵横比；V_P^mea为测量的纵波速度；V_P^cal(K，μ，V_sh，V_coal，ρ，ϕ，S_w，a)为实际计算得到的纵波速度；K、μ分别为各组分的体积模量和剪切模量；S_w为含水饱和度。模拟退火算法跳出了局部搜索的局限性，可在解空间内趋于全局最优。本文反演仅是对实际测井解释物性参数进行修正，是分参数逐个反演的，以孔隙度和泥质含量为主。各物性参数进行反演时，参数每扰动一次，就会随机产生一种新状态，以Metropolis准则判断是否接受新状态，最终通过大量次数迭代得到满足终止条件的解(图 11)。

将本文方法应用于东海盆地M区，对测井Z进行纵波速度和横波速度估算。模型中孔隙度类型设定为泥质孔隙(软孔)和砂质孔隙(硬孔)两种，对应的孔隙度为ϕ_c=ϕV_sh、ϕ_s=ϕ-ϕ_c(其中ϕ_c为泥质孔隙，ϕ_s为砂质孔隙)，输入参数包括泥质含量、含水饱和度、纵波速度、横波速度、孔隙度、密度、煤层占比、不同矿物组分及所含流体的模量等，对比并分析模型预测效果。

在纵波速度约束下，利用模拟退火算法对孔隙度、泥质含量、孔隙纵横比、煤层占比等物性参数做修正。其中，硬孔和软孔的孔隙纵横比是根据实际工区储层发育特征预先分别给定的，如0.220和0.025。对比图 12中的修正曲线(红线)与实际数据(蓝线)，可见除孔隙度和泥质含量曲线进行了轻微修正之外，其余数据与原始曲线高度近似，也进一步证明了所构建模型的有效性。

图 13为测井Z纵、横波速度预测结果与实际测井曲线的对比。可看出纵波速度和横波速度的计算结果(红线)与真实测井曲线(蓝线)的变化趋势基本一致，且在煤层低值处具有较高的吻合度。

图 14为纵、横波速度相对误差分析。可见纵波速度误差在0.05以内，横波速度误差集中在0.10以内，整体速度误差较小且呈正态分布，证明了所构建的地震岩石物理模型在煤系地层的适用性。

4 结论

地层含煤对岩石弹性模量的影响较大，在建模过程中不可忽略。本文根据实际地震资料及地质背景，构建了一种煤系地层地震岩石物理模型；充分考虑地层含煤对地震响应的影响，分析微观物性参数对弹性参数的作用。理论分析和实例研究均表明：剪切模量对煤层变化敏感，可作为煤层的指示因子；泊松阻抗随煤层变化幅度小，但受泥质含量影响大，可为后续储层反演做贡献。

根据相关地震岩石物理理论计算，可先后得到背景矿物、干岩石骨架及含流体孔隙介质的体积模量和剪切模量。在实际测井资料约束下，构建反演目标泛函，形成一种煤系地层岩石横波速度非线性计算方法，提高了储层岩石弹性模量预测的精度，完善了不同储层类型的横波速度计算方法。将该方法应用于实际工区测井的速度预测，预测结果与测井数据吻合度较高，进一步验证了本文方法的有效性和实用性。