火山岩储层地震勘探一直是地球物理界关注的重点之一, 伴随着油气勘探领域的不断拓展, 火山岩储层的勘探开发逐渐受到重视。张新荣等^[1]分析了火山岩储层作为复杂特殊油气储层在其勘探开发中面临的诸多难题; 宋宗平^[2]和张芝铭等^[3]应用基于宽带约束的模拟退火反演方法并结合地震数据、测井数据等获取火山岩纵波速度、横波速度、密度等弹性参数, 从而为地质研究人员提供重要的解释依据, 它已成为火山岩储层勘探开发的重要研究内容; 陈永波等^[4]通过古地貌恢复研究火山岩沉积; 朱超等^[5]通过叠后波阻抗反演预测火山岩优质储层平面展布特征; 杨晓光等^[6]采用分频属性反演和基于神经网络的储层参数反演有效预测了火山岩储层分布; 毛庆鑫^[7]认为火山岩结构复杂、岩性横向变化快、具有较强的各向异性特征, 在地震剖面上表现为同相轴横向连续性差、反射呈现低频和强振幅特征, 可见空白反射或杂乱反射结构。

然而, 基于叠后地震数据的地震属性(相干、波形聚类等)分析和叠后波阻抗反演等技术手段对火山岩机构内部期次划分、岩性和物性的预测存在不足, 越来越不能够满足火山岩勘探开发的需求^[8]。针对这一问题, 何辉等^[9]应用叠前方位各向异性法, 通过优选衰减起始频率属性, 综合预测了佳木河组火山岩储层裂缝分布特征; 张新培等^[10]利用测井响应和多井约束下的波阻抗反演方法实现了优质火山岩储层的识别和预测; 印兴耀等^[11]认为, 目前应用纵波速度、横波速度和密度参数等结合井参与的低频背景模型进行叠前同时反演来识别火山岩机构内部储层效果不佳的原因是, 地质模型与实际地质情况差异较大而导致反演结果准确率不高。如何弱化低频背景模型的影响, 使反演结果既符合地质规律又有较高的准确度是目前火山岩机构内部储层预测的关键所在。

马尔可夫随机场反演方法在储层建模领域早有应用, TARANTOLA^[12]最早提出贝叶斯框架下的概率反演方法, 基于马尔可夫链蒙特卡洛(markov chain monte carlo, MCMC)统计方法获得待反演参数的后验概率密度函数, 统计分析反演结果, 获取最优的反演解, 并对反演结果进行不确定性分析; SPIKES^[13]将多种蒙特卡洛抽样方法引入地质统计建模, 联合估算阻抗等弹性参数和孔隙度等储层属性; MARIT等^[14]使用马尔可夫随机场建立岩性/流体分类概率分布的先验模型, 基于MCMC方法获取弹性参数、孔隙度以及岩性/流体分类后验概率分布的预测和模型参数的不确定性评估; MICHEL等^[15]将分级贝叶斯模型引入MCMC方法; ODD等^[16]使用MCMC随机反演方法反演纵波阻抗; 陈裕华^[17]应用纵波、转换波振幅和旅行时数据, 反演孔隙度等储层物性参数; 田玉昆^[18]采用Huber-Markov随机场约束边缘保护方法, 求解叠前纵横波速度和密度三参数反演问题。可见, 在学者们的不断深入研究下, 马尔可夫随机场反演方法在最优解、去除边界效应等方面不断完善, 但大部分的研究重点是改进算法和模型试算上, 只有少部分研究应用到实际数据, 而且是碎屑岩和碳酸盐岩数据中。

针对火山岩地层内部地震反射波横向变化快、呈现中弱反射和非均质性强的特点, 目前仍没有较好的反演方法。冯玉辉等^[19]系统地进行了火山岩地震响应特征描述; 代春萌等^[20]进行了火山岩相控反演; 有研究者应用体序域建模识别火山岩储层, 对于没有井钻遇的火山岩反演结果与体序域模型具有较强的相似性。

本文将马尔可夫随机场引入到火山岩地震叠前同时反演算法中, 针对火山岩岩性多样且测井曲线数值叠合严重的特点, 绕开火山岩岩性的预测, 利用测井数据划分岩相, 对定义的每个岩相都计算出低频背景模型, 并将每个岩相应用马尔可夫随机场求取先验分布函数, 从而有效避免了反演结果由于地震能量差异和频率不足等引起的假象, 提高了预测精度。

1 叠前同时反演面临的问题 1.1 归一化问题

在应用叠前确定性同时反演方法进行岩石物理交会分析时, 需要确定纵波阻抗和横波阻抗的线性关系。图 1为典型的横波阻抗与纵波阻抗交会图。分析图 1可知, 图中线性关系需要3条线才能够表述清楚, 两条通过倾向砂岩的蓝色数据团(含水和含烃砂岩)以及一条通过倾向泥岩的橙色数据团。

1.2 低频模型问题

目前, 建立低频模型通常是将井的波阻抗曲线进行低通滤波和内插, 这种方法仅适用于较平缓的地层, 对于复杂的地下地质情况该方法主要存在以下两个方面的不足:

1) 对于高净毛比值(有效储层厚度与储层总厚度的比值)的井和低净毛比值的井, 由于无法确定某个位置的准确净毛比, 因此期望在两者之间内插获得可靠的低频背景模型不现实;

2) 低频模型的实现算法只包含针对平行于顶层、平行于底层、顶底层之间等比例模式内插的3种算法, 适用于平行地层、超覆和尖灭地层, 但是对于碳酸盐缝洞、火山岩内幕、复杂断块等复杂地层问题适用性较差, 反演结果会出现较大误差。

1.3 物理特性识别

图 2给出了纵波阻抗模型及其合成地震记录与叠前同时反演的纵波阻抗结果。图 2中黑色曲线代表纵波阻抗真实模型(横波阻抗和密度也是一样, 未显示), 红色曲线表示纵波阻抗叠前同时反演结果, 黑色虚线是计算纵波阻抗用到的低频背景模型, 右侧是零入射角的合成地震记录。分析图 2可知:合成地震记录中同相轴最大振幅值位置对应的是岩相变换的地方, 叠前同时反演结果只能在这些岩相变换位置附近提供可靠的结果, 远离这些岩相变换的位置, 纵波阻抗结果趋向于低频模型(如箭头所示)。

2 方法原理

在理想情况下, 希望用一种不依赖于低频背景模型的内插方法来得到阻抗信息。采用迭代方法, 先用地震数据反演出相(给定一个均匀的阻抗模型), 然后由反演的相得到地震数据的阻抗信息, 再在这些阻抗信息基础上, 由地震数据再反演出相, 迭代至收敛。为了实现这一反演策略, 引入马尔可夫随机场反演理论。

2.1 马尔可夫随机场原理

马尔可夫随机场(MRF)方法建立在马尔可夫随机模型(MR)和Bayes理论的基础上, 马尔可夫随机场模型提供了不确定性描述与先验知识联系的纽带, 并利用观测图像根据统计决策和估计理论中的最优准则, 确定分割问题的目标函数, 求解满足这些条件或目标函数的最大可能分布, 从而将分割问题转化为最优化问题^[21-22]。

MRF是基于图像中像素的灰度信息, 利用像素之间的互相作用与空间的关系, 实现图像分割^[23]。MRF模型由于自身的特点, 既能反映图像的随机性, 又能反映图像的空间结构, 从而能够有效描述特定图像的空间性质, 保证了目标函数有唯一的最小点; 在求解MRF描述的不确定性问题时, 将图像先验知识转化为先验概率分布, 基于贝叶斯理论, 采用先验分布函数和似然函数来描述最大后验估计的图像分布, 同时保留了边缘结构, 便于边缘识别。应用数学公式推导出的参数具有明确的物理含义, 而常规线性模型中的参数一般只作为被拟合了的参数而出现, 较少具有真实含义。

2.2 马尔可夫随机场的地震反演

将地震反演看作一个统计问题, 贝叶斯定理可以写成:

$ \pi \left( {\mathit{\boldsymbol{Z}}\mid {\mathit{\boldsymbol{S}}_{{\rm{real }}}}} \right) \approx L\left( {{\mathit{\boldsymbol{S}}_{{\rm{real }}}}\mid \mathit{\boldsymbol{Z}}} \right)p(\mathit{\boldsymbol{Z}}) $

(1)

式中:Z是波阻抗; S_real是实际地震数据; π(Z|S_real)是后验分布函数; L(S_real|Z)是似然函数; p(Z)是先验分布。

将分布表示为多正态形式, 通常比较合理, 并更容易进行数学处理。先验分布P(Z)可以从纵波阻抗的线性深度趋势、纵波阻抗与横波阻抗、纵波阻抗与密度之间的交会结果得到, 并通过井数据进行不确定性评估。这些趋势拟合可以表示为一个多正态先验分布P(Z), 形式为:

$ P(\mathit{\boldsymbol{Z}}) \approx \frac{{\exp \left[ { - \frac{{{{\left( {\mathit{\boldsymbol{Z}} - {\mathit{\boldsymbol{Z}}_0}} \right)}^{\rm{T}}}\mathit{\boldsymbol{C}}_{\rm{P}}^{ - 1}\left( {\mathit{\boldsymbol{Z}} - {\mathit{\boldsymbol{Z}}_0}} \right)}}{2}} \right]}}{{{{\left| {{\mathit{\boldsymbol{C}}_{\rm{P}}}} \right|}^{\frac{1}{2}}}}} $

(2)

式中:C_P是描述阻抗方差和相关性的协方差矩阵; Z₀是波阻抗的初始模型。

似然函数L(S_real|Z)可以表示为:

$ \begin{array}{*{20}{l}} {L\left( {{\mathit{\boldsymbol{S}}_{{\rm{real }}}}\mid \mathit{\boldsymbol{Z}}} \right) \cong }\\ {\qquad \frac{{\exp \left\{ { - \frac{{{{\left[ {{\mathit{\boldsymbol{S}}_{{\rm{real }}}} - \mathit{\boldsymbol{F}}(\mathit{\boldsymbol{Z}})} \right]}^{\rm{T}}}\mathit{\boldsymbol{C}}_{\rm{d}}^{ - 1}\left[ {{\mathit{\boldsymbol{S}}_{{\rm{real }}}} - \mathit{\boldsymbol{F}}(\mathit{\boldsymbol{Z}})} \right]}}{2}} \right\}}}{{{{\left| {{\mathit{\boldsymbol{C}}_{\rm{d}}}} \right|}^{\frac{1}{2}}}}}} \end{array} $

(3)

式中:F(Z)是从阻抗Z导出合成地震记录的函数; C_d是表示地震噪声的协方差矩阵。

理想情况下, 为了准确描述地震反演问题的物理特性, 应对每个岩相和阻抗进行反演, 并且每个相的岩石物理模型和低频模型都应该使用, 这样做在数学上的要求很高, 同时也不能用标准的优化公式来解决, 但可以使用迭代方法解决这一问题:

首先给定一个均匀的阻抗模型从地震数据中反演出相, 然后从这些相中反演出阻抗信息, 再在阻抗信息基础上从地震数据中再反演出相, 迭代至收敛^[18-19]。在这种情况下, 公式(1)可以写成:

$ \pi \left( {\mathit{\boldsymbol{Z}}, \mathit{\boldsymbol{F}}\mid {\mathit{\boldsymbol{S}}_{{\rm{real }}}}} \right) \approx L\left( {{\mathit{\boldsymbol{S}}_{{\rm{real }}}}\mid \mathit{\boldsymbol{Z}}} \right)p(\mathit{\boldsymbol{Z}}\mid \mathit{\boldsymbol{F}})p(\mathit{\boldsymbol{F}}) $

(4)

式中:π(Z, F|S_real)是后验分布函数; p(Z|F)为岩相的先验分布; p(F)为岩相的概率密度。

先验分布p(Z|F)与(2)式非常相似, 唯一的区别在于先验均值Z₀和协方差矩阵C_P取决于:

$ \begin{array}{*{20}{l}} {p(\mathit{\boldsymbol{Z}}\mid \mathit{\boldsymbol{F}}) = }\\ {\quad \frac{{\exp {{\left\{ { - \frac{{\left[ {\mathit{\boldsymbol{Z}} - {\mathit{\boldsymbol{Z}}_0}(\mathit{\boldsymbol{F}})} \right]}}{2}} \right\}}^{\rm{T}}}{{\left[ {{\mathit{\boldsymbol{C}}_{\rm{P}}}(\mathit{\boldsymbol{F}})} \right]}^{ - 1}}\left[ {\mathit{\boldsymbol{Z}} - {\mathit{\boldsymbol{Z}}_0}(\mathit{\boldsymbol{F}})} \right]}}{{{{\left[ {{\mathit{\boldsymbol{C}}_{\rm{P}}}(\mathit{\boldsymbol{F}})} \right]}^{\frac{1}{2}}}}}} \end{array} $

(5)

对于岩相的概率密度p(F), 应用离散马尔可夫随机场求取^[24-25]。简单地说, 任何三维地震网格都是由许多像素组成, 定义一系列像素对, 在反演时窗内的每个像素能直接与相邻像素相连接, 因此每个像素属于6个像素对(边缘、角落除外)。然后, 整个网格的配置F的概率由所有像素对上的势能之和定义, 即所谓的吉布斯分布:

$ P(\mathit{\boldsymbol{F}}) \approx {\rm{exp}}\left[ { - \sum {{V_{\rm{p}}}({\mathit{\boldsymbol{F}}_{\rm{p}}})} } \right] $

(6)

其中, V_p表示每组像素对上的相系列F_p的“势能”。

对于每个像素团, 可以将V_p写为V_p=β·I(F_centre, F_neighbour), 其中, F_centre为中心相, F_neighbour为临近相。如果F_centre和F_neighbour相同, 则离散指标函数I为0, 如果它们不相同则为1, 并且β为一个正连续参数。所以, 公式(6)可改写为:

$ P(\mathit{\boldsymbol{F}}) \approx {\rm{exp}}\left[ { - \sum {\beta \mathit{{I}}({F_{{\rm{cengtre}}}}, {F_{{\rm{neighbour}}}}} } \right] $

(7)

因此, 只要两个相邻相不同, 就会减少概率P(F)。由于水平方向(主测线和联络线方向)和垂直方向的连续性不同, 因此, 需要对其赋予不同的β值。然而, 由于地质体很少是完全水平的, 使用年代地层的年龄概念来确定的相同年龄的邻点, 有可能在时间坐标上并不是邻近的像素。马尔可夫随机场可以阻止不合理的岩相组合, 即可以确保不会在含气砂岩的顶部得到含水砂岩。与常规叠前同时反演方法的不同之处在于增加了每个相的深度趋势和交会图进行约束, 而没有使用变差函数。

3 模型数据分析

基于马尔可夫随机场的反演方法是贝叶斯同时反演方法的一种扩展。两种方法之间的最大区别在于:在贝叶斯形式的同时反演中, 纵波阻抗与深度是一个线性趋势; 而在马尔可夫随机场的反演中可以反演两个或者多个岩相, 所以可以指定两个或更多这样的趋势。这意味着在反演时窗内, 阻抗结果与每个相趋势之间都可以建立联系。

为了验证基于马尔可夫随机场反演算法的有效性, 建立了一个如图 3a所示砂岩侵入泥岩的楔状模型。模型频宽0~60Hz, 并应用25Hz雷克子波进行反演, 为了更好地进行对比分析, 只考虑纵波速度的变化。图 3b给出了当入射角为0时的地震正演结果。对正演结果进行常规叠前同时反演和基于马尔可夫随机场反演, 两种反演结果分别如图 3c和图 3d所示。分析可知:在常规叠前同时反演中使用了一个低频背景模型, 其纵波速度值在泥岩纵波速度值和砂岩纵波速度值之间; 基于马尔可夫随机场的反演使用了两个低频趋势, 其纵波速度值分别为泥岩纵波速度值和砂岩纵波速度值, 从反演结果对比来看, 基于马尔可夫随机场的反演结果与模型更加吻合, 预测准确率更高。

4 实例分析

将同时反演方法和基于马尔可夫随机场的反演方法应用于长岭龙凤山地区中基性火山岩油气储层的预测。由测井数据分析可知, 研究区目的层系储层以凝灰岩为主, 因此定义了3个相来进行反演:将孔隙度大于5%且泥质含量小于40%的凝灰岩对应的波阻抗定义为凝灰岩高孔隙相; 将孔隙度小于5%且泥质含量小于40%的凝灰岩对应的波阻抗定义为凝灰岩干层相; 将泥质含量大于40%的岩性对应的波阻抗定义为泥岩相。

由密度概率分布图(图 4)可以看出, 密度能够对储层进行有效识别, 因此选取密度参数作为敏感参数参与到后续反演之中。

研究区目的层段反演结果如图 5所示。其中, 图 5a为常规叠前同时反演得到的密度数据体; 图 5b为基于马尔可夫随机场反演得到的岩相体与地震剖面叠合显示结果, 图 5b中红色代表凝灰岩高孔隙相, 黄色代表凝灰岩干层相, 灰色代表泥岩相; 图 5c为基于马尔可夫随机场反演得到的密度数据体。图 5a中反演结果在井眼处与密度曲线对应较好, 但同一火山岩机构之间的两口井横向连续性和整体期次差异较大; 由于划分了3个岩相, 进行了3个低频模型约束, 图 5b中反演的岩相体在横向展布上与地震资料更加吻合, 在纵向油气水配置关系上更符合石油地质理论和实际地质情况; 图 5c中的结果不仅在井眼处与密度曲线对应较好, 而且同一火山岩机构之间的两口井横向连续性更好, 能够较好地识别火山岩期次, 同时也能够对储层物性进行较好的预测, 可指导下一步井位部署。本文方法在松辽盆地龙凤山地区火山岩的勘探开发中发挥了重要作用, 依据预测结果部署了3口直井均获得工业气流, 4口水平井均获得了高产气流, 以预测结果为依据在该地区提交火山岩天然气探明储量102×10⁸m³。

5 结论

本文将马尔可夫随机场引入到火山岩地震叠前反演算法中, 其优点在于:采用火山岩多条测井数据交会划分岩相, 且测井数据既不直接参与低频模型的建立也不直接参与反演, 该方法对定义的每个岩相都计算出低频背景趋势, 使反演结果有效避免了由于地震能量差异、频率不足等引起的假象。

从马尔可夫随机场反演结果来看:

1) 马尔可夫随机场反演方法不依赖传统的低频模型, 反演结果更遵从于地震数据;

2) 反演得到的岩相体符合石油地质学规律, 对地质研究具有很大帮助;

3) 反演得到的密度体与井震标定后的测井曲线有较好的匹配关系, 可以用于该地区火山岩机构内部储层的物性预测。

该反演方法适合于火山岩机构内部储层物性及含气性预测, 对其它类似火山岩储层的地区具有推广意义。