0 前言
在页岩储层裂缝参数及含气性预测中，叠前地震反演是常用技术，该技术通过估计岩石弹性性质或借助岩石物理模型建立储层物性参数之间的定量联系，实现储层物性参数的直接反演[6-9]。Zong等[10]提出了一种新的正交各向异性介质反射率模型参数化方法，结合贝叶斯理论开发了实用的AVAZ反演方法，该方法增强了正交各向异性介质反演的稳定性和可靠性，为含垂直或近垂直裂缝的页岩储层裂缝属性估算及油气预测提供了有效技术支持。陈珂磷等[11]使用模型参数化的方法推导了HTI介质反射系数近似表达式，并基于贝叶斯理论构建AVAZ反演方法，采用柯西概率分布和高斯概率分布作为模型参数与似然函数的先验信息，结合奇异值分解进行参数去相关处理，引入平滑背景约束以提高反演的稳定性和精度。通过与现有近似方程的精度对比及不同噪声条件下的实际工区应用，此方法在页岩裂缝型储层反演预测中具有可行性、高鲁棒性与稳定性，能够为裂缝发育区域的精准表征提供有效技术支撑。为降低反演过程中独立参数的个数并提高反演的稳定性，李红梅等[12]通过变换域方法挖掘宽频地震数据的低频信息构建各向异性参数初始模型，采用方位振幅差异反演策略将六参数直接反演转化为二次三参数反演，并构建了综合矿物脆性和裂缝特征的各向异性脆性指示因子，提升了HTI介质反演的稳定性，为页岩储层的脆性预测和甜点评价提供了有效技术支撑。然而，含气填充的岩石中大都存在几何形状各异的孔隙，造成储层物性参数与地震反射系数方程间存在较强的非线性关系，上述储层弹性参数与物性参数反演的过程中未能充分考虑孔隙形状的复杂性，使得储层裂缝检测及流体预测容易出现较大的误差。
为了进一步在储层预测的过程中考虑孔隙形状对地震振幅的影响，潘辉[13]推导了包含等效孔隙纵横比与流体体积模量的反射系数方程，建立了地震数据、孔隙度、等效孔隙纵横比以及流体识别敏感参数间的数学关系，基于贝叶斯框架的MCMC方法，实现了等效孔隙纵横比、孔隙度及流体体积模量的地震反演预测。该方法使用简化参数描述孔隙形状的特征，既能提高反射系数方程在复杂储层应用的合理性，也为基于地震资料预测储层的孔隙空间奠定了基础。然而，该方法未对储层裂缝发育时存在各向异性的情况进行有效解释，因此难以对裂缝型页岩储层进行裂缝参数的预测。
针对上述裂缝型含气页岩储层预测过程中对孔隙形状影响考虑不足以及储层参数预测稳定性的问题，本次研究将重点考虑孔隙形状这一关键因素，旨在构建一套高精度的页岩储层裂缝预测与流体识别方法。首先，通过测井数据优选适用于目标靶区含气页岩储层刻画所需的流体识别因子。之后，结合岩石物理建模与各向异性介质理论，推导考虑孔隙形状的各向异性反射系数方程，建立地震反射特征与储层裂缝参数与流体识别因子的定量关系。在此基础上，依托贝叶斯框架，发展了一种裂缝参数与流体识别因子的两步叠前地震各向异性反演方法用于提高待反演参数的准确性与稳定性。最后，将所提出方法应用于合成记录测试与实际工区测试，验证了方法的有效性与适用性。
1 目标靶区流体识别因子优选
以多口探井的主力产气层段为对象，开展测井岩石物理参数的交会分析，以优选适用于目标靶区的流体识别因子。图1展示了多种含气页岩储层识别因子的交会情况。图1中的Gassmann流体项的定义为 $f=\rho {V}_{p}^{2}-{\gamma }_{dry}^{2}\rho {V}_{s}^{2}$ 。其中的Vp、Vs和ρ分别代表地层纵波速度、横波速度以及密度， ${\gamma }_{dry}^{2}$ 代表地层干岩石骨架的孔隙纵横比。由图1可知，在本次研究的靶区，Gassmann流体项f [14]在区分含气页岩与含水页岩的效果显著优于其他3个常用的流体识别因子（纵波阻抗、泊松比以及Vp/Vs）：图1(a)中的黑圈内的气层位置对应了f的低值异常，且f在气层和水层分别所对应的参数范围具有较为明显的差异。因此，本文采用Gassmann流体项f作为含气页岩的识别因子，开展考虑孔隙形状影响的页岩储层流体识别及裂缝预测研究。
2 考虑孔隙形状影响的页岩储层各向异性反射系数方程
考虑到目标靶区的页岩储层多发育高角度裂缝，因此，该地层可以近似成为具有水平对称轴的横向各向同性（horizontal transverse isotropy, HTI）介质。考虑到HTI 介质的各向异性项由裂缝法向弱度和切向弱度所表征，因此基于该参数建立的反射系数方程，可描述地震波在反射界面两侧的能量分配规律。当纵波入射、纵波反射时，反射系数 ${R}_{pp}^{HTI}\left(\theta,\varphi \right)$ 的表达式如下[15]：
式中： ${R}_{pp}^{HTI-iso}\left(\theta \right)$ 与 ${R}_{pp}^{HTI-ani}\left(\theta,\varphi \right)$ 分别表示各向同性背景项与各向异性裂缝项；θ与 $\varphi $ 分别表示入射角和方位角；Mb、μb和ρb分别表示岩石各向同性背景的纵波模量、剪切模量和密度； ${\delta }_{N}$ 和 ${\delta }_{T}$ 分别表示裂缝的法向弱度和切向弱度。考虑到Δ表示各参数在上、下界面处的差值，为了方便后续在方程的各向同性背景项中引入流体敏感参数，将式(1)所示的反射系数方程在反射界面两侧地层性质变化较缓的假设条件下，可以进一步拓展为：
式中：下标1和2分别表示反射系数方程中反射界面的上层介质与下层介质中的各参数。
为了进一步在反射系数方程的各向同性背景项中考虑孔隙形状影响并引入流体识别因子，依托岩石物理建模，建立考虑孔隙形状的VTI介质岩石物理模型，进一步表征岩石弹性参数与流体识别因子之间的定量关系。岩石物理建模方法如下：
第一步，确定岩石基质的模量。考虑到组成目标靶区页岩岩石基质的矿物主要为黏土矿物和石英，因此岩石基质的模量可由矿物混合物的总模量来定义，并由V-R-H平均[16]计算如下：
式中：Km与μm分别代表岩石基质的体积模量与剪切模量；C代表泥质含量；Kc、Kq、μc和μq则分别代表黏土体积模量、石英体积模量、黏土剪切模量和石英剪切模量。
第二步，计算考虑孔隙形状影响的干岩石骨架模量。考虑到实际储层中发育孔隙的扁度（纵横比）不同，可将其干岩石骨架等效为具有不同孔隙纵横比的多孔介质。为了进一步量化孔隙形状对干岩石骨架模量的影响，基于单孔介质岩石物理理论，将不同孔隙比的孔隙等效为具有单一孔隙纵横比的单孔介质，并计算考虑孔隙形状的干岩石骨架模量[17]如下：
式中：Kdry和μdry分别为干岩石骨架的体积模量和剪切模量； $\varphi $ 为孔隙度；p*和q*分别为反映孔隙形状的几何因子。为了进一步明确表征孔隙形状的几何因子对干岩石骨架模量的影响，本次研究绘制了如图2所示的不同p*和q*情况下的干岩石骨架模量随孔隙度的变化情况。可以看出，随着孔隙度的增加，干岩石骨架的模量逐渐降低，该现象是由于岩石孔隙含量的增加会破坏岩石基质的刚性结构，使得岩石更容易受外力作用而变形。在此基础上，随着孔隙形状几何因子的增大，干岩石骨架模量呈现降低趋势，孔隙结构因子的变化间接反映了孔隙纵横比的变化。随着孔隙结构因子增大，等效单一孔隙的纵横比变小，孔隙更容易被压缩，因此干岩石骨架刚性变弱，其模量减小。分析结果进一步证明了考虑孔隙形状的必要性与合理性，为叠前地震反演提供了更为可靠的岩石物理理论基础。
式中：Ksat与μsat分别为饱和岩石的体积模量与剪切模量；f为Gassmann流体项；Kf为孔隙流体体积模量。
式中： ${M}_{b}={K}_{sat}+\frac{4}{3}{\mu }_{sat}$ ， μb=μsat， λb=Mb-2 μb， χ=λb/Mb=1-2μb/Mb， g=μb/Mb。
经过上述四个岩石物理建模步骤，即可得到考虑孔隙形状的饱和HTI介质岩石物理模型。该模型阐明了HTI介质孔隙形状与岩石刚度之间的关系，为基于岩石物理反演求取孔隙形状几何因子奠定了基础。
为了验证岩石物理模型的合理性，从研究区的一口探井中选取了一段进行岩石物理建模正演验证，并计算孔隙形状几何因子，结果如图3所示。可以看出，预测的纵波速度和横波速度与实测曲线具有较好的吻合度，验证了所建立的岩石物理模型能够较好地建立储层速度与物性参数之间的关系，具有较好的合理性。岩石物理模型中包含的超参数（p*和q*）的计算将通过岩石物理建模反演得到。模拟退火算法用于模型中超参数的优化求解，该算法通过迭代调整超参数的数值来最小化由岩石物理模型正演预测得到的纵波速度与实际测井解释的纵波速度之间的误差。当误差在迭代后达到一定程度时，算法将停止迭代，并输出孔隙形状几何因子。上述研究结果为所建立的考虑孔隙形状影响的HTI岩石物理模型的准确性和实用性提供了强有力的理论验证，证实了该模型在地震反演中的适用性与合理性。
各向同性背景岩石模量与岩石物性参数以及流体识别因子之间的换算关系可由式(8)至式(13)来获得。将式(8)至式(11)带入式(12)与式(13)中，可得：
考虑到 ${M}_{b}={K}_{sat}+\frac{4}{3}{\mu }_{sat}$ ，同时，μb=μsat，最终可得各向同性岩石背景纵波模量、横波模量、岩石物性参数和流体识别因子之间的定量关系如下：
将式(17)与式(18)带入式(7)，以下角标1和2分别代表反射界面的上层介质与下层介质，即可得到考虑孔隙形状影响的页岩储层各向异性反射系数方程，式(19)建立了HTI介质地震反射系数与入射角、方位角、流体识别因子及裂缝参数之间的非线性直接联系，为叠前地震反演预测地层流体分布以及裂缝发育情况奠定了理论基础。
表1 双层裂缝型页岩与泥岩模型参数Table 1 Model parameters of dual-layer fractured shale and mudstone |
| f/（1013 Pa） | | C | ρb/(kg·m-3) | | | |
|---|---|---|---|---|---|---|
| 含气页岩 | 1.65 | 0.15 | 0.30 | 2 612 | 0.18 | 0.11 |
| 含水页岩 | 2.06 | 0.15 | 0.30 | 2 642 | 0.14 | 0.08 |
| 泥岩 | 2.37 | 0.05 | 0.90 | 2 704 | 0.10 | 0.05 |
3 叠前地震各向异性反演
式(19)具有6个独立的参数（f、 $\varphi $ 、C、ρb、 ${\delta }_{N}$ 和 ${\delta }_{T}$ ），因此，需要采用方位角作差的反演策略降低多个参数同时反演中独立参数的数量，以提高反演结果的准确性与稳定性。反演目标函数的构建将基于贝叶斯框架，以实现综合考虑先验信息的约束以及参数之间的相关性信息。当似然函数服从高斯分布 ${p}_{Gauss}\left(S\text{'}\left|m\text{'}\right.\right)$ ，先验信息服从柯西分布pCauchy(m')时，后验概率 $p\left(m\text{'}\right|S\text{'})$ 可以计算如下：
在叠前地震各向异性反演的实际应用中，根据贝叶斯反演理论，似然函数旨在通过背景噪声的概率密度分布来描述反演结果m'与实际地震数据S'之间的匹配程度，先验分布用于描述待反演的模型参数m'的先验信息，后验分布则是将似然函数与先验分布的概率密度同时考虑，用于建立先验模型约束下的叠前反演目标泛函。根据式(20)的描述，似然函数表达式如下：
式中： ${\sigma }_{n}^{2}$ 为噪声的方差；G'为正演算子。同时，先验分布的表达式描述如下：
式中： ${\sigma }_{m}^{}$ 、 ${\sigma }_{m}^{2}$ 分别为待反演参数的标准差、方差；N为待反演模型采样点的个数。结合式(21)与(22)，考虑多个方位角度作差情况，建立贝叶斯框架下的第一步反演裂缝参数 ${\delta }_{N}$ 与 ${\delta }_{T}$ 的目标函数如下：
式中：u代表第u个入射角；U代表参与反演的入射角的总个数；l代表第l个方位角作差；L则代表方位角作差的总个数。为了进一步提高反演结果的稳定性，将式(23)加入平滑背景约束对地震数据的低频缺失作以补充，可得最终的第一步反演的目标函数如下：
在第一步反演目标函数建立的基础上，同样基于贝叶斯框架并结合平滑背景约束，进而构建第二步反演流体识别因子等各向同性背景项的目标函数如下：
式(23)至(30)中：W表示子波矩阵；m代表待反演参数；ΔS与S分别代表方位角度作差之后的地震数据和原始地震数据； ${\Sigma }_{\Delta S}^{}$ 与 ${\Sigma }_{S}^{}$ 分别表示方位角度作差之后和原始地震数据的方差。为求解所建立的非线性目标函数，本文采用启发式算法-遗传算法[22-23]对目标函数进行迭代求解。遗传算法通过模拟自然界的选择与遗传过程实现优化求解，具体涉及复制、交叉和突变等操作流程。该算法以初始种群作为运算起点，借助遗传算子的迭代操作生成更适应环境的个体集合，促使种群在搜索空间中逐步向最优群体进化。在此过程中，种群通过不断繁殖与演化，最终收敛为对环境适应性最强的群体，进而得到目标函数收敛时所对应的最优解（反演结果）。相比于常规梯度下降等算法，启发式算法在目标函数寻优时，不需要计算复杂的目标函数的梯度项，同时，由于其每次迭代更新的随机性，因此能够避免寻优过程中陷入局部极小值，对复杂非线性问题的求解具有更强的适应性。
4 合成记录测试
为系统测试本文所提出的方法的合理性与可行性，基于实际页岩作业靶区的测井解释结果，建立了如图5所示的含储层裂缝参数与流体识别因子的待反演地层参数模型，开展合成记录测试工作。
以图5所示的模型作为输入数据，基于式(19)所示的反射系数方程计算不同方位角度与不同入射角度下的反射系数方程，在此基础上结合褶积模型进一步计算并输出如图6所示的无噪声合成地震记录。该地震记录反映了地层参数与地震振幅信息之间的非线性关联，能够为使用叠前地震各向异性反演地层参数提供数据基础。为了进一步在合成记录测试中评价反演方法的抗噪性与稳定性，将图6所示的合成地震记录添加了信噪比（signal to noise ratio, SNR）为10的随机噪声以开展含噪声合成地震记录下的反演测试工作。含信噪比为10的随机噪声的合成地震记录如图7所示。分别将图6与图7所示的合成地震记录作为观测数据输入，基于方位角度作差的两步叠前地震各向异性反演方法，分别输出了如图8中蓝色曲线与青色曲线所示的无噪声和SNR=10的反演结果。
图7 含信噪比为10的随机噪声的合成地震记录Fig. 7 The Synthetic seismic record containing random noise with a SNR of 10 |
观察反演结果不难发现，由于采用了两步地震反演降低了每一步反演过程中独立参数的个数，并且采用了贝叶斯框架，同时考虑了先验信息约束与待反演参数的相关性信息，更重要的是反演结果考虑了孔隙形状对地层的影响，使得反演结果具有较高的合理性与稳定性，反演能够对地层的裂缝参数与流体识别因子的纵向变化进行较为准确的预测。对比无噪声反演结果与SNR=10合成地震记录下的反演结果，虽然SNR=10的反演结果的准确性有所降低，但仍然能够保持较高的精度，具有良好的抗噪声能力，能够为该方法在实际页岩作业靶区的应用提供技术支持。为了定量化表征反演结果与模型值之间的差异，图9展示了图8中的反演结果（蓝色曲线与青色曲线）与模型值（红色曲线）之间的误差绝对值，可以看出，无噪声反演结果具有较小的误差绝对值，SNR=10合成地震记录下的反演结果的误差绝对值受观测数据准确性的影响有所增加，但仍然保持在合理的范围内。误差绝对值的展示从量化的角度进一步验证了反演结果具有较好的合理性与准确性。
合成记录测试验证了该方法的合理性与理论可行性，为本方法应用于实际页岩作业区开展裂缝预测与流体识别提供了理论依据和方法支撑。
5 实际工区测试
选取中国黔北地区道真向斜五峰组—龙马溪组地震资料，开展考虑孔隙形状影响的储层裂缝预测与流体识别方法的应用实践。
图10为用于实际工区测试的实际页岩气田作业区的过井地震剖面。其中，井1作为已知井井参与反演并提供反演所需的先验信息，井2则作为盲井用与反演结果的验证。在该地震剖面所示的工区中，工区目的层下部（约4 km处）发育一套主要的含气页岩，且高角度裂缝发育。地震数据经过了严格的精细处理，保证了振幅随入射角度和方位角度的精确变化关系，为叠前各向异性地震反演的开展提供了良好的数据基础。孔隙形状因子可由岩石物理模型结合实测测井曲线进行逆推预测得到。以图10所示的地震数据作为观测数据输入，首先进行方位地震数据作差并开展第一步裂缝参数（ ${\delta }_{N}$ 和 ${\delta }_{T}$ ）的反演，输出的考虑孔隙形状影响和未考虑孔隙形状影响的反演结果分别如图11a和图11b所示。其中，分图中的测井位置处，本文采用了变密度的形式展示了裂缝参数的滤波后测井曲线。
可以看出，裂缝参数的反演结果与已知井和盲井的测井解释结果都具有较高的吻合程度，考虑孔隙形状影响的反演结果相比于未考虑孔隙形状的反演结果具有更好的准确性与可靠性。考虑孔隙形状影响的反演结果在纵向上具有较高的分辨率，横向上具有较强的连续性，能够满足地层性质刻画的需要。通过观察黑色圆圈位置处的裂缝法向弱度和切向弱度的高值异常可以发现，裂缝参数的反演结果能够对含气页岩目的层的高角度裂缝的发育状态进行有效地预测，并对该套目的层展布特征进行精细地刻画，预测结果具有较好的合理性。而相应位置处未考虑孔隙形状的反演结果与测井解释结果的吻合度则相对较低，对于高角度裂缝的精细刻画效果则相对较差。
在此基础上，分别将图11a和图11b所示的裂缝参数反演结果作为已知信息，同样以图10所示的地震数据作为观测数据输入，进行第二步地层各向同性背景参数f、 $\varphi $ 、C和ρb的反演。考虑孔隙形状影响和未考虑孔隙形状影响的预测结果分别如图12a和图12b所示。由于分步反演的方法减少了第二步反演过程中的独立参数的个数，且贝叶斯框架引入了先验信息约束和待反演参数的相关性信息，与此同时，反射系数方程纳入了孔隙形状的影响，因此基于该方程得到的流体识别因子f等各向同性背景参数的反演结果与井位置处滤波后的测井解释结果吻合程度更高。而且反演结果的分辨率能够满足描述地层性质的纵向变化和横向分布状态，使得含气储层的刻画具有更高的分辨率。在图12a黑圈内的含气页岩主力目的层，考虑孔隙形状的流体识别因子f相比于未考虑孔隙形状的流体识别因子f具有更为明显的低值异常，可见该方法能够对含气页岩储层的展布进行合理、有效的预测，并对后续的页岩气开发提供参考价值。实际工区测试结果进一步证实了本次研究所提出方法在实际页岩储层勘探作业区的适用性，为页岩气的勘探提供了新的理论方法和技术支撑。
6 结论
本文围绕孔隙形状影响的页岩储层裂缝预测与流体识别展开系统研究，通过理论推导、方法构建与实例验证，取得了主要结论如下：
（1）流体识别因子优选结果分析表明Gassmann流体项在目标靶区区分含气与含水页岩的效果优于纵波阻抗、泊松比等常用敏感因子，其低值异常可有效指示气层位置，为含气页岩的预测提供了可靠依据。
（2）基于岩石物理建模，推导了考虑孔隙形状影响的HTI 介质各向异性反射系数方程。该方程建立了地震反射系数与入射角、方位角、流体识别因子f及裂缝参数的非线性关系，具有良好的入射角与方位角敏感性，为叠前各向异性反演提供了合理的正演算子。
（3）所提出的考虑孔隙形状的页岩储层裂缝预测与流体识别方法在合成记录测试中展现了较高精度与较强抗噪性，对黔北地区道真向斜五峰组—龙马溪组页岩气储层的裂缝参数反演预测结果与测井解释吻合度较高，可精准刻画高角度裂缝发育特征，流体识别因子反演能有效指示含气页岩目的层展布，在实际工区具有较强的适用性与有效性，为页岩气勘探开发提供了重要的技术支撑。
