一种粘弹性介质中的球面PP波反射系数建模方法

技术领域

本发明属于石油地球物理勘探领域，具体涉及一种粘弹性介质中的球面PP波反射系数建模方法，从储层参数出发，基于White模型和球面波的平面波分解算法计算粘弹性介质中球面PP波反射系数的建模。

背景技术

传统的AVO(振幅随偏移距)或AVA(振幅随入射角变化)是基于平面波(几何地震学)的假设。然而，当震源或检波器离反射界面很近或者需要研究临界角附近的反射波场的时候，几何地震学近似就变得不准确。这是由于实际地震勘探中使用的是点源，激发的是球面波。(1961)最早提出了球面波反射系数的概念，并通过解析解的结果分析了球面波反射波场(1959；1960)以及临界角附近的反射波和首波的动力学特征(1962；1965)。基于球面波的平面波分解算法，Haase(2004)给出了层状均匀介质中球面波反射系数的计算方法，并讨论了第1类和第3类AVO对应的球面PP和PS波的反射系数振幅和相位特点。Ursenbach(2007)通过引入一种特殊形式的子波，使得在计算球面波反射系数过程中对频率的积分可以解析求得，大大提高了计算效率。Ayzenberg等人(2007，2009)给出了利用有效反射系数(ERCs)表示的球面波反射系数形式。

传统AVO(AVA)技术的另一个假设就是地下介质完全弹性假设。而VSP记录(Hauge，1981；Hedl in等，2001)、测井数据(Schmitt，1999)以及实验室中的岩石物理测量(Behura等，2007；Winkler和Nur，1982)均显示地震波在传播过程中会出现衰减和频散现象，尤其是对于含有碳氢化合物的地层，地震波的衰减更为明显。根据地震波传播理论，国外部分学者将其归因于地下介质 强吸收的反射系数。White(1975)，De Hoop(1991)，Ursin(2002)等在关于如何计算吸收反射系数方面给出了一系列理论上的讨论，Innanen和Weglein(2007)将一般的吸收逆散射问题简化为反演单一界面处的反射系数，Innanen(2011)详细推导了地震波由弹性非衰减介质入射到吸收介质时反射及透射系数的表达式。实际上，完整描述或反演地下介质中吸收的反射系数，既需要考虑振幅随角度的变化(AVA)，同时也要考虑振幅随频率的变化关系(AVF)。

目前，国际上尚没有综合考虑了以上两种情况的关于粘弹性介质中球面波反射系数特征的研究。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种粘弹性介质中的球面PP波反射系数建模方法，以实现粘弹性介质中球面波反射系数建模为目标，将White模型和球面波的平面波分解算法有效结合，实现了粘弹性介质中球面PP波反射系数建模。该建模方法物理意义明确且更加忠于实际地震勘探情况，对于更加精确地描述含油气储层有一定的指导意义。

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

一种粘弹性介质中的球面PP波反射系数建模方法，包括：

(1)基于White模型计算纵波相速度υ_p和品质因子Q^-1：

从储层参数出发，基于White模型计算得到纵波相速度υ_p和品质因子Q^-1；所述储层参数包括岩石骨架特性参数和孔隙流体特性参数；

(2)求取频散孔隙介质中平面PP波反射系数；

将第(1)步得到的纵波相速度υ_p和品质因子Q^-1代入到频散介质的Zoeppritz方程中，计算得到频散孔隙介质中的平面PP波的反射系数

(3)求取球面波反射系数；

在第(2)步中求得平面PP波反射系数之后，根据球面波的平面波分解 算法对频散孔隙介质中的球面PP波反射系数进行建模，求得频散孔隙介质中的球面PP波反射系数

所述步骤(1)是利用公式(12)求得纵波相速度υ_p，利用公式(13)求得品质因子Q^-1的：

${\upsilon }_p={\left[\mathrm{Re}\left(\frac{1}{\upsilon }\right)\right]}^{-1},---\left(12\right)$

$Q^{-1}={\arctan }^{-1}\left[\frac{\mathrm{Im}\left({\upsilon }^2\right)}{\mathrm{Re}\left({\upsilon }^2\right)}\right],---\left(13\right)$

其中，Re和Im分别代表实部和虚部；

υ为复速度，根据关系求得；

其中为平均密度，ρ₁和ρ₂分别为含气和含水孔隙介质的密度，p₁和p₂分别为权重因子；p₁=d₁/(d₁+d₂)，p₂=d₂/(d₁+d₂)，每个周期内包含两个孔隙介质层1和2，对应小标1和2，层厚为d_l，l=1,2；

$E={[\frac{1}{E_0}+\frac{2{(r_2-r_1)}^2}{\mathrm{i\omega }(d_1+d_2)(I_1+I_2)}]}^{-1},---\left(1\right)$

其中，

$E_0={[\frac{p_1}{E_{G1}}+\frac{p_2}{E_{G2}}]}^{-1},---\left(2\right)$

其中，$E_G=K_G+\frac{4}{3}{\mu }_{\mathrm{dry}},---\left(3\right)$

式中μ_dry为干岩石的剪切模量，K_G为Gassmann体积模量：

K_G=K_dry+α²M,(4)

其中

$\alpha =1-\frac{K_{\mathrm{dry}}}{K_s},---\left(5\right)$

$M=\frac{K_s}{1-\phi -K_{\mathrm{dry}}/K_s+{\mathrm{\phi K}}_s/K_f},---\left(6\right)$

α为Biot系数，K_s是岩石矿物体积模量，K_f是流体体积模量，K_dry是干岩石体积模量，φ是孔隙度。

$r=\frac{\mathrm{\alpha M}}{E_G},---\left(7\right)$

$I=\frac{\eta }{\mathrm{\kappa k}}\coth \left(\frac{\mathrm{kd}}{2}\right),---\left(8\right)$

其中η为流体的粘滞性系数，κ为岩石骨架的渗透率，

$k=\sqrt{\frac{\mathrm{i\omega \eta }}{{\mathrm{\kappa K}}_E}},---\left(9\right)$

$K_E=\frac{E_{\mathrm{dry}}M}{E_G},---\left(10\right)$

$E_{\mathrm{dry}}=K_{\mathrm{dry}}+\frac{4}{3}{\mu }_{\mathrm{dry}},---\left(11\right).$

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

所述频散介质的Zoeppritz方程如下：

$\left(\begin{array}{cccc}\sin {\theta }_{p1}& \cos {\theta }_{s1}& -\sin {\theta }_{p2}& \cos {\theta }_{s2}\\ \cos {\theta }_{p1}& -\sin {\theta }_{s1}& \cos {\theta }_{p2}& \sin {\theta }_{s2}\\ \frac{{\rho }_1{c}_{s1}^2}{c_{p1}}\sin {2\theta }_{p1}& {\rho }_1{c}_{s1}\cos {2\theta }_{s1}& \frac{{\rho }_2{c}_{s2}^2}{c_{p2}}\sin {2\theta }_{p2}& -{\rho }_2{c}_{s2}\cos {2\theta }_{s2}\\ {\rho }_1{c}_{p1}\cos {2\theta }_{s1}& -{\rho }_1{c}_{s1}\sin {2\theta }_{s1}& -{\rho }_2{c}_{p2}\cos {2\theta }_{s2}& -{\rho }_2{c}_{s2}\sin {2\theta }_{s2}\end{array}\right)$

$\times \left(\begin{array}{c}{R}_{\mathrm{pp}}^{*}\\ R_{\mathrm{ps}}^{*}\\ T_{\mathrm{pp}}^{*}\\ T_{\mathrm{ps}}^{*}\end{array}\right)=\left(\begin{array}{c}-\sin {\theta }_{p1}\\ \cos {\theta }_{p1}\\ \frac{{\rho }_1{c}_{s1}^2}{c_{p1}}\sin {2\theta }_{p1}\\ -{\rho }_1{c}_{p1}\cos {2\theta }_{s1}\end{array}\right),---\left(14\right)$

其中c是速度，θ是角度，脚标p和s分别代表P波和S波，脚标1和2分别代表泥岩盖层和砂岩储层。依赖于频率的P波速度可以表示为：

$\frac{1}{c}=\frac{1}{V}-\frac{\mathrm{i\alpha }}{\omega },---\left(15\right)$

V为相速度，通过公式(12)求得，α为吸收系数，其表达式为：

$\alpha =(\sqrt{Q^2+1}-Q)\frac{\omega }{V},---\left(16\right)$

将步骤(1)得到的纵波相速度υ_p和品质因子Q^-1代入公式(15)和公式(16)，然后求解方程(14)，得到频散孔隙介质中的平面PP波的反射系数是入射角θ_p1和频率ω的函数。

所述步骤(3)中所述根据球面波的平面波分解算法对频散孔隙介质中的球面PP波反射系数进行建模是通过下式实现的：

$R_{\mathrm{pp}}^{\mathrm{spherical}}\left({\theta }_i\right)=\underset{\Gamma }{\int }W\left({S,\theta ,\theta }_i\right)R_{\mathrm{pp}}^{*}(\theta ,\omega )d\left(\cos \theta \right)---\left(19\right)$

其中，

$W\left({S,\theta ,\theta }_i\right)=\frac{[-J_1(\sin \theta \sin {\theta }_i/S)\sin \theta \sin {\theta }_i+i_1{J}_1(\sin \theta \sin {\theta }_i/S)\cos \theta \cos {\theta }_i]}{S(1-\mathrm{iS})\exp [i(1-\cos \theta \cos {\theta }_i)/S]}---\left(20\right)$

其中，S＝V_p1/(Rω)，V_p1为上层介质中的P波速度(Vp是已知的)，R=(z+h)/cosθ_i，其中θ_i是球面波的入射角，θ是平面波的入射角，J₁代表一阶Bessel函数，i是虚数单位；

h为震源到反射界面的深度，z为检波器到反射界面的深度；

通过求解公式(19)即得到频散孔隙介质中的球面PP波反射系数公式(19)中等号右边的所有参数都是已知的。

与现有技术相比，本发明的有益效果是：本发明相对于已有技术最大的创新点就是引入了地层的粘弹性信息。不仅考虑了地震波在地下传播过程中的球面波效应，而且考虑了含油气储层的频散和衰减效应，能够更加精确地反映储层的物理特性。

附图说明

图1 储层模型示意图。

图2 震源和接收点间界面反射射线的几何示意图。

图3-1 P波相速度随频率的变化关系。

图3-2 品质因子随频率的变化关系。

图4-1 平面PP波反射系数振幅随入射角和频率变化。

图4-2 平面PP波反射系数相位随入射角和频率变化。

图5-1 球面PP波反射系数振幅随入射角和频率变化。

图5-2 球面PP波反射系数相位随入射角和频率变化。

图6 本发明方法的步骤框图。

具体实施方式

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

本发明将White模型和球面波的平面波分解算法进行结合，综合考虑了粘弹性介质的衰减和频散效应以及地震波传播的球面波效应。首先根据White模型和粘弹性介质的Zeoppritz方程得到频散孔隙介质中的平面波反射系数，然后基于球面波的平面波分解算法得到频散孔隙介质中球面波反射系数的AVAF(振幅随频率和角度变化)特征。

如图6所示，本方法具体包括：

(1)基于White模型计算纵波相速度υ_p和品质因子Q^-1

从储层参数出发，基于White模型的建模方法，可以得到纵波相速度υ_p和品质因子Q^-1。储层参数包括岩石骨架特性参数和孔隙流体特性参数，所述岩石骨架特性参数包括岩石矿物体积模量K_s，干岩石体积模量K_dry，干岩石的剪切模量μ_dry，孔隙度φ，岩石骨架渗透率κ，所述孔隙流体特性参数包括流体体积模量K_f，流体粘滞性系数η，流体密度ρ。

White模型的具体建模方法如下：

考虑周期性薄互层，每个周期内包含两个孔隙介质层1和2(该方法的使用不受孔隙介质层数的限制)，孔隙内分别含气和含水，层厚为d_l，l=1,2(如图1所示)。根据White模型，Carcione和Picotti给出该模型的P波复模量表达式为：

$E={[\frac{1}{E_0}+\frac{2{(r_2-r_1)}^2}{\mathrm{i\omega }(d_1+d_2)(I_1+I_2)}]}^{-1},---\left(1\right)$

其中，

$E_0={[\frac{p_1}{E_{G1}}+\frac{p_2}{E_{G2}}]}^{-1},---\left(2\right)$

式(2)中p₁=d₁/(d₁+d₂)，l=1,2，忽略脚标，对于每层介质可以得到：

$E_G=K_G+\frac{4}{3}{\mu }_{\mathrm{dry}},---\left(3\right)$

式中μ_dry为干岩石的剪切模量，K_G为Gassmann体积模量：

K_G=K_dry+α²M,(4)

其中

$\alpha =1-\frac{K_{\mathrm{dry}}}{K_s},---\left(5\right)$

$M=\frac{K_s}{1-\phi -K_{\mathrm{dry}}/K_s+{\mathrm{\phi K}}_s/K_f},---\left(6\right)$

α也称Biot系数，K_s是岩石矿物体积模量，K_f是流体体积模量，K_dry是干岩石体积模量，φ是孔隙度。令

$r=\frac{\mathrm{\alpha M}}{E_G},---\left(7\right)$

为快P波流体张力与总正应力的比值，

$I=\frac{\eta }{\mathrm{\kappa k}}\coth \left(\frac{\mathrm{kd}}{2}\right),---\left(8\right)$

是与慢P波相关的阻抗，其中η为流体的粘滞性系数，κ为岩石骨架的渗透率，

$k=\sqrt{\frac{\mathrm{i\omega \eta }}{{\mathrm{\kappa K}}_E}},---\left(9\right)$

是慢P波的复波数，

$K_E=\frac{E_{\mathrm{dry}}M}{E_G},---\left(10\right)$

是有效模量，其中

$E_{\mathrm{dry}}=K_{\mathrm{dry}}+\frac{4}{3}{\mu }_{\mathrm{dry}},---\left(11\right)$

若令ρ为体积密度，则根据关系可得到复速度υ，其中(P是权重因子，具体意义在公式(2)之后给予了说明，ρ₁和ρ₂分别为含气和含水孔隙介质的密度)是平均密度。这样，相速度和品质因子可表示为：

${\upsilon }_p={\left[\mathrm{Re}\left(\frac{1}{\upsilon }\right)\right]}^{-1},---\left(12\right)$

$Q^{-1}={\mathrm{atc}\tan }^{-1}\left[\frac{\mathrm{Im}\left({\upsilon }^2\right)}{\mathrm{Re}\left({\upsilon }^2\right)}\right],---\left(13\right)$

Re和Im分别代表实部和虚部。

(2)频散孔隙介质中平面PP波反射系数求取；

将第(1)步得到的纵波相速度υ_p和品质因子Q^-1代入到频散介质的Zoeppritz方程中，可以得到频散孔隙介质中的平面PP波的反射系数

在White模型中，剪切波的吸收特性不受含气性影响，因此在大多数情况下可以忽略。由Hooke定律和界面处的连续性条件，可以得到频散介质对应的Zoeppritz方程：

$\left(\begin{array}{cccc}\sin {\theta }_{p1}& \cos {\theta }_{s1}& -\sin {\theta }_{p2}& \cos {\theta }_{s2}\\ \cos {\theta }_{p1}& -\sin {\theta }_{s1}& \cos {\theta }_{p2}& \sin {\theta }_{s2}\\ \frac{{\rho }_1{c}_{s1}^2}{c_{p1}}\sin {2\theta }_{p1}& {\rho }_1{c}_{s1}\cos {2\theta }_{s1}& \frac{{\rho }_2{c}_{s2}^2}{c_{p2}}\sin {2\theta }_{p2}& -{\rho }_2{c}_{s2}\cos {2\theta }_{s2}\\ {\rho }_1{c}_{p1}\cos {2\theta }_{s1}& -{\rho }_1{c}_{s1}\sin {2\theta }_{s1}& -{\rho }_2{c}_{p2}\cos {2\theta }_{s2}& -{\rho }_2{c}_{s2}\sin {2\theta }_{s2}\end{array}\right)$

$\times \left(\begin{array}{c}{R}_{\mathrm{pp}}^{*}\\ R_{\mathrm{ps}}^{*}\\ T_{\mathrm{pp}}^{*}\\ T_{\mathrm{ps}}^{*}\end{array}\right)=\left(\begin{array}{c}-\sin {\theta }_{p1}\\ \cos {\theta }_{p1}\\ \frac{{\rho }_1{c}_{s1}^2}{c_{p1}}\sin {2\theta }_{p1}\\ -{\rho }_1{c}_{p1}\cos {2\theta }_{s1}\end{array}\right),---\left(14\right)$

其中c是速度，θ是角度，脚标p和s分别代表P波和S波，脚标1和2分别代表泥岩盖层和砂岩储层。依赖于频率的P波速度可以表示为：

$\frac{1}{c}=\frac{1}{V}-\frac{\mathrm{i\alpha }}{\omega },---\left(15\right)$

V为相速度(就是公式(12)中的υ_p：)，α为吸收系数，其表达式为

$\alpha =(\sqrt{Q^2+1}-Q)\frac{\omega }{V},---\left(16\right)$

P波的相速度和品质因子可以通过式(12)和式(13)计算得到。

通过求解方程(14)，可以得到频散孔隙介质中的平面PP波反射系数是入射角θ_p1和频率ω的函数。

(3)基于平面波分解算法的球面波反射系数求取；

在第(2)步中求得频散孔隙介质中平面PP波反射系数之后，根据球面波的平面波分解算法，最终就可以建模得到频散孔隙介质中的球面PP波反射系数具体方法如下：

根据球面波的平面波分解理论，反射球面谐波的位移势函数可以表示为 (Haase，2004)：

$\phi \left(\omega \right)=\mathrm{Ai\omega exp}(-\mathrm{i\omega t})\underset{0}{\overset{\infty }{\int }}\frac{p}{\xi }{R}_{\mathrm{pp}}\left(p\right)J_0\left(\mathrm{\omega pr}\right)\exp \left[\mathrm{i\omega \xi }(z+h)\right]\mathrm{dp}---\left(17\right)$

其中r，z，h的物理意义如图2所示，r为偏移距大小，h为震源到反射界面的深度，z为检波器到反射界面的深度，ω为球面谐波的角频率，p和ξ分别为上层介质中的水平和垂直慢度，J_o为零阶贝塞尔函数，R_pp为平面波的反射系数。在粘弹性介质中，平面波的反射系数R_pp不仅与入射角度有关，还与频率有关，所以将通过解方程(14)得到的代入式(17)：

$\phi \left(\omega \right)=\mathrm{Ai\omega exp}(-\mathrm{i\omega t})\underset{0}{\overset{\infty }{\int }}\frac{p}{\xi }{R}_{\mathrm{pp}}^{*}(p,\omega )J_0\left(\mathrm{\omega pr}\right)\exp \left[\mathrm{i\omega \xi }(z+h)\right]\mathrm{dp}---\left(18\right)$

式(18)即为粘弹性介质中球面谐波势函数的表达式。为得到球面谐波的反射系数形式，对式(18)做如下处理：

①在检波点位置处求位移势函数在波传播方向的偏导数，得到位移函数；

②用均匀介质中的球面波场对反射波场校正，消除球面扩散对反射波振幅的影响；

③将积分变量由射线参数p变为cosθ。

最终得到频散孔隙介质中球面谐波反射系数的形式为：

$R_{\mathrm{pp}}^{\mathrm{spherical}}\left({\theta }_i\right)=\underset{\Gamma }{\int }W\left({S,\theta ,\theta }_i\right)R_{\mathrm{pp}}^{*}(\theta ,\omega )d\left(\cos \theta \right)---\left(19\right)$

其中，

$W\left({S,\theta ,\theta }_i\right)=\frac{[-J_1(\sin \theta \sin {\theta }_i/S)\sin \theta \sin {\theta }_i+i_1{J}_1(\sin \theta \sin {\theta }_i/S)\cos \theta \cos {\theta }_i]}{S(1-\mathrm{iS})\exp [i(1-\cos \theta \cos {\theta }_i)/S]}---\left(20\right)$

上式中S＝V_p1/(Rω)，V_p1为上层介质中的P波速度，R=(z+h)/cosθ_i。

(19)式就是本发明方法最后得到的表达式，通过(19)式就完成了本发明的建模。

下面通过一个理论储层模型来说明本发明的效果。模型为具有低阻抗的浅层疏松含气砂岩储层，上覆有泥岩盖层。泥岩盖层的特性参数为：Vp=2190m/s，Vs=820m/s，ρ=2.16g/cm³。砂岩储层的特性参数为：K_dry＝1.56GPa，μ_dry＝1.10GPa，φ=0.33，K＝2.00darcy，ρ_g＝0.15g/cm³，η_g＝0.01cP，K_w=2.42GPa，ρ_w=1.00g/cm³，η_w＝1.00cP，d1=2m，d2=0.5m，h=500m。砂岩岩石颗粒的体积模量和密度分别为K_s=38GPa和ρ_s＝2.65g/cm³。

图3-1和图3-2分别为由White模型得到的砂岩储层的P波相速度和品质因子随频率的变化关系。图4-1和图4-2分别为频散孔隙介质中平面PP波反射系数AVAF特征，即砂岩储层定界面对应的平面PP波反射系数(振幅和相位)随入射角和频率的变化关系。图5-1和图5-2为使用本发明的方法得到的砂岩储层定界面对应的球面PP波反射系数(振幅和相位)随入射角和频率的变化关系，即频散孔隙介质中球面PP波反射系数AVAF特征。

AVO技术作为一种较为成熟的烃类直接检测手段，在实际生产应用中已经取得了良好的效果。传统的AVO技术是基于平面波理论建立起来的，并总是假设地下介质是完全弹性的。然而，地震勘探中使用的是爆炸源，得到的是球面波记录，当震源或检波器离反射界面很近或者需要研究临界角附近的反射波场的时候，平面波近似就变得不准确。而且VSP资料、测井资料以及实验室中的岩石物理测量均显示地震波在实际传播过程中会发生衰减以及速度频散现象，尤其对于含有碳氢化合物的区域，衰减非常明显。所以，忽略地震波在地下传播的球面波效应及地下介质的频散和衰减效应会对AVO分析及储层预测带来巨大的风险。

本发明将White模型和球面波的平面波分解算法有效结合，综合考虑了地震波传播的球面波效应及介质的频散和衰减效应，得到了一种粘弹性介质中球面PP波反射系数的建模方法，并对其AVAF特征进行了研究。这项发明物理意义明确且更加忠于实际地震勘探情况，能够更加真实地反映含油气储层的物理 特性，对于进行精确地储层预测工作具有一定的指导意义。

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