CN103576195A - 一种随压力变化的裂隙介质横波速度预测方法 - Google Patents

Abstract

一种随压力变化的裂隙介质横波速度预测方法，采集测井数据，建立随压力变化的孔隙介质横波速度预测模型；将校正过的各个参数代入随压力变化的孔隙介质横波速度预测模型得到孔隙介质的与压力相关的纵横波速度；根据孔隙介质下的与压力相关的纵波速度、横波速度得到裂隙介质对称方向的与压力相关的纵波速度、横波速度并代入Thomson公式中，得到平行于裂隙方向的与压力相关的横波速度和垂直于裂隙方向的与压力相关的横波速度。本发明既考虑了压力对于介质横波速度的影响，也考虑到了储层中含有裂隙，更加符合实际。根据本发明得到的裂隙介质横波速度，可以更好的评估油气藏岩石属性，包括孔隙度、密度、岩性与流体含量，可以建立AVO模。

Description

一种随压力变化的裂隙介质横波速度预测方法

技术领域

本发明属于地震勘探中岩石物理领域，具体涉及一种随压力变化的裂隙介质横波速度预测方法。

背景技术

人们通常用纵横波速度比为常数来代替横波速度，但是对于不同的介质纵横波速度比应该是不同的，这样的假设不合理。因此很多横波速度的预测方法随之而出，有经验公式，也有基于岩石物理理论的。

比较著名的经验公式包括Castagna等人在1985年提出的著名的泥岩线为Vs=0.862Vp-1.172。Gardner于1974年给出了不同的岩性之间的速度与密度的关系，其中它的平均变换式为ρ=0.23V^0.25，这个平均变换式是对所有岩性的速度与密度关系的最佳拟合，它适合于所有岩性，不仅仅适用用某种岩性。而Castagna于1993年又对Gardner的公式进行了扩充，得到了不同岩性的速度与密度之间的关系：对于砂岩有

对于页岩有

对于石灰岩有

对于白云岩有

对硬石膏则是

Wyllie等人在1958年和1963年陆续提出了充满盐水的孔隙介质的孔隙度与速度之间的经验关系：1/V=(1-φ)/V_ma+φ/V_f1，其中V为岩石的整体速度，V_ma为岩石骨架的速度，V_f1为孔隙流体的速度，φ为孔隙度。当某些测井曲线缺失的时候，或者地震振幅异常都可以应用这些期限进行质量监控，但是这些公式对于岩性的依赖很强，而且依赖于局部条件，而Mavko等人多次在他们的书中提到：“这些关系公式都是经验公式，因而严格来说它们只能用在当时研究的岩石上”，因此这些经验公式不具有普遍性。

随着岩石物理理论的完善，基于岩石物理理论的横波速度预测方法逐渐成为研究的主流。如Greenberg和Castagna于1992年利用Biot-GassmannTheory（BGT）进行横波速度的预测，也就是假设在纵横波速度之间存在一个稳固的关系同时假设固体岩石成分之间的混合定律是线性的。Xu和White将Kuster和

于1974年建立的理论与差分有效介质理论结合，进行岩石弹性模量的计算，具体表现为利用孔隙纵横比来表征砂泥成分之间的关系，然后利用Gassmann方程来预测横波速度，该方法后来又被引申到预测裂隙介质的横波速度预测。Nolen-Hoeksema与Wang Zhijing于1996年根据实验室测得的干岩石的横波速度，利用Gassmann方程计算出干岩石的体变模量与切变模量，进而用到流体饱和岩石的横波速度预测中。在2006年的时候Lee提出的用固结参数联系基质弹性模量与骨架弹性模量之间的关系，通过实测纵波速度与预测纵波速度的比较得出固结参数，然后利用固结参数计算横波速度。2008年孙福利等人利用实际的数据对Lee的方法进行了验证，并且提出了固结参数的取值范围。

但是现有的这些方法没有考虑裂隙介质中横波速度分裂问题，更没有考虑到压力对裂隙介质属性的影响。而在实际生产中，储层存在裂隙，或者随着油田开发过程中的注水压裂等而产生裂隙，这些裂隙显然是受储层压力而发生变化的。在油田注水开发及CCS（碳捕捉与封存）技术中，注入井点的压力很高，随着水或者CO₂的高压注入，注入井附近可能产生新裂缝或者将原有裂缝撑大。而注CO2驱油时，生产井点则由于控制孔隙压力而可能导致裂隙的闭合。此时，注入井与生产井点储层内压力诱导下的快慢横波速度会发生变化。如何预测出不同孔隙压力下裂隙介质的快慢波速度，是利用四维地震监测CO₂在地下封存安全性的基础。无论是利用四维地震解释裂隙介质还是四维方位AVO反演、各向异性反演，以及四维地震转换波解释，都需要利用随孔隙压力变化的快慢横波速度才能进行。而同时考虑了压力变化和裂隙的快慢横波速度的地震解释也才更加符合实际。

发明内容

本发明的目的在于提供一种随压力变化的孔隙加裂隙介质的横波速度预测方法，预测得到的横波速度更加符合实际。

为实现上述目的，本发明采用如下的技术方案：

本发明包括以下步骤：

1）采集测井数据：采集总孔隙度φ_z、岩石的体积密度ρ以及实际纵波速度Vp、流体饱和度；采集岩石骨架的体变模量K_ma，岩石骨架的切变模量μ_ma，差异压力p，测量配位数C_p'，变形之前接触区域的半径a与颗粒的半径R；

分别对总孔隙度φ_z、岩石的体积密度ρ进行校正得到孔隙介质的孔隙孔隙度φ_p、孔隙介质的体积密度ρ_p，利用流体饱和度得到混合流体的体变模量K_f；计算岩石骨架的泊松比ν_ma；

对测量配位数C_p'进行线性拟合并且加权，得到含有加权系数W的配位数C_p，然后利用Digby、Mindlin公式以及Gassmann方程建立随压力变化的孔隙介质横波速度预测模型；

2）对岩心进行岩石物理测试，得到岩心的纵波速度、岩心的横波速度、岩心的密度以及岩心的孔隙度，然后利用随压力变化的孔隙介质横波速度预测模型对岩心的纵波速度、岩心的横波速度、岩心的密度以及岩心的孔隙度进行计算得到岩心的加权系数W_core，然后利用岩心的加权系数W_core对加权系数W进行校正，得到校正后的各向同性介质的加权系数W_isotropy，如公式（14）所示：

$W_{\mathrm{isotropy}}=W*(1-\frac{W_{\log }-W_{\mathrm{core}}}{W_{\log }})---\left(14\right)$

其中，W_log为岩心对应层位的加权系数；

3）将校正后的各向同性介质的加权系数W_isotropy、孔隙介质的孔隙孔隙度φ_p以及孔隙介质下的密度ρ_p，代入Gassmann方程的变形公式，得到孔隙介质下随压力变化的纵波速度α以及横波速度β；

4）根据孔隙介质下随压力变化的纵波速度α以及横波速度β，得到对称方向的弹性介质的纵波速度α₀、横波速度β₀；将对称方向的弹性介质的纵波速度α₀、横波速度β₀代入公式（15）-（17）中，得到平行于裂隙方向的横波速度V_s1和垂直于裂隙方向的横波速度V_s2；

$V_p^2\left(\mathrm{\θ}\right)={\mathrm{\α}}_0^2[1+2{\mathrm{\δ}\sin }^2\mathrm{\θ}{\cos }^2\mathrm{\θ}+2\mathrm{\ϵ}{\sin }^4\mathrm{\θ}]---\left(15\right)$

$V_{s1}^2\left(\mathrm{\θ}\right)={\mathrm{\β}}_0^2[1+2\frac{{\mathrm{\α}}_0^2}{{\mathrm{\β}}_0^2}{(\mathrm{\ϵ}-\mathrm{\δ})\sin }^2\mathrm{\θ}{\cos }^2\mathrm{\θ}]---\left(16\right)$

$V_{s2}^2\left(\mathrm{\θ}\right)={\mathrm{\β}}_0^2[1+2\mathrm{\γ}{\sin }^2\mathrm{\θ}]---\left(17\right)$

其中θ为波前与对称轴的夹角，ε、γ和δ为各向异性参数，Vp(θ)为预测纵波速度；

5）根据步骤4）得到的平行于裂隙方向的横波速度V_s1和垂直于裂隙方向的横波速度V_s2，建立各向异性介质的AVO模型，预测油气藏岩石属性。

所述步骤1）中若取岩心的井与进行采集测井数据的井是两口井，则标定岩心在进行采集数据的井中的层位。

所述步骤1）中随压力变化的孔隙介质横波速度预测模型建立的具体过程为：

a）对测量配位数C_p'进行线性拟合，得到测量配位数C_p'与孔隙度φ的关系为：C_p'=11.759e^1-φ-12.748          （1）

对测量配位数C_p'进行进行加权，得到配位数C_p，如公式（2）所示，

C_p=W*C_p'          （2）

b）在均匀的孔隙介质的条件下，利用Gassmann方程的变形进行纵横波速度的预测；公式（3）-（5）为Gassmann方程的变形为：

$V_{p_{\mathrm{sat}}}=\sqrt{\frac{K_{\mathrm{dry}}+\frac{{(1-\frac{K_{\mathrm{dry}}}{K_{\mathrm{ma}}})}^2}{\frac{\mathrm{\φ}}{K_f}+\frac{1-\mathrm{\φ}}{K_{\mathrm{ma}}}-\frac{K_{\mathrm{dry}}}{K_{\mathrm{ma}}^2}}+\frac{4}{3}{\mathrm{\μ}}_{\mathrm{sat}}}{\mathrm{\ρ}}}---\left(3\right)$

μ_sat=μ_dry（4）

$V_{s_{\mathrm{sat}}}=V_{\mathrm{dry}}=\sqrt{\frac{{\mathrm{\μ}}_{\mathrm{sat}}}{\mathrm{\ρ}}}=\sqrt{\frac{{\mathrm{\μ}}_{\mathrm{dry}}}{\mathrm{\ρ}}}---\left(5\right)$

其中

分别为预测的纵波速度、横波速度，μ_sat为孔隙介质的切变模量，K_dry为干岩石的体变模量、μ_dry为干岩石的切变模量，K_ma为岩石骨架的体变模量；K_f为混合流体的体变模量，ρ为岩石的体积密度；

将配位数C_p代入Digby公式中得到公式（7）：

$K_{\mathrm{dry}}=\frac{W*{C_p}^{\′}(1-\mathrm{\φ}){\mathrm{\μ}}_{\mathrm{ma}}b}{3\mathrm{\πR}(1-v)}---\left(7\right)$

将配位数C_p代入Mindlin公式中得到公式（9）：

${\mathrm{\μ}}_{\mathrm{dry}}=\frac{W*{C_p}^{\′}(1-\mathrm{\φ})}{20\mathrm{\πR}}(\frac{4{\mathrm{\μ}}_{\mathrm{ma}}b}{1-v_{\mathrm{ma}}}+\frac{12{\mathrm{\μ}}_{\mathrm{ma}}a}{2-v_{\mathrm{ma}}})---\left(9\right)$

其中，ν_ma为岩石骨架的泊松比，μ_ma为岩石骨架的切变模量；φ为孔隙度；C_p'为测量配位数；a为岩石颗粒变形之前接触区域的半径，b为岩石颗粒变形之后接触区域的半径，R为岩石颗粒的半径；

将干岩石的体变模量K_dry与干岩石的切变模量μ_dry代入Gassmann方程的变形公式中计算预测纵波速度

利用预测的纵波速度值等于实际测量的纵波速度值，得到加权系数W；

c）将加权系数W代入公式（2），得到配位数C_p，再把配位数C_p值代入计算随压力变化的干岩石切变模量的Mindlin公式中得到干岩石的切变模量μ_dry，再将干岩石的切变模量μ_dry代入Gassmann方程计算得到预测横波速度

即建立起了随压力变化的孔隙介质横波速度预测模型。

所述步骤b）中Digby公式为：

$K_{\mathrm{dry}}=\frac{C_p(1-\mathrm{\φ}){\mathrm{\μ}}_{\mathrm{ma}}b}{3\mathrm{\πR}(1-v)}---\left(6\right).$

所述步骤b）和步骤c）中Mindlin公式为：

${\mathrm{\μ}}_{\mathrm{dry}}=\frac{C_p(1-\mathrm{\φ})}{20\mathrm{\πR}}(\frac{4{\mathrm{\μ}}_{\mathrm{ma}}b}{1-v}+\frac{12{\mathrm{\μ}}_{\mathrm{ma}}a}{2-v})---\left(8\right).$

所述步骤b）中

$\frac{b}{R}={[d^2+{\left(\frac{a}{R}\right)}^2]}^{\frac{1}{2}}---\left(11\right)$

$d^3+\frac{3}{2}{\left(\frac{a}{R}\right)}^{2}d-\frac{3\mathrm{\π}(1-v)p}{2C_p(1-\mathrm{\φ}){\mathrm{\μ}}_{\mathrm{ma}}}=0---\left(12\right).$

其中p为差异压力。

所述步骤3）中将总孔隙度φ_z校正为孔隙介质下的孔隙孔隙度φ_p的具体过程为：将总孔隙度φ_z减去裂隙孔隙度φ_c得到孔隙孔隙度φ_p。

所述步骤1）中对岩石的体积密度ρ进行校正，得到孔隙介质情况下的孔隙介质的密度ρ_p的具体过程如下：

首先利用公式ρ_d=ρ^*(1-φ_z)计算得到干岩石的体积密度，其中ρ_d为干岩石的体积密度，然后利用公式ρ_p=ρ_d/(1-φ_p)得到孔隙介质的体积密度ρ_p。

所述步骤4）中对称方向的弹性介质的纵波速度α₀、横波速度β₀是通过以下过程得到的：

根据公式（18）和（19），计算得到对称方向的弹性介质的纵波速度α₀、横波速度β₀；

$\mathrm{\α}={\mathrm{\α}}_0*{\left[\frac{1+2*{(1-v)}^2/(1-2*v)}{1+2*\mathrm{\ϵ}}\right]}^{\frac{1}{2}}---\left(18\right)$

$\mathrm{\β}={\mathrm{\β}}_0*{(1+2*\mathrm{\γ})}^{\frac{1}{2}}---\left(19\right)$

其中α、β分别为孔隙介质下随压力变化的纵波速度、横波速度；ν为孔隙介质的泊松比，ε以及γ为各向异性参数。

所述步骤4）中各向异性参数ε、γ和δ的表达式分别为：

$\mathrm{\ϵ}=\frac{8}{3}(1-\frac{K_f}{K_{\mathrm{ma}}})D_{\mathrm{ci}}\left[\frac{(1-v_d^2)E}{(1-v^2)E_d}\right]{\mathrm{\η}}_c---\left(20\right)$

$\mathrm{\γ}=\frac{8}{3}\left(\frac{1-v_d}{2{-v}_d}\right){\mathrm{\η}}_c---\left(21\right)$

$\mathrm{\δ}=2(1-v)\mathrm{\ϵ}-2\left(\frac{1-2v}{1-v}\right)\mathrm{\γ}---\left(22\right)$

公式（20）和（21）中的 ${\mathrm{\η}}_c=\frac{3}{4\mathrm{\Π}}\frac{{\mathrm{\φ}}_c}{(c/a)}---\left(26\right)$

公式（20）中Dci为流体影响因子，在低频和高频声波测量时分别如下：

公式（23）中 $A_c\left(v_d\right)=\frac{16}{9}\left(\frac{1-v_d^2}{1-2v_d}\right)---\left(25\right)$

其中，K_ma为岩石骨架的体变模量，K_f为混合流体的体变模量，ν为孔隙介质的泊松比，E为杨氏模量，ν_d和E_d为相应的各向同性孔隙干岩石的泊松比和杨氏模量，ν_ma则为岩石骨架的泊松比，φ_z为总孔隙度，D_ci(lo)为低频声波测量时的流体影响因子，D_ci(mh)为高频声波测量时的流体影响因子。

本发明的有益效果是：本发明针对现有技术中利用孔隙介质横波速度预测模型预测所有储层横波速度与实际情况不符的问题，本发明考虑了储层中发育裂隙的情况，以建立的随压力变化孔隙介质横波速度预测模型为基础，利用岩石物理测试资料对预测的孔隙介质纵横波速度等参数进行校正，然后利用Thomson裂隙介质理论进行横波速度的预测，这样更加符合存在裂隙的储层。由于横波速度对压力是敏感的，本发明既考虑了压力对于介质横波速度的影响，同时也考虑到了储层中含有裂隙，更加符合实际。根据本发明得到的裂隙介质横波速度，可以建立AVO模型，进行转换波分析，通过AVO分析，可以更好地评估油气藏岩石属性，包括孔隙度、密度、岩性与流体含量。

附图说明

图1为C语言编程流程图。

图2为岩心为干岩石进行校正，得到的快横波速度与实测快横波速度的对比，其中实线为实测值，虚线为预测值。

图3为岩心为干岩石进行校正，得到的慢横波速度与实测慢横波速度的对比，其中实线为实测值，虚线为预测值。

图4为岩心为油水饱和岩石进行校正，得到的快横波速度与实测快横波速度的对比，其中实线为实测值，虚线为预测值。

图5为岩心为油水饱和岩石进行校正后得到的慢横波速度与实测慢横波速度的对比，其中实线为实测值，虚线为预测值。

具体实施方式

下面结合附图对本发明做详细描述。本发明中的*表示乘号。本发明利用采集数据，得到测井资料，进行储层横波速度的预测，具体涉及一种利用Thomson裂隙介质理论结合随压力变化的孔隙介质模型以及岩石物理测试资料进行含裂隙的储层横波速度预测的方法。

本发明具体包括以下步骤：

1）本发明利用考虑压力变化的Digby、Mindlin公式以及Gassmann方程的变形建立的随压力变化孔隙介质横波速度预测模型来预测储层介质的配位数，然后利用岩石物理测试资料来校正该配位数，得到更为准确的孔隙介质的配位数，以及纵横波速度。以孔隙介质的纵横波速度为基础，利用Thomson裂隙介质理论进行横波速度的预测。

需要说明一下孔隙度的问题，在裂隙介质中总孔隙度为裂隙孔隙度和孔隙孔隙度的和；而下面建立的是孔隙介质模型，介质为孔隙介质，所以只有一个孔隙度；建立的随压力变化孔隙介质横波速度预测模型的具体过程为：

a）采集测井数据：采集总孔隙度φ_z、岩石的体积密度ρ以及实际纵波速度Vp、饱和度；采集岩石骨架的体变模量K_ma，岩石骨架的切变模量μ_ma,差异压力p，测量配位数C_p'，岩石颗粒变形之前接触区域的半径a与岩石颗粒的半径R；根据公式

计算得到岩石骨架的泊松比ν_ma，其中K_ma为岩石骨架的体变模量，μ_ma为岩石骨架的切变模量。

根据测量配位数C_p'与e^1-φ成正比的结论，进行线性拟合得到测量配位数与孔隙度的关系为：C_p'=11.759e^1-φ-12.748          （1）

其中C_p'为测量配位数，φ为孔隙度；

在实际情况中，对测量配位数C_p'进行加权，得到配位数C_p，如公式（2）所示，其中W为加权系数：

C_p=W*C_p'         （2）；

b）假设介质为均匀的孔隙介质，利用Gassmann方程的变形进行纵横波速度的预测；公式（3）-（5）为Gassmann方程的变形：

$V_{p_{\mathrm{sat}}}=\sqrt{\frac{K_{\mathrm{dry}}+\frac{{(1-\frac{K_{\mathrm{dry}}}{K_{\mathrm{ma}}})}^2}{\frac{\mathrm{\φ}}{K_f}+\frac{1-\mathrm{\φ}}{K_{\mathrm{ma}}}-\frac{K_{\mathrm{dry}}}{K_{\mathrm{ma}}^2}}+\frac{4}{3}{\mathrm{\μ}}_{\mathrm{sat}}}{\mathrm{\ρ}}}---\left(3\right)$

μ_sat=μ_dry            （4）

$V_{s_{\mathrm{sat}}}=V_{\mathrm{dry}}=\sqrt{\frac{{\mathrm{\μ}}_{\mathrm{sat}}}{\mathrm{\ρ}}}=\sqrt{\frac{{\mathrm{\μ}}_{\mathrm{dry}}}{\mathrm{\ρ}}}---\left(5\right)$

其中

与

分别为预测的纵波速度、横波速度，μ_sat为孔隙介质的切变模量，K_dry与μ_dry分别为干岩石的体变模量和切变模量，K_ma为岩石骨架的体变模量；φ为岩石的孔隙度；K_f为混合流体的体变模量，ρ为岩石的体积密度；其中岩石骨架的体变模量K_ma、孔隙度φ以及岩石的体积密度ρ均从采集的测井数据资料中获得，混合流体的体变模量K_f由流体的饱和度得到；

计算随压力变化的压实紧密的干岩石的体变模量的Digby公式为：

$K_{\mathrm{dry}}=\frac{C_p(1-\mathrm{\φ}){\mathrm{\μ}}_{\mathrm{ma}}b}{3\mathrm{\πR}(1-v)}---\left(6\right).$

将公式（2）中的配位数C_p公式代入计算随压力变化的压实紧密的干岩石的体变模量的Digby公式得到公式（7）：

$K_{\mathrm{dry}}=\frac{W*{C_p}^{\′}(1-\mathrm{\φ}){\mathrm{\μ}}_{\mathrm{ma}}b}{3\mathrm{\πR}(1-v)}---\left(7\right)$

计算随压力变化的干岩石切变模量的Mindlin公式为：

${\mathrm{\μ}}_{\mathrm{dry}}=\frac{C_p(1-\mathrm{\φ})}{20\mathrm{\πR}}(\frac{4{\mathrm{\μ}}_{\mathrm{ma}}b}{1-v}+\frac{12{\mathrm{\μ}}_{\mathrm{ma}}a}{2-v})---\left(8\right).$

将公式（2）中的的配位数C_p公式计算随压力变化的干岩石切变模量的Mindlin公式中得到公式（9）：

${\mathrm{\μ}}_{\mathrm{dry}}=\frac{W*{C_p}^{\′}(1-\mathrm{\φ})}{20\mathrm{\πR}}(\frac{4{\mathrm{\μ}}_{\mathrm{ma}}b}{1-v_{\mathrm{ma}}}+\frac{12{\mathrm{\μ}}_{\mathrm{ma}}a}{2-v_{\mathrm{ma}}})---\left(9\right)$

其中， $\frac{b}{R}={[d^2+{\left(\frac{a}{R}\right)}^2]}^{\frac{1}{2}}---\left(11\right)$

$d^3+\frac{3}{2}{\left(\frac{a}{R}\right)}^{2}d-\frac{3\mathrm{\π}(1-v)p}{2C_p(1-\mathrm{\φ}){\mathrm{\μ}}_{\mathrm{ma}}}=0---\left(12\right)$

φ为孔隙度；C_p'为测量配位数；p为差异压力；a为岩石颗粒变形之前接触区域的半径，b为岩石颗粒变形之后接触区域的半径，R为岩石颗粒的半径；差异压力p的变化会体现在

中，即公式（11），无论是在公式（11）中求得差异压力p值还是在Digby、Mindlin公式中利用差异压力p值，在实际计算中

是作为一个整体的；

将K_dry与μ_dry代入Gassmann方程的变形公式中计算预测纵波速度

由于实际采集的数据中有纵波速度，预测的纵波速度值要接近实际测量的纵波速度值，因此利用公式（13）即实际测量纵波速度值Vp_measured减去预测纵波速度值Vp_sat(W)等于零（从理论上来说预测的横波速度要无限接近实际的横波速度，而在实际上，无限接近无法实际操作，因此假设实际测量纵波速度值Vp_measured与预测纵波速度值Vp_sat(W)完全相等），得到加权系数W；

Vp_measured-Vp_por(W)→0                （13）

c）将加权系数W代入公式（2），得到配位数C_p，得到配位数C_p之后，把配位数C_p值代入计算随压力变化的干岩石切变模量的Mindlin公式中得到干岩石的切变模量μ_dry，进而将干岩石的切变模量μ_dry代入Gassmann方程即代入公式（5）计算得到预测横波速度

即建立起了随压力变化的孔隙介质横波速度预测模型。

本发明利用已经建立好的随压力变化的孔隙介质横波速度预测模型，预测实际资料的加权系数W、配位数C_p以及预测随压力变化的孔隙介质横波速度，预测的加权系数W见表1。

2）由于在孔隙介质模型中，假设储层为孔隙介质，因此所用到的实际的采集的数据也假设为孔隙介质的测量结果，而当储层含有裂隙时这与实际情况是不符合的，因此本方法利用岩石物理测试资料对上述预测结果进行校正。

进行岩石物理测试的岩心相对于整个储层来说是各向同性的；如果该岩心所在的井与测井所在的井为同一口井，那么不需要标定层位，如果是两口井（即取岩心的井与进行数据采集的井不是同一口井），则需要根据地质资料标定岩心在采集测井数据所在井中的层位。对岩心进行岩石物理测试，得到岩心的纵波速度、横波速度、密度以及孔隙度参数，然后利用随压力变化的孔隙介质横波速度预测模型对这些数据进行计算得到岩心的加权系数W_core，然后利用岩心的加权系数W_core对加权系数W进行校正，得到校正后的各向同性介质的加权系数W_isotropy，同时，在利用岩石物理测试数据的时候，使岩心的差异压力与储层的差异压力要相同。标定好层位之后，利用公式（14）对加权系数W进行校正，得到校正后各向同性介质的加权系数W_isotropy：

$W_{\mathrm{isotropy}}=W*(1-\frac{W_{\log }-W_{\mathrm{core}}}{W_{\log }})---\left(14\right)$

其中W_isotropy为校正后的各向同性介质的加权系数，W_core为计算得到的岩心的加权系数，W_log为岩心对应层位的加权系数。校正之前的加权系数W利用随压力变化的孔隙介质模型预测出来。利用干岩石物理测试结果校正后得到的校正后各向同性介质的加权系数W_isotropy与利用饱和岩石物理测试结果校正后得到的W_isotropy见表1。表1为校正之前的加权系数W以及利用干岩石以及油水饱和岩石物理测试结果校正得到的加权系数W。校正之前的加权系数为利用孔隙介质模型预测的加权系数W值。该储层进行岩心岩石物理测试的时候，既对干岩石进行了测试，也对油水饱和岩石进行了测试。因此分别利用干岩石的预测结果和油水饱和岩石的预测结果进行校正。

表1预测的加权系数W

3）利用裂隙孔隙度将总孔隙度校正为孔隙孔隙度，同时对岩石的体积密度进行校正，得到孔隙介质的体积密度，具体过程为：将由测井数据得到的总孔隙度φ_z减去裂隙孔隙度φ_c得到孔隙孔隙度φ_p；然后对体积密度进行校正，方法如下：首先利用下述公式ρ_d=ρ^*(1-φ_z)计算得到干岩石的体积密度，其中ρ_d为干岩石的体积密度，ρ为测井数据资料得到的体积密度，然后利用公式ρ_p=ρ_d/(1-φ_p)得到孔隙介质的体积密度，其中ρ_p为校正后得到的孔隙介质的体积密度，φ_p为孔隙孔隙度；然后将得到的校正后的各向同性介质的加权系数W_isotropy，以及校正过的孔隙度，体积密度代入孔隙介质横波速度预测模型，得到孔隙介质下随压力变化的纵波速度α以及横波速度β；

4）然后利用Thomson裂隙介质理论进行横波速度的预测；上述孔隙介质纵横波速度的预测以及Thomson裂隙介质预测横波速度的实现是通过C语言编程实现的，具体流程见图1。Thomsen于1986年提出了裂隙为垂向裂隙且有一个对称轴的弱各向异性的弹性介质的相速度公式如下：

$V_p^2\left(\mathrm{\θ}\right)={\mathrm{\α}}_0^2[1+2{\mathrm{\δ}\sin }^2\mathrm{\θ}{\cos }^2\mathrm{\θ}+2\mathrm{\ϵ}{\sin }^4\mathrm{\θ}]---\left(15\right)$

$V_{s1}^2\left(\mathrm{\θ}\right)={\mathrm{\β}}_0^2[1+2\frac{{\mathrm{\α}}_0^2}{{\mathrm{\β}}_0^2}{(\mathrm{\ϵ}-\mathrm{\δ})\sin }^2\mathrm{\θ}{\cos }^2\mathrm{\θ}]---\left(16\right)$

$V_{s2}^2\left(\mathrm{\θ}\right)={\mathrm{\β}}_0^2[1+2\mathrm{\γ}{\sin }^2\mathrm{\θ}]---\left(17\right)$

其中α₀、β₀分别为对称方向的孔隙的纵波速度、横波速度，可以利用上述利用孔隙介质模型预测的纵波速度、横波速度以及公式（18）和公式（19）获得；ε、γ和δ为各向异性参数；θ为波前与对称轴的夹角。V_s1是平行于裂隙方向的横波速度，而V_s2则是垂直于裂隙方向的横波速度；这两个横波速度不同（即快横波速度与慢横波速度，在裂隙介质中横波会发生分裂，产生两个横波，按照速度划分，就可以成为快横波与慢横波，相应的速度就是快横波速度与慢横波速度），实际测井时，是先后接收到两个横波速度，因此称其为快横波和慢横波。

$\mathrm{\α}={\mathrm{\α}}_0*{\left[\frac{1+2*{(1-v)}^2/(1-2*v)}{1+2*\mathrm{\ϵ}}\right]}^{\frac{1}{2}}---\left(18\right)$

$\mathrm{\β}={\mathrm{\β}}_0*{(1+2*\mathrm{\γ})}^{\frac{1}{2}}---\left(19\right)$

其中α，β为上述孔隙介质模型得到的随压力变化的孔隙介质的纵横波速度，α₀、β₀分别为储层对称方向的纵波速度、横波速度，ν为孔隙介质的泊松比，ε以及γ为各向异性参数。

其中各向异性参数分别为：

$\mathrm{\ϵ}=\frac{8}{3}(1-\frac{K_f}{K_{\mathrm{ma}}})D_{\mathrm{ci}}\left[\frac{(1-v_d^2)E}{(1-v^2)E_d}\right]{\mathrm{\η}}_c---\left(20\right)$

$\mathrm{\γ}=\frac{8}{3}\left(\frac{1-v_d}{2{-v}_d}\right){\mathrm{\η}}_c---\left(21\right)$

$\mathrm{\δ}=2(1-v)\mathrm{\ϵ}-2\left(\frac{1-2v}{1-v}\right)\mathrm{\γ}---\left(22\right)$

其中Dci为流体影响因子，在低频和高频声波测量时分别如下：

$D_{\mathrm{ci}}\left(\mathrm{lo}\right)={[1-\frac{K_f}{K_{\mathrm{ma}}}+\frac{K_f}{K_d{\mathrm{\φ}}_z}((1-\frac{K_d}{K_f})+A_c\left(v_d\right){\mathrm{\η}}_c)]}^{-1}---\left(23\right)$

$A_c\left(v_d\right)=\frac{16}{9}\left(\frac{1-v_d^2}{1-2v_d}\right)---\left(25\right)$

${\mathrm{\η}}_c=\frac{3}{4\mathrm{\Π}}\frac{{\mathrm{\φ}}_c}{(c/a)}---\left(26\right)$

在上述公式中，K_ma为岩石骨架的体变模量，K_f为裂隙中混合流体的体变模量，ν为孔隙介质的泊松比，E为杨氏模量，ν_d和E_d为相应的各向同性孔隙干岩石的泊松比和杨氏模量，ν_ma则为岩石骨架的泊松比，φ_z为总孔隙度也就是孔隙孔隙度φ_p与裂隙孔隙度φ_c的和。

由于θ是未知的，公式（15）中Vp(θ)为预测纵波速度，由于预测纵波速度是趋近与实测纵波速度Vp，可以利用Vp与Vp(θ)相减为0，得到θ。因此首先利用实测纵波速度和公式（15）预测出θ值，然后将所得到的θ代入公式（16）和公式（17）得到两个横波速度。

上述孔隙介质纵横波速度的预测、角度θ的求取以及Thomson裂隙介质预测横波速度的实现是通过C语言编程实现的，具体流程见图1，具体过程为：输入差异压力p以及校正过的参数，包括孔隙介质的孔隙孔隙度φ_p、孔隙介质的体积密度ρ_p以及校正过的加权系数W_isotropy，利用利用考虑压力变化的Digby、Mindlin公式以及Gassmann方程的变形公式计算随压力变化的孔隙介质的纵波速度α、横波速度β；然后利用一个循环来求出波前与对称轴的夹角θ，首先设定一个初值即定义波前与对称轴的夹角θ=1（角度），将该θ值和孔隙介质的孔隙孔隙度φ_p、孔隙介质的密度ρ_p以及孔隙介质的纵波速度Vp_measured代入预测纵波速度Vp(θ)，得到预测纵波速度Vp(θ)，然后将该波前与对称轴的夹角θ加上1得到新的θ值，然后判断新的θ值是否小于90度，如果小于90度，则继续将该新的θ值和孔隙介质的孔隙孔隙度φ_p、孔隙介质的密度ρ_p以及孔隙介质的纵波速度Vp_measured代入预测纵波速度Vp(θ)，然后在将θ值加上1判断新的θ值是否小于90度，若小于90度重复上述步骤，即将大于0度小于90度的每一个整数角度下的预测纵波速度Vp(θ)计算出来；若新的θ值不小于90度，则再继续下一步，即得到这89个预测纵波速度之后，将它们分别与实际纵波速度Vp相减，得到误差，找到误差最小的那个角度下的预测纵波速度Vp(θ)，这个角度即为波前与对称轴的夹角θ；然后将得到的角度以及孔隙介质的孔隙孔隙度φ_p、孔隙介质的密度ρ_p、随压力变化的孔隙介质的纵波速度α、横波速度β代入Thomson公式，得到裂隙介质的横波速度并且输出，C语言程序结束。

图2是岩心为干岩石进行校正即利用表1中校正过的W_isotropy值，得到的快横波速度与实测快横波速度的对比，其中实线为实测值，虚线为预测值，其平均误差为7.742%。图3是岩心为干岩石进行校正即利用表1中校正过的加权系数W值，得到的慢横波速度与实测慢横波速度的对比，其中实线为实测值，虚线为预测值，其平均误差为7.633%。图4是岩心为饱和岩石进行校正即利用表1中校正过的W值，得到的快横波速度与实测快横波速度的对比，其中实线为实测值，虚线为预测值，其平均误差为7.478%。图5中岩心为干岩石进行校正即利用表1中校正过的W值，得到的慢横波速度与实测慢横波速度的对比，其中实线为实测值，虚线为预测值，其平均误差为7.274%。

5）根据步骤4）得到的平行于裂隙方向的横波速度V_s1和垂直于裂隙方向的横波速度V_s2，建立各向异性介质的AVO模型，预测油气藏岩石属性，包括孔隙度、密度、岩性与流体含量。

通过AVO分析，地球物理学家可以更好地评估油气藏岩石属性，包括孔隙度、密度、岩性与流体含量，而横波速度是建立AVO模型、转换波分析过程中的一个不可缺少的弹性参数，而在多数情况下研究工区是没有横波速度数据的，因此本发明给出了一个横波速度预测的方法。储层中还可能发育裂隙，那么如果利用孔隙介质横波速度预测模型与实际情况不符，而本发明以已经建立好的随压力变化孔隙介质横波速度预测模型为基础，利用岩石物理测试资料对预测的孔隙介质纵横波速度等参数进行校正，然后利用Thomson裂隙介质理论进行横波速度的预测，这样更加符合存在裂隙的储层。横波速度对压力是敏感的，本发明既考虑了压力对于介质横波速度的影响，同时也考虑到了储层中含有裂隙，更加符合实际。

本发明采集测井数据，建立起随压力变化的孔隙介质横波速度预测模型；利用由岩心的测试资料得到的加权系数对实际采集到的数据得到的加权系数进行校正，得到校正后的各向同性介质的加权系数；将校正过的各个参数代入随压力变化的孔隙介质横波速度预测模型得到孔隙介质的与压力相关的纵横波速度；根据孔隙介质下的与压力相关的纵波速度、横波速度得到裂隙介质对称方向的与压力相关的纵波速度、横波速度并带入Thomson提出的公式中，得到平行于裂隙方向的与压力相关的横波速度V_s1和垂直于裂隙方向的与压力相关的横波速度V_s2。

Claims (10)

1.一种随压力变化的裂隙介质横波速度预测方法，其特征在于，包括以下步骤：

1）采集测井数据：采集总孔隙度φ_z、岩石的体积密度ρ以及实际纵波速度Vp、流体饱和度；采集岩石骨架的体变模量K_ma，岩石骨架的切变模量μ_ma，差异压力p，测量配位数C_p'，岩石颗粒变形之前接触区域的半径a与岩石颗粒的半径R；

分别对总孔隙度φ_z、岩石的体积密度ρ进行校正得到孔隙介质的孔隙孔隙度φ_p、孔隙介质的体积密度ρ_p，利用流体饱和度得到混合流体的体变模量K_f；计算岩石骨架的泊松比ν_ma；

对测量配位数C_p'进行线性拟合并且加权，得到含有加权系数W的配位数C_p，然后利用Digby、Mindlin公式以及Gassmann方程的变形公式建立随压力变化的孔隙介质横波速度预测模型；

2）对岩心进行岩石物理测试，得到岩心的纵波速度、岩心的横波速度、岩心的密度以及岩心的孔隙度，然后利用随压力变化的孔隙介质横波速度预测模型对岩心的纵波速度、岩心的横波速度、岩心的密度以及岩心的孔隙度进行计算得到岩心的加权系数W_core，然后利用岩心的加权系数W_core对加权系数W进行校正，得到校正后的各向同性介质的加权系数W_isotropy，如公式（14）所示：

$W_{\mathrm{isotropy}}=W*(1-\frac{W_{\log }-W_{\mathrm{core}}}{W_{\log }})---\left(14\right)$

其中，W_log为岩心对应层位的加权系数；

3）将校正后的各向同性介质的加权系数W_isotropy、孔隙介质的孔隙孔隙度φ_p以及孔隙介质下的密度ρ_p代入Gassmann方程的变形公式，得到孔隙介质下随压力变化的纵波速度α以及横波速度β；

4）根据孔隙介质下随压力变化的纵波速度α以及横波速度β，得到对称方向的弹性介质的纵波速度α₀、横波速度β₀；将对称方向的弹性介质的纵波速度α₀、横波速度β₀代入公式（15）-（17）中，得到平行于裂隙方向的横波速度V_s1和垂直于裂隙方向的横波速度V_s2；

$V_p^2\left(\mathrm{\θ}\right)={\mathrm{\α}}_0^2[1+2{\mathrm{\δ}\sin }^2\mathrm{\θ}{\cos }^2\mathrm{\θ}+2\mathrm{\ϵ}{\sin }^4\mathrm{\θ}]---\left(15\right)$

$V_{s1}^2\left(\mathrm{\θ}\right)={\mathrm{\β}}_0^2[1+2\frac{{\mathrm{\α}}_0^2}{{\mathrm{\β}}_0^2}{(\mathrm{\ϵ}-\mathrm{\δ})\sin }^2\mathrm{\θ}{\cos }^2\mathrm{\θ}]---\left(16\right)$

$V_{s2}^2\left(\mathrm{\θ}\right)={\mathrm{\β}}_0^2[1+2\mathrm{\γ}{\sin }^2\mathrm{\θ}]---\left(17\right)$

其中θ为波前与对称轴的夹角，ε、γ和δ为各向异性参数，Vp(θ)为预测纵波速度；

5）根据步骤4）得到的平行于裂隙方向的横波速度V_s1和垂直于裂隙方向的横波速度V_s2，建立各向异性介质的AVO模型，预测油气藏岩石属性。

2.一种随压力变化的裂隙介质横波速度预测方法，其特征在于，所述步骤1）中若取岩心的井与进行采集测井数据的井是两口井，则标定岩心在进行采集数据的井中的层位。

3.根据权利要求1所述的一种随压力变化的裂隙介质横波速度预测方法，其特征在于，所述步骤1）中随压力变化的孔隙介质横波速度预测模型建立的具体过程为：

a）对测量配位数C_p'进行线性拟合，得到测量配位数C_p'与孔隙度φ的关系为：C_p'=11.759e^1-φ-12.748               （1）

对测量配位数C_p'进行加权，得到配位数C_p，如公式（2）所示，

C_p=W*C_p'              （2）

b）在均匀的孔隙介质的条件下，利用Gassmann方程的变形进行纵横波速度的预测；公式（3）-（5）为Gassmann方程的变形为：

$V_{p_{\mathrm{sat}}}=\sqrt{\frac{K_{\mathrm{dry}}+\frac{{(1-\frac{K_{\mathrm{dry}}}{K_{\mathrm{ma}}})}^2}{\frac{\mathrm{\φ}}{K_f}+\frac{1-\mathrm{\φ}}{K_{\mathrm{ma}}}-\frac{K_{\mathrm{dry}}}{K_{\mathrm{ma}}^2}}+\frac{4}{3}{\mathrm{\μ}}_{\mathrm{sat}}}{\mathrm{\ρ}}}---\left(3\right)$

μ_sat=μ_dry（4）

$V_{s_{\mathrm{sat}}}=V_{\mathrm{dry}}=\sqrt{\frac{{\mathrm{\μ}}_{\mathrm{sat}}}{\mathrm{\ρ}}}=\sqrt{\frac{{\mathrm{\μ}}_{\mathrm{dry}}}{\mathrm{\ρ}}}---\left(5\right)$

其中

分别为预测的纵波速度、横波速度，μ_sat为孔隙介质的切变模量，K_dry为干岩石的体变模量、μ_dry为干岩石的切变模量，K_ma为岩石骨架的体变模量；K_f为混合流体的体变模量，ρ为岩石的体积密度；

将配位数C_p代入Digby公式中得到公式（7）：

$K_{\mathrm{dry}}=\frac{W*{C_p}^{\′}(1-\mathrm{\φ}){\mathrm{\μ}}_{\mathrm{ma}}b}{3\mathrm{\πR}(1-v)}---\left(7\right)$

将配位数C_p代入Mindlin公式中得到公式（9）：

${\mathrm{\μ}}_{\mathrm{dry}}=\frac{W*{C_p}^{\′}(1-\mathrm{\φ})}{20\mathrm{\πR}}(\frac{4{\mathrm{\μ}}_{\mathrm{ma}}b}{1-v_{\mathrm{ma}}}+\frac{12{\mathrm{\μ}}_{\mathrm{ma}}a}{2-v_{\mathrm{ma}}})---\left(9\right)$

其中，ν_ma为岩石骨架的泊松比，μ_ma为岩石骨架的切变模量；φ为孔隙度；C_p'为测量配位数；a为岩石颗粒变形之前接触区域的半径，b为岩石颗粒变形之后接触区域的半径，R为岩石颗粒的半径；

将干岩石的体变模量K_dry与干岩石的切变模量μ_dry代入Gassmann方程的变形公式中计算预测纵波速度

利用预测的纵波速度值等于实际测量的纵波速度值，得到加权系数W；

c）将加权系数W代入公式（2），得到配位数C_p，再把配位数C_p值代入计算随压力变化的干岩石切变模量的Mindlin公式中得到干岩石的切变模量μ_dry，再将干岩石的切变模量μ_dry代入Gassmann方程的变形公式中计算得到预测横波速度即建立起了随压力变化的孔隙介质横波速度预测模型。

4.根据权利要求3所述的一种随压力变化的裂隙介质横波速度预测方法，其特征在于，所述步骤b）中Digby公式为：

$K_{\mathrm{dry}}=\frac{C_p(1-\mathrm{\φ}){\mathrm{\μ}}_{\mathrm{ma}}b}{3\mathrm{\πR}(1-v)}---\left(6\right).$

5.根据权利要求3所述的一种随压力变化的裂隙介质横波速度预测方法，其特征在于，所述步骤b）和步骤c）中Mindlin公式为：

${\mathrm{\μ}}_{\mathrm{dry}}=\frac{C_p(1-\mathrm{\φ})}{20\mathrm{\πR}}(\frac{4{\mathrm{\μ}}_{\mathrm{ma}}b}{1-v}+\frac{12{\mathrm{\μ}}_{\mathrm{ma}}a}{2-v})---\left(8\right).$

6.根据权利要求3所述的一种随压力变化的裂隙介质横波速度预测方法，其特征在于，所述步骤b）中

$\frac{b}{R}={[d^2+{\left(\frac{a}{R}\right)}^2]}^{\frac{1}{2}}---\left(11\right)$

$d^3+\frac{3}{2}{\left(\frac{a}{R}\right)}^{2}d-\frac{3\mathrm{\π}(1-v)p}{2C_p(1-\mathrm{\φ}){\mathrm{\μ}}_{\mathrm{ma}}}=0---\left(12\right)$

其中p为差异压力。

7.根据权利要求1所述的一种随压力变化的裂隙介质横波速度预测方法，其特征在于，所述步骤3）中将总孔隙度φ_z校正为孔隙介质下的孔隙孔隙度φ_p的具体过程为：将总孔隙度φ_z减去裂隙孔隙度φ_c得到孔隙孔隙度φ_p。

8.根据权利要求1所述的一种随压力变化的裂隙介质横波速度预测方法，其特征在于，所述步骤1）中对岩石的体积密度ρ进行校正，得到孔隙介质情况下的孔隙介质的密度ρ_p的具体过程如下：

首先利用公式ρ_d=ρ^*(1-φ_z)计算得到干岩石的密度，其中ρ_d为干岩石的密度，然后利用公式ρ_p=ρ_d/(1-φ_p)得到孔隙介质的密度ρ_p。

9.根据权利要求1所述的一种随压力变化的裂隙介质横波速度预测方法，其特征在于，所述步骤4）中对称方向的弹性介质的纵波速度α₀、横波速度β₀是通过以下过程得到的：

根据公式（18）和（19），计算得到对称方向的弹性介质的纵波速度α₀、横波速度β₀；

$\mathrm{\α}={\mathrm{\α}}_0*{\left[\frac{1+2*{(1-v)}^2/(1-2*v)}{1+2*\mathrm{\ϵ}}\right]}^{\frac{1}{2}}---\left(18\right)$

$\mathrm{\β}={\mathrm{\β}}_0*{(1+2*\mathrm{\γ})}^{\frac{1}{2}}---\left(19\right)$

其中α、β分别为孔隙介质下随压力变化的纵波速度、横波速度；ν为孔隙介质的泊松比，ε以及γ为各向异性参数。

10.根据权利要求1所述的一种随压力变化的裂隙介质横波速度预测方法，其特征在于，所述步骤4）中各向异性参数ε、γ和δ的表达式分别为：

$\mathrm{\ϵ}=\frac{8}{3}(1-\frac{K_f}{K_{\mathrm{ma}}})D_{\mathrm{ci}}\left[\frac{(1-v_d^2)E}{(1-v^2)E_d}\right]{\mathrm{\η}}_c---\left(20\right)$

$\mathrm{\γ}=\frac{8}{3}\left(\frac{1-v_d}{2{-v}_d}\right){\mathrm{\η}}_c---\left(21\right)$

$\mathrm{\δ}=2(1-v)\mathrm{\ϵ}-2\left(\frac{1-2v}{1-v}\right)\mathrm{\γ}---\left(22\right)$

公式（20）和（21）中的 ${\mathrm{\η}}_c=\frac{3}{4\mathrm{\Π}}\frac{{\mathrm{\φ}}_c}{(c/a)}---\left(26\right)$

公式（20）中Dci为流体影响因子，在低频和高频声波测量时分别如下：

公式（23）中 $A_c\left(v_d\right)=\frac{16}{9}\left(\frac{1-v_d^2}{1-2v_d}\right)---\left(25\right)$

其中，K_ma为岩石骨架的体变模量，K_f为混合流体的体变模量，ν为孔隙介质的泊松比，E为杨氏模量，ν_d和E_d为相应的各向同性孔隙干岩石的泊松比和杨氏模量，ν_ma则为岩石骨架的泊松比，φ_z为总孔隙度，D_ci(lo)为低频声波测量时的流体影响因子，D_ci(mh)为高频声波测量时的流体影响因子。

Images