CN104570065A - 一种利用地震波阻抗定量反演孔隙度的方法 - Google Patents

Abstract

本发明提供了一种利用地震波阻抗定量反演孔隙度的方法，属于石油与天然气勘探开发领域。本方法以地震波阻抗为输入数据，利用岩石物理方程和地震波传播速度公式构建孔隙度反演的误差函数，并在孔隙度的取值范围内构建解的可行域空间，在该可行域空间内采用非线性全局寻优算法进行迭代反演，得到最终的孔隙度。本发明基于岩石物理方程和弹性参数进行地震孔隙度反演，物理意义明确，适应性强，将孔隙度的取值范围加入反演过程作为约束，提高了反演的稳定性和精度，可保持叠前弹性参数反演的优点，同时又具有对油气分布有利区带直接预测的优势，可直接用于实际地震资料的孔隙度反演。

Description

一种利用地震波阻抗定量反演孔隙度的方法

技术领域

本发明属于石油与天然气勘探开发领域，具体涉及一种利用地震波阻抗定量反演孔隙度的方法。

背景技术

地震勘探技术是油气勘探中应用最为广泛的一种地球物理学方法，在地表或者井中人工激发地震波，利用地震波在不同介质中传播的速度、振幅、频率、相位、波形等参数的变化来分析、预测油气储层分布范围及储层物性特征。随勘探开发程度逐渐深入，地震勘探已经从认识地下构造形态的构造勘探，逐渐发展成直接应用地震信息判断岩性、分析岩相、定量计算岩层物性参数的岩性勘探。地震性质及反射率受控于地下介质岩性、流体类型、孔隙度、压力、温度、矿物、岩理等，而岩石物理正是研究岩石物理属性(物性参数)与地球物理属性之间的桥梁，即岩石物理属性可由地球物理观测数据及其特征(属性)来表征，岩石物理是建立地震资料预测工具及解释反演结果的物理基础。近几年来，以Gassmann方程为代表的岩石物理理论研究在许多方面都取得了很大的进展。利用岩石物理技术研究的成果，可以为含油气分析和储层识别提供敏感性岩石物理参数，能够有效地指导储层识别和预测。孔隙度和饱和度是描述储层性质的重要参数，利用地震资料定量反演储层参数具有横向连续性好的优点，对储层精细描述具有重要意义。

目前地震孔隙度反演方法主要有三类：

①利用速度求孔隙度，根据孔隙度与地震波传播速度的关系求取孔隙度，例如利用Wyllie公式和改进型的Raymer公式由地震纵波速度等参数求取孔隙度。该类方法的优点是计算方便，但Wyllie公式未考虑岩石的体积形变，难以在复杂地质条件下应用。后来的研究表明，它们也不适用于含气砂岩孔隙度的预测；

②统计型反演方法，由于地震声学反演和弹性反演的结果与岩石孔隙度无直接联系，因此在将反演结果转换成孔隙度时，主要利用了反演结果和孔隙度的统计关系，例如速度-孔隙度交汇图以及密度、波阻抗、泊松比等参数和孔隙度的交汇图。该类方法目前应用较普遍，但要获得较好的效果需要大量的地质、测井、岩石物理、地震多属性参数资料，并构建适当的地质模型；

③第三类方法是基于Biot-Gassmann方程的孔隙度计算方法，如HeXilei与He Zhenhua(2012)提出的岩石骨架与地震孔隙度反演。张应波(1994)推导了孔隙度的解析式，该方法的优点是获得了计算孔隙度的解析式，有坚实的理论基础，缺点是需要预先提供的参数很多，除了Gassmann方程中孔隙度以外的4个参数外，还要求提供应力、孔隙压力、流体黏滞系数和地震波的衰减系数等。上述参数的获得是一个复杂的过程，使得基于地震数据的地震孔隙度反演的应用受到了限制。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种利用地震波阻抗定量反演孔隙度的方法，针对地震资料解释中难以获得孔隙度剖面及常规地震孔隙度反演方法存在的反演精度低、输入数据繁杂及难以适用复杂地质条件，应用受到限制等不足，本发明基于岩石物理方程和弹性参数进行地震孔隙度反演，物理意义明确，适应性强，将孔隙度的取值范围加入反演过程作为约束，提高了反演的稳定性和精度，可保持叠前弹性参数反演的优点，同时又具有对油气分布有利区带直接预测的优势，可直接用于实际地震资料的孔隙度反演。

本发明是通过以下技术方案实现的：

一种利用地震波阻抗定量反演孔隙度的方法，以地震波阻抗为输入数据，利用岩石物理方程和地震波传播速度公式构建孔隙度反演的误差函数，并在孔隙度的取值范围内构建解的可行域空间，在该可行域空间内采用非线性全局寻优算法进行迭代反演，得到最终的孔隙度。

所述方法包括：

(1)输入终止条件的值ε，在孔隙度的取值范围内随机给定孔隙度的初始值；

(2)输入固体矿物的体积模量K₀和剪切模量μ₀，输入固结系数c，计算干岩石弹性模量；

(3)计算饱和岩石弹性模量；

(4)由所述饱和岩石弹性模量计算波阻抗，所述波阻抗包括纵波阻抗I_p和横波阻抗I_s；

(5)输入比重因子λ和由叠前地震反演得到的纵波阻抗I_pi和横波阻抗I_si’建立孔隙度反演的误差函数；

(6)利用非线性全局寻优算法对所述孔隙度反演的误差函数进行迭代反演，直到所述孔隙度反演的误差函数满足终止条件；

(7)输出满足终止条件的孔隙度。

所述步骤(2)计算干岩石弹性模量是这样实现的：

所述干岩石弹性模量包括体积模量和剪切模量，分别由公式(1)和(2)计算得到：

K_dry=K₀(1-φ)／(1+cφ)  (1)

μ_dry=μ₀(1-φ)／(1+1.5cφ)  (2)

其中，K_dry和μ_dry分别为干岩石体积模量和剪切模量，φ为孔隙度。

所述步骤(3)是这样实现的：

利用岩石物理方程计算饱和岩石弹性模量，所述饱和岩石弹性模量包括体积模量和剪切模量，所述岩石物理方程为公式(3)和公式(4)：

$K_{\mathrm{sat}}=K_{\mathrm{dry}}+{(1-K_{\mathrm{dry}}/K_0)}^2/[\mathrm{\φ}/K_f+(1-\mathrm{\φ})/K_0-K_{\mathrm{dry}}/K_0^2]---\left(3\right)$

μ_sat=μ_dry  (4)

其中，Kf为孔隙饱含流体的体积模量，K_sat和μ_sat为饱和岩石体积模量和剪切模量。

所述步骤(4)是这样实现的：

由地震波在均匀、各向同性、弹性介质中的传播速度，利用地震波传播速度公式得到纵波阻抗和横波阻抗，所述地震波传播速度公式为公式(5)和公式(6)：

$I_p=\sqrt{{\mathrm{\ρ}}_{\mathrm{sat}}(K_{\mathrm{sat}}+4{\mathrm{\μ}}_{\mathrm{sat}}/3)}---\left(5\right)$

$I_s=\sqrt{{\mathrm{\ρ}}_{\mathrm{sat}}{\mathrm{\μ}}_{\mathrm{sat}}}---\left(6\right)$

其中，ρ_sat为饱和岩石密度。

所述步骤(5)中建立的孔隙度反演的误差函数如下：

f(φ)=λ[I_p-I_pi]²+(1-λ)[I_s-I_si]²  (7)

其中，0≤λ≤1，由用户给定。

所述步骤(6)包括：

(61)首先由孔隙度的取值范围建立可行域空间，如公式(8)：

φ∈[φ_min，φ_max]  (8)

其中，φ_min为孔隙度取值范围的最小值，φ_max为孔隙度取值范围的最大值；

(62)在所述可行域空间内生成一个孔隙度可行解，并重复步骤(2)-(5)得到误差函数的值，然后判断该误差函数的值是否满足终止条件，如果否，则返回步骤(62)，如果是，则停止。

在所述可行域空间内生成一个孔隙度可行解是利用下式实现的：

φ_r=r(φ_max-φ_min)  (10)

其中，φ_r为孔隙度可行解，r为随机数，其范围是0≤r≤1。

所述终止条件如下：

f(φ)≤ε  (9)。

与现有技术相比，本发明的有益效果是：基于地震资料的孔隙度反演在油气／流体识别方面比弹性三参数(纵、横波速度和岩石密度)反演更直接、有效。本发明基于岩石物理方程和弹性参数进行地震孔隙度反演，物理意义明确，适应性强，将孔隙度的取值范围加入反演过程作为约束，提高了反演的稳定性和精度，可保持叠前弹性参数反演的优点，同时又具有对油气分布有利区带直接预测的优势，可直接用于实际地震资料的孔隙度反演。

附图说明

图1是干岩石弹性模量随孔隙度变化曲线，上部为体积模量，下部为剪切模量，固结系数c取值为2。

图2是干岩石弹性模量随孔隙度变化曲线，上部为体积模量，下部为剪切模量，固结系数c取值为5。

图3是干岩石弹性模量随孔隙度变化曲线，上部为体积模量，下部为剪切模量，固结系数c取值为10。

图4是饱和岩石弹性模量随孔隙度变化曲线，上部为体积模量，下部为剪切模量，固结系数c取值为2。

图5是饱和岩石弹性模量随孔隙度变化曲线，上部为体积模量，下部为剪切模量，固结系数c取值为5。

图6是饱和岩石弹性模量随孔隙度变化曲线，上部为体积模量，下部为剪切模量，固结系数c取值为10。

图7是波阻抗随孔隙度变化曲线，从上至下，固结系数c取值为2、5、10。

图8是波阻抗剖面，在本发明中该数据作为输入数据给定。

图9是本发明方法反演的孔隙度剖面，与测井数据温和良好。

图10是本发明方法的步骤框图。

具体实施方式

下面结合附图对本发明作进一步详细描述：

孔隙度是描述储层性质的重要参数，因为流体是直接赋存在岩石孔隙中的，所以基于地震资料的孔隙度反演在油气／流体识别方面有时比弹性三参数(纵、横波速度和岩石密度)反演更直接、有效。

本发明针对地震资料解释中难以获得孔隙度剖面及常规地震孔隙度反演方法存在的反演精度低、输入数据繁杂及难以适用复杂地质条件的不足，提供了一种利用地震波阻抗定量反演孔隙度的方法，以地震波阻抗为输入数据，利用岩石物理方程和地震波传播速度公式构建了地震孔隙度反演的误差函数，并在孔隙度取值范围的基础上构建解的可行域空间，在该可行域空间内采用非线性全局寻优算法进行迭代反演，将孔隙度的取值范围加入反演过程作为约束，提高了反演的稳定性和精度。

本发明是通过以下技术方案实现的：

一种利用地震波阻抗定量反演孔隙度的方法，首先利用Biot-Gassmann岩石物理方程和地震波在地下介质中的传播速度公式建立了地震孔隙度反演的误差函数，然后在孔隙度取值范围的基础上构建解的可行域空间，最后在该可行域空间内采用非线性全局寻优算法进行迭代反演。

如图10所示，所述方法包括以下步骤：

(1)计算干岩石的弹性模量，包括体积模量和剪切模量，分别由公式(1)和(2)计算得到：

K_dry=K₀(1-φ)／(1+cφ)  (1)

μ_dry=μ₀(1-φ)／(1+1.5cφ)  (2)

上面两个公式的右边，除开φ，其它都是已知的数据。

其中，K₀和μ₀分别为固体矿物的体积模量和剪切模量(已知的，是从岩石物理实验室得到。)，K_dry和μ_dry分别为干岩石体积模量和剪切模量(这两个参数是在给定固体矿物的体积模量和剪切模量K₀和μ₀及孔隙度φ后，由公式(1)和(2)计算得到的)，φ为孔隙度，c是固结系数(该参数对于砂岩有个取值范围，在该取值范围内由用户任意取值)，表示岩石的固结程度，对于砂岩，通常取2≤c≤20。

图1-图3是利用步骤(1)计算的干岩石弹性模量随孔隙度变化情况，固结系数依次取为2、5、10。

(2)由Biot-Gassmann岩石物理方程计算饱和岩石的弹性模量，包括体积模量和剪切模量，分别由公式(3)和(4)计算得到：

$K_{\mathrm{sat}}=K_{\mathrm{dry}}+{(1-K_{\mathrm{dry}}/K_0)}^2/[\mathrm{\φ}/K_f+(1-\mathrm{\φ})/K_0-K_{\mathrm{dry}}/K_0^2]---\left(3\right)$

μ_sat=μ_dry  (4)

其中，Kf为孔隙饱含流体的体积模量(岩石孔隙中包括流体有油、气、水等，这个参数取为水的体积模量，是已知的)，K_sat和μ_sat为饱和岩石体积模量和剪切模量。

图4-6是利用步骤(2)计算的饱和岩石弹性模量随孔隙度变化情况，固结系数依次取为2、5、10。

(3)由地震波在均匀、各向同性、弹性介质中的传播速度得到纵波阻抗和横波阻抗：

$I_p=\sqrt{{\mathrm{\ρ}}_{\mathrm{sat}}(K_{\mathrm{sat}}+4{\mathrm{\μ}}_{\mathrm{sat}}/3)}---\left(5\right)$

$I_s=\sqrt{{\mathrm{\ρ}}_{\mathrm{sat}}{\mathrm{\μ}}_{\mathrm{sat}}}---\left(6\right)$

其中I_p为纵波阻抗，I_s为横波阻抗，ρ_sat为饱和岩石密度(这个值由测井数据得到)。

图7是利用步骤(3)计算的波阻抗随孔隙度变化情况，固结系数依次取为2、5、10。

(4)由上述步骤(3)计算得到的纵波阻抗和横波阻抗建立地震孔隙度反演的误差函数，由公式(7)计算得到：

f(φ)=λ[I_p-I_pi]²+(1-λ)[I_s-I_si]²  (7)

所述公式(7)中，I_p和I_s为分别为由上述步骤(3)计算纵波阻抗和横波阻抗，I_pi和I_si为为由叠前地震反演得到的纵波阻抗和横波阻抗，这里作为已知数据给定。λ为调节纵波阻抗和横波阻抗在误差函数中的比重因子，其取值范围为0≤λ≤1，当λ=0时，只有横波阻抗参与计算；当λ=1时，只有纵波阻抗参与计算，λ默认取值为0.5(这个值由用户在取值范围内任意取值。λ为调节纵波阻抗和横波阻抗在误差函数中的比重。λ越大，就意味着纵波阻抗在误差函数中的比重大，当λ取为1时，那么就只有纵波阻抗参与反演；当λ取为0时，只有横波阻抗参与反演。λ取为0.5，是两者在误差函数中占有同等的比重。)。

(5)由孔隙度的取值范围建立可行域空间，由公式(8)计算：

φ∈[φ_min，φ_max]  (8)

所述公式(8)中φ_min为孔隙度取值范围的最小值，φ_max为孔隙度取值范围的最大值，φ_min和φ_max由测井数据得到(这两个值是在孔隙度测井曲线中选出来的。孔隙度测井曲线由常规测井得到，在孔隙度测井曲线中选出最大值和最小值，选出的最小值和最大值就是φ_min和φ_max)。

(6)在上述步骤(5)中的可行域空间中利用非线性全局寻优算法对上述步骤(4)中的误差函数进行迭代反演：在孔隙度可行域空间内产生孔隙度可行解，继而利用步骤(1)-(4)计算岩石弹性模量、纵横波阻抗及误差函数值，重复上述过程，直到产生的孔隙度可行解使得误差函数值满足终止条件为止。孔隙度可行解生成方式和终止条件分别由下式(9)和(10)给出：

f(φ)≤ε  (9)

φ_r=r(φ_max-φ_min)  (10)

所述公式(9)和(10)中，ε为较小的数(由用户在最开始的时候给定的)，如ε=0.001。φ_r为孔隙度可行解，r为随机数，其范围是0≤r≤1。

本发明具有如下特点：

由Biot-Gassmann岩石物理方程和地震波的传播速度公式推导出用于地震孔隙度反演的误差函数，该误差函数中仅需要纵、横波阻抗作为输入数据，不需要其他解释成果数据。纵横波阻抗可由叠前地震反演得到。该误差函数只反演孔隙度一个参数，避免了超高维问题所固有的多解性问题，并且提高了反演的效果。

本发明以孔隙度的取值范围为基础建立了孔隙度可行域空间，在该可行域空间内进行反演，将孔隙度的取值范围作为约束加入到反演的过程中，提高了反演的稳定性和精度。采用非线性全局寻优算法在可行域空间内进行反演。

本发明中，流体直接赋存在岩石孔隙中的，基于地震资料的孔隙度反演在油气／流体识别方面比弹性参数(纵、横波速度和岩石密度)反演更直接、有效，物理意义明确，适应性强。

下面结合附图及实际数据反演为例，说明本发明的效果。

由于流体直接赋存在基于地震资料的孔隙度反演在油气／流体识别方面比弹性三参数(纵、横波速度和岩石密度)反演更直接、有效。针对地震资料解释中难以获得孔隙度剖面及常规地震孔隙度反演方法存在的反演精度低、输入数据繁杂及难以适用复杂地质条件，应用受到限制等不足，本发明基于岩石物理方程和弹性参数进行地震孔隙度反演，物理意义明确，适应性强，将孔隙度的取值范围加入反演过程作为约束，提高了反演的稳定性和精度，可保持叠前弹性参数反演的优点，同时又具有对油气分布有利区带直接预测的优势。

下面以某地区的实际地震数据为例，利用本发明方法进行孔隙度反演，进而说明本发明的效果。

实施例1

本实施例是本发明用于某地区的实际地震数据的孔隙度反演实施实例，一种利用地震波阻抗定量反演孔隙度的方法的步骤如下：

①给定算法终止条件，在孔隙度的取值范围内随机给定孔隙度的初始值；

②给定固体矿物的体积模量K₀和剪切模量μ₀，给定固结系数，利用本发明步骤(1)计算干岩石弹性模量；

③利用本发明步骤(2)计算饱和岩石弹性模量；

④由上述步骤计算的饱和岩石弹性模量和本发明步骤(3)计算波阻抗，包括纵波阻抗和横波阻抗；

⑤给定比重因子和波阻抗剖面，由本发明步骤(4)建立孔隙度反演的误差函数；

⑥利用非线性全局寻优算法对上述步骤的误差函数进行迭代反演，重复上述步骤①-⑤，直到误差函数满足算法终止条件；

⑦输出满足算法终止条件的孔隙度。

图8是本实施例中实际地震数据反演得到的波阻抗，在本发明中波阻抗作为输入数据给定。图9是利用本发明方法反演得到的孔隙度剖面，可以看出，其与测井数据吻合良好。

本发明可用于石油地球物理勘探中的储层预测及评价技术领域。本发明从地震资料中直接反演孔隙度剖面，覆盖面广，并有较高的横向分辨率。本发明方法以地震波阻抗数据为基础，首先利用Biot-Gassmann岩石物理方程和地震波在地下介质中的传播速度公式建立了地震孔隙度反演的误差函数，然后基于孔隙度的取值范围构建了解的可行域空间，最后在该可行域空间内采用非线性全局寻优算法进行迭代反演。本发明基于岩石物理方程和弹性参数进行地震孔隙度反演，物理意义明确，适应性强，将孔隙度的取值范围加入反演过程作为约束，提高了反演的稳定性和精度，可保持叠前弹性参数反演的优点，同时又具有对油气分布有利区带直接预测的优势，可直接用于实际地震资料的孔隙度反演。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (9)

1.一种利用地震波阻抗定量反演孔隙度的方法，其特征在于：所述方法以地震波阻抗为输入数据，利用岩石物理方程和地震波传播速度公式构建孔隙度反演的误差函数，并在孔隙度的取值范围内构建解的可行域空间，在该可行域空间内采用非线性全局寻优算法进行迭代反演，得到最终的孔隙度。

2.根据权利要求1所述的利用地震波阻抗定量反演孔隙度的方法，其特征在于：所述方法包括：

(1)输入终止条件的值ε，在孔隙度的取值范围内随机给定孔隙度的初始值；

(2)输入固体矿物的体积模量K₀和剪切模量μ₀，输入固结系数c，计算干岩石弹性模量；

(3)计算饱和岩石弹性模量；

(4)由所述饱和岩石弹性模量计算波阻抗，所述波阻抗包括纵波阻抗I_p和横波阻抗I_s；

(5)输入比重因子λ和由叠前地震反演得到的纵波阻抗I_pi和横波阻抗I_si’建立孔隙度反演的误差函数；

(6)利用非线性全局寻优算法对所述孔隙度反演的误差函数进行迭代反演，直到所述孔隙度反演的误差函数满足终止条件；

(7)输出满足终止条件的孔隙度。

3.根据权利要求2所述的利用地震波阻抗定量反演孔隙度的方法，其特征在于：所述步骤(2)计算干岩石弹性模量是这样实现的：

所述干岩石弹性模量包括体积模量和剪切模量，分别由公式(1)和(2)计算得到：

K_dry=K₀(1-φ)／(1+cφ)  (1)

μ_dry=μ₀(1-φ)／(1+1.5cφ)  (2)

其中，K_dry和μ_dry分别为干岩石体积模量和剪切模量，φ为孔隙度。

4.根据权利要求3所述的利用地震波阻抗定量反演孔隙度的方法，其特征在于：所述步骤(3)是这样实现的：

利用岩石物理方程计算饱和岩石弹性模量，所述饱和岩石弹性模量包括体积模量和剪切模量，所述岩石物理方程为公式(3)和公式(4)：

$K_{\mathrm{sat}}=K_{\mathrm{dry}}+{(1-K_{\mathrm{dry}}/K_0)}^2/[\mathrm{\φ}/K_f+(1-\mathrm{\φ})/K_0-K_{\mathrm{dry}}/K_0^2]---\left(3\right)$

μ_sat=μ_dry  (4)

其中，K_f为孔隙饱含流体的体积模量，K_sat和μ_sat为饱和岩石体积模量和剪切模量。

5.根据权利要求4所述的利用地震波阻抗定量反演孔隙度的方法，其特征在于：所述步骤(4)是这样实现的：

由地震波在均匀、各向同性、弹性介质中的传播速度，利用地震波传播速度公式得到纵波阻抗和横波阻抗，所述地震波传播速度公式为公式(5)和公式(6)：

$I_p=\sqrt{{\mathrm{\ρ}}_{\mathrm{sat}}(K_{\mathrm{sat}}+4{\mathrm{\μ}}_{\mathrm{sat}}/3)}---\left(5\right)$

$I_s=\sqrt{{\mathrm{\ρ}}_{\mathrm{sat}}{\mathrm{\μ}}_{\mathrm{sat}}}---\left(6\right)$

其中，ρ_sat为饱和岩石密度。

6.根据权利要求5所述的利用地震波阻抗定量反演孔隙度的方法，其特征在于：所述步骤(5)中建立的孔隙度反演的误差函数如下：

f(φ)=λ[I_p-I_pi]²+(1-λ)[I_s-I_si]²  (7)

其中，0≤λ≤1，由用户给定。

7.根据权利要求6所述的利用地震波阻抗定量反演孔隙度的方法，其特征在于：所述步骤(6)包括：

(61)首先由孔隙度的取值范围建立可行域空间，如公式(8)：

φ∈[φ_min，φ_max]  (8)

其中，φ_min为孔隙度取值范围的最小值，φ_max为孔隙度取值范围的最大值；

(62)在所述可行域空间内生成一个孔隙度可行解，并重复步骤(2)-(5)得到误差函数的值，然后判断该误差函数的值是否满足终止条件，如果否，则返回步骤(62)，如果是，则停止。

8.根据权利要求7所述的利用地震波阻抗定量反演孔隙度的方法，其特征在于：在所述可行域空间内生成一个孔隙度可行解是利用下式实现的：

φ_r=r(φ_max-φ_min)  (10)

其中，φ_r为孔隙度可行解，r为随机数，其范围是0≤r≤1。

9.根据权利要求8所述的利用地震波阻抗定量反演孔隙度的方法，其特征在于：所述终止条件如下：

f(φ)≤ε  (9)。