CN104570064A - 一种砂岩地层横波速度计算方法 - Google Patents

Abstract

本发明提供了一种砂岩地层横波速度计算方法，属于岩石物理学及地球物理测井领域。本方法包括：(1)输入砂岩地层的实测纵波速度V_p、密度、孔隙度φ和泥质含量V_sh曲线；给定固结系数c的初始值和间隔；(2)计算干岩石弹性模量；(3)计算饱和岩石弹性模量，并计算得到计算纵波速度；(4)建立实测纵波速度V_p与所述计算纵波速度的误差函数；(5)在固结系数c的取值范围内以所述间隔改变固结系数c的值，每改变一次，重复一次步骤(2)至步骤(4)，得到每个新固结系数c对应的误差函数，然后比较所有固结系数c对应的误差函数，找到使得误差函数最小的固结系数c′，利用该使得误差函数最小的固结系数c′计算得到的横波速度即为所求的横波速度。

Description

一种砂岩地层横波速度计算方法

技术领域

本发明属于岩石物理学及地球物理测井领域，具体涉及一种砂岩地层横波速度计算方法。

背景技术

油气勘探目标已由简单的以构造油气藏勘探为主转入以复杂的地层、岩性油气藏或构造-地层岩性复合油气藏勘探为主，勘探难度大大增加。油气勘探的地区范围已由东部陆地平原为主转入西部山区、沙漠、黄土区和近海海域，在这些地区获取合格的地质、地球物理资料的难度大、成本高。油气勘探目的层的埋深愈来愈大，有的超过了7000m，属于深层或超深层勘探，使地球物理勘探的一些有效的技术难以应用。目前为降低复杂地区和深层的勘探成本，避免钻探到干井和水井的风险，油气工业界希望尽可能正确区分和预测储层中含流体(油、气、水)的状况。目前采用的叠前AVO分析技术、叠前反演技术在流体预测中起到了重要作用。这些叠前地震资料AVO分析及叠前反演等技术均需要准确的纵波、横波速度测井曲线，然而在实际测井资料中往往缺乏横波速度信息，比如开发井中只有极少数井有横波测井资料，而开发时间较长的老油田则横波测井资料更少，这给储层预测工作带来了一定的困难。为此，许多地球物理工作者对横波速度预测进行了大量的方法研究工作，其中有统计经验公式法和理论模型方法。如Han(1986)、K1imentors(1991)、Greenberg(1992)、Castagna(1993)、Goldberg(1998)给出的岩石速度、孔隙度及岩石其它参数之间的经验公式，其缺陷是受地域性影响大，精度较低。目前在老油田应用井震联合技术进行储层精细描述，剩余油分布预测以及利用地震数据反演储层参数时，都需要有纵、横波测井曲线和密度测井曲线来建立地层模型。此外在地震岩石属性反演研究中也需要使用横波速度测井资料。但是应用上述经验公式计算的横波速度误差较大，会导致错误的解释结果。比如人们利用简化的Zoepprits方程来计算AVO响应(Aki and Lichards，1980；Shuey，1985)，而在这些简化的方程中所用到的横波速度又是通过经验公式或者纵横波速度比为常数的方法获得。这些方式所计算出的横波速度往往与实际地层的横波速度不吻合，从而导致计算出的AVO有较大的误差。由于常规统计经验公式法存在的受地域性影响大，精度较低等问题，使得横波速度的应用受到限制。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种砂岩地层横波速度计算方法，针对常规横波速度计算方法存在的受地域性影响大、计算误差大及运行效率低，不适于测井资料的实时解释及叠前地震反演和叠前地震属性分析等工作，应用受到限制等缺点，本发明方法利用岩石物理学原理从测井数据中计算得到横波速度。本发明方法所需要的参数少，只需要纵波速度、孔隙度和泥质含量，就可以预测出比较准确的横波速度，可以避免传统的经验公式中因统计关系而存在的不确定性，本发明应用范围广。具有计算速度快、稳定性好的优点，可直接用于测井数据解释及其他储层预测工作。

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

一种砂岩地层横波速度计算方法，包括：

(1)输入砂岩地层的实测纵波速度V_p、密度、孔隙度φ和泥质含量V_sh曲线；给定固结系数c的初始值和间隔；

(2)计算干岩石弹性模量；

(3)计算饱和岩石弹性模量，并计算得到计算纵波速度；

(4)建立实测纵波速度V_p与所述计算纵波速度的误差函数；

(5)在固结系数c的取值范围内以所述间隔改变固结系数c的值，每改变一次，重复一次步骤(2)至步骤(4)，得到每个新固结系数c对应的误差函数，然后比较所有固结系数c对应的误差函数，找到使得误差函数最小的固结系数c′，利用该使得误差函数最小的固结系数c′计算得到的横波速度即为所求的横波速度。

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

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

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

其中，K₀和μ₀分别为固体矿物的体积模量和剪切模量，K_dry和μ_dry分别为干岩石体积模量和剪切模量；

所述公式(1)中，固体矿物的体积模量K₀由下述公式(3)计算得到：

K₀=[V_shK_sh+V_sK_s+1/(V_sh/K_sh+V_s/K_s)]/2     (3)

式中，K₀、K_s、K_sh分别代表固体矿物的体积模量、砂岩的体积模量、泥岩的体积模量；

V_sh和V_s分别为泥岩和砂岩的体积百分含量，由测井数据的泥质含量得到；

砂岩的体积模量K_s和泥岩的体积模量K_sh由岩石物理实验得到；

所述公式(1)中，固体矿物的剪切模量μ0由下述公式(4)计算得到

μ₀=[V_shμ_sh+V_sμ_s+1/(V_sh/μ_sh+V_s/μ_s)]/2      (4)

式中，μ₀、μ_s、μ_sh分别代表固体矿物的剪切模量、砂岩的剪切模量、泥岩的剪切模量，其中砂岩的剪切模量μ_s和泥岩的剪切模量μ_sh由岩石物理实验得到。

所述步骤(3)中的所述饱和岩石弹性模量包括体积模量和剪切模量，由公式(5)计算得到：

$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],{\mathrm{\μ}}_{\mathrm{sat}}={\mathrm{\μ}}_{\mathrm{dry}}---\left(5\right)$

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

所述步骤(3)中的计算纵波速度由公式(6)计算得到：

$V_p^{\mathrm{est}}=\sqrt{(K_{\mathrm{ast}}+{4\mathrm{\μ}}_{\mathrm{sat}}/3)/{\mathrm{\ρ}}_{\mathrm{sat}}},V_s^{\mathrm{est}}=\sqrt{{\mathrm{\μ}}_{\mathrm{sat}}/{\mathrm{\ρ}}_{\mathrm{sat}}}---\left(6\right)$

其中，为纵波速度，为横波速度，ρ_sat为饱含流体的岩石密度。

所述步骤(4)中的误差函数由公式(7)确定：

$f\left(c\right)={\left|\right|V_p-V_p^{\mathrm{est}}\left|\right|}_2---\left(7\right).$

所述步骤(5)中，对于砂岩，固结系数c的取值范围取2≤c≤20。

所述步骤(5)中利用该使得误差函数最小的固结系数c′计算得到的横波速度即为所求的横波速度是这样实现的：

利用所述使得误差函数最小的固结系数c′，重复步骤(2)和步骤(3)得到横波速度

与现有技术相比，本发明的有益效果是：本发明利用常规测井数据计算砂岩地层横波速度，所需要的参数少，只需要纵波速度、孔隙度和泥质含量，就可以预测出比较准确的横波速度，避免了传统的经验公式中因统计关系而存在的不确定性，应用范围广。具有计算速度快、稳定性好的优点，可直接用于测井数据解释工作。

附图说明

图1是实施例1中的纵波速度、孔隙度、泥质含量、密度测井曲线。

图2是实施例1中的纵波速度、横波速度、孔隙度和泥质含量。从左至右依次为实测纵波速度与预测的纵波速度的拟合情况、实测横波速度与预测横波速度的拟合情况、孔隙度、泥质含量。

图3是实施例1中的实测横波速度与预测横波速度交会图。

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

具体实施方式

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

本发明的目的就是针对现有横波速度估算方法在测井资料应用中存在的不足，围绕现有横波速度估算方法存在的计算误差大、运行效率低、适用范围受地域限制等问题，提供一种砂岩地层横波速度计算方法，本发明方法以常规测井资料为基础，所需要的参数少，只需要纵波速度、孔隙度和泥质含量，就可以预测出比较准确的横波速度，可以避免传统的经验公式中因统计关系而存在的不确定性，应用范围广，能够为叠前反演等储层预测工作提供基础资料。

为了从常规测井数据中计算横波速度，本发明方法基于常规测井资料，首先利用Pride公式和Biot-Gassmann方程建立了砂岩地层横波速度预测的误差函数，然后在固结系数取值范围内以一维搜索的方式来调整固结系数，直到误差函数最小为止，利用使得误差函数取得最小值的固结系数计算横波速度。

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

一种砂岩地层横波速度计算方法，首先利用Pride公式和Biot-Gassmann方程建立了砂岩地层横波速度预测的误差函数，然后在固结系数取值范围内以一维搜索的方式来调整固结系数，得到一个使得误差函数最小的固结系数，最后利用该固结系数计算横波速度。

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

(1)输入纵波速度V_p、孔隙度φ和泥质含量V_sh曲线。以上所述的曲线由测井数据得到。

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

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

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

其中，K₀和μ₀分别为固体矿物的体积模量和剪切模量，K_dry和μ_dry分别为干岩石体积模量和剪切模量，φ为孔隙度，c是固结系数，表不岩石的固结程度，对于砂岩，通常取2≤c≤20。

所述公式(1)中，固体矿物的体积模量K₀由下述公式(3)计算得到：

K₀=[V_shK_sh+V_sK_s+1/(V_sh/K_sh+V_s/K_s)]/2      (3)

式中，K₀、K_s、K_sh分别代表固体矿物的体积模量、砂岩的体积模量、泥岩的体积模量。V_sh和V_s分别为泥岩和砂岩的体积百分含量，由测井数据的泥质含量得到。砂岩的体积模量K_s和泥岩的体积模量K_sh为输入数据，由岩石物理实验得到，见表1。

所述公式(1)中，固体矿物的剪切模量μ₀计算方式与体积模量K₀计算方式相同。

μ₀=[V_shμ_sh+V_sμ_s+1/(V_sh/μ_sh+V_s/μ_s)]/2        (4)

式中，μ₀、μ_s、μ_sh分别代表固体矿物的剪切模量、砂岩的剪切模量、泥岩的剪切模量，其中砂岩的剪切模量μ_s和泥岩的剪切模量μ_sh由为输入数据，由岩石物理实验得到，见表1。

(3)由Biot-Gassmann岩石物理方程计算饱和岩石的弹性模量，并计算地层的速度。包括体积模量和剪切模量、纵波速度和横波速度，分别由公式(5)和(6)计算得到：

$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],{\mathrm{\μ}}_{\mathrm{sat}}={\mathrm{\μ}}_{\mathrm{dry}}---\left(5\right)$

$V_p^{\mathrm{est}}=\sqrt{(K_{\mathrm{ast}}+{4\mathrm{\μ}}_{\mathrm{sat}}/3)/{\mathrm{\ρ}}_{\mathrm{sat}}},V_s^{\mathrm{est}}=\sqrt{{\mathrm{\μ}}_{\mathrm{sat}}/{\mathrm{\ρ}}_{\mathrm{sat}}}---\left(6\right)$

其中，K_f为孔隙饱含流体的体积模量(孔隙流体体积模量K_f为输入参数，见表1中水的体积模量)，K_sat和μ_sat为饱含流体岩石体积模量和剪切模量。

其中，为纵波速度，为横波速度，ρ_sat为饱含流体的岩石密度(由测井数据的密度曲线得到)。

(4)建立砂岩地层实测纵波速度与上述步骤理论计算纵波速度的误差函数，由公式(7)确定：

$f\left(c\right)={\left|\right|V_p-V_p^{\mathrm{est}}\left|\right|}_2---\left(7\right).$

其中V_p为实测纵波速度，即测井数据中的纵波速度(就是步骤(1)中输入的)。为理论计算的纵波速度。c是步骤(2)中的固结系数。

(5)在固结系数c的取值范围内以一定的间隔依次计算上述公式(6)中的误差函数，取得使得误差函数最小的那个固结系数c′。这个过程需要重复上述步骤(2)-(4)。然后利用该固结系数c′计算横波速度。

所述的间隔可以取0.1，也可以取其他值，如0.02，0.05，0.1，0.15，0.2，0.25，0.3等等。

下面以某地区的实际测井数据为例，利用本发明方法进行横波速度的计算，进而说明本发明的效果。

实施例1本实施例是本发明用于某地区的实际测井数据进行横波速度计算的实施实例，一种快速求取地层横波速度的计算方法的步骤如下：

(1)输入纵波速度V_p、密度、孔隙度φ和泥质含量V_sh曲线；

(2)根据表1中的弹性参数及上述步骤输入的泥质含量，由公式(3)计算固体矿物的弹性模量K₀和μ₀；

(3)给定固结系数c初始值及间隔(一般取0.1)，由上述步骤(2)计算的固体矿物的弹性模量和公式(1)及(2)计算干岩石弹性模量K_dry和μ_dry；

(4)利用上述步骤(3)中计算的干岩石弹性模量K_dry和μ_dry，由Gassmann公式计算流体饱和岩石的弹性模量K_sat和μ_sat；

(5)由公式(5)和步骤(1)输入的密度计算饱和流体岩石的纵波速度，即理论计算的纵波速度

(6)根据公式(6)建立理论计算的纵波速度与实测纵波速度V_p的误差函数；

(7)在固结系数的取值范围内，以一维搜索的方式求取使得误差函数取得最小值时的固结系数。按照间隔(一般取0.1)依次计算固结系数在各个取值情况下误差函数值，将下一次取值情况的误差函数值与当前的误差函数值进行比较，保留截止当前使得目标函数最小的固结系数，循环进行，直到计算完成固结系数的取值范围。在这个过程需要对上述步骤(3)-(6)进行迭代。

具体实施时，可以采用如图4所示的f(c)≤ε.其中ε是大于0的一个实数。目标函数f(c)是平方和后再开根号，因此f(c)的理论最小值是0，但是通常不需要这么小的值就可以，比如f(c)的值比0.01小了，这个时候就可以停止算法迭代，不需要再去与固结系数取其他值时的目标函数进行比较了。输出满足f(c)≤ε时的横波速度。

(8)将步骤(7)中求取的固结系数代入步骤(3)-(4)，再利用公式(5)计算得到横波速度。

图1为某地区的实际测井数据，包括纵波速度、孔隙度、泥质含量、密度测井曲线。

表1为常用岩石的弹性参数，包括砂岩、泥岩及水的弹性模量。

图2为本发明方法计算的纵波速度和横波速度与实测纵波速度和横波速度的拟合情况。预测的横波速度与实测横波速度相对误差均值为10％左右。

图3为本发明方法计算的横波速度与实测横波速度的交会图，标准偏差为68％。

岩性及弹性参数	数值
		砂岩的体积模量	38GPa
砂岩的剪切模量	44GPa
		泥岩的体积模量	20.9GPa
泥岩的剪切模量	6.85GPa
		水的体积模量	2.29GPa

表1

本发明公开一种砂岩地层横波速度计算方法，属于岩石物理学及地球物理测井技术领域。本发明方法基于常规测井资料，首先利用Pride公式和Biot-Gassmann方程建立了砂岩地层横波速度预测的误差函数，然后在固结系数取值范围内以一维搜索的方式来调整固结系数，直到误差函数最小为止。本发明方法所需要的参数少，只需要纵波速度、孔隙度和泥质含量，就可以预测出比较准确的横波速度，可以避免传统的经验公式中因统计关系而存在的不确定性，应用范围广。具有计算速度快、稳定性好的优点，可直接用于测井数据解释工作。

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

Claims (7)

1.一种砂岩地层横波速度计算方法，其特征在于：所述方法包括：

(1)输入砂岩地层的实测纵波速度V_p、密度、孔隙度φ和泥质含量V_sh曲线；给定固结系数c的初始值和间隔；

(2)计算干岩石弹性模量；

(3)计算饱和岩石弹性模量，并计算得到计算纵波速度；

(4)建立实测纵波速度V_p与所述计算纵波速度的误差函数；

(5)在固结系数c的取值范围内以所述间隔改变固结系数c的值，每改变一次，重复一次步骤(2)至步骤(4)，得到每个新固结系数c对应的误差函数，然后比较所有固结系数c对应的误差函数，找到使得误差函数最小的固结系数c′，利用该使得误差函数最小的固结系数c′计算得到的横波速度即为所求的横波速度。

2.根据权利要求1所述的砂岩地层横波速度计算方法，其特征在于：所述步骤(2)中的干岩石弹性模量包括体积模量和剪切模量，分别由公式(1)和(2)计算得到：

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

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

其中，K₀和μ₀分别为固体矿物的体积模量和剪切模量，K_dry和μ_dry分别为干岩石体积模量和剪切模量；

所述公式(1)中，固体矿物的体积模量K₀由下述公式(3)计算得到：

K₀=[V_shK_sh+V_sK_s+1/(V_sh/K_sh+V_s/K_s)]/2     (3)

式中，K₀、K_s、K_sh分别代表固体矿物的体积模量、砂岩的体积模量、泥岩的体积模量；

V_sh和V_s分别为泥岩和砂岩的体积百分含量，由测井数据的泥质含量得到；

砂岩的体积模量K_s和泥岩的体积模量K_sh由岩石物理实验得到；

所述公式(1)中，固体矿物的剪切模量μ₀由下述公式(4)计算得到：

μ₀=[V_shμ_sh+V_sμ_s+1/(V_sh/μ_sh+V_s/_μs)]/2       (4)

式中，μ₀、μ_s、μ_sh分别代表固体矿物的剪切模量、砂岩的剪切模量、泥岩的剪切模量，其中砂岩的剪切模量μ_s和泥岩的剪切模量μ_sh由岩石物理实验得到。

3.根据权利要求2所述的砂岩地层横波速度计算方法，其特征在于：所述步骤(3)中的所述饱和岩石弹性模量包括体积模量和剪切模量，由公式(5)计算得到：

$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],{\mathrm{\μ}}_{\mathrm{sat}}={\mathrm{\μ}}_{\mathrm{dry}}---\left(5\right)$

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

4.根据权利要求3所述的砂岩地层横波速度计算方法，其特征在于：所述步骤(3)中的计算纵波速度由公式(6)计算得到：

$V_p^{\mathrm{est}}=\sqrt{(K_{\mathrm{ast}}+{4\mathrm{\μ}}_{\mathrm{sat}}/3)/{\mathrm{\ρ}}_{\mathrm{sat}}},V_s^{\mathrm{est}}=\sqrt{{\mathrm{\μ}}_{\mathrm{sat}}/{\mathrm{\ρ}}_{\mathrm{sat}}}---\left(6\right)$

其中，为纵波速度，为横波速度，ρ_sat为饱含流体的岩石密度。

5.根据权利要求4所述的砂岩地层横波速度计算方法，其特征在于：所述步骤(4)中的误差函数由公式(7)确定：

$f\left(c\right)={\left|\right|V_p-V_p^{\mathrm{est}}\left|\right|}_2---\left(7\right).$

6.根据权利要求5所述的砂岩地层横波速度计算方法，其特征在于：所述步骤(5)中，对于砂岩，固结系数c的取值范围取2≤c≤20。

7.根据权利要求6所述的砂岩地层横波速度计算方法，其特征在于：所述步骤(5)中利用该使得误差函数最小的固结系数c′计算得到的横波速度即为所求的横波速度是这样实现的：

利用所述使得误差函数最小的固结系数c′，重复步骤(2)和步骤(3)得到横波速度