基于岩石物理模型的含气砂岩储层地震响应数值模拟方法与流程

技术特征：

1.基于岩石物理模型的含气砂岩储层地震响应数值模拟方法，其特征在于：岩石物理模型为空间周期排列球状斑块饱和模型；

该方法包括以下步骤：

根据Wood定律计算低频情况下孔隙流体的有效体积模量，从而利用Gassmann方程得到岩石等效体积模量；高频下首先根据Gassmann方程计算不同流体的体积模量，再利用Hill理论计算岩石体积模量；利用Johnson理论给出动态体积模量，进而计算地震纵波复速度；根据Carcoine衰减理论得到地震P波相速度和品质因子的频散关系，并将其引入频率波数域得到复波数，结合利用频变反射系数公式与子波褶积得到的初始波场，代入二维黏滞弥散波动方程进行数值模拟。

2.根据权利要求1所述的基于岩石物理模型的含气砂岩储层地震响应数值模拟方法，其特征在于具体实现过程如下：

基于空间周期排列球状斑块饱和模型，球体内部饱含气，外部饱含水，当地震波传播到斑块饱和介质时，由于气和水的存在导致产生了不同的孔隙压力，利用Biot理论解释相应的孔隙压力动态平衡过程，并给出相应的地震波扩散长度：

$L=\sqrt{F_i/\mathrm{\ω}}---\left(1\right)$

其中ω为角频率，F_i为扩散系数，i＝g,w；其中g代表气，w代表水：

其中，κ为斑块饱和模型的渗透率，η为孔隙流体的黏滞系数，M_c、K_av根据含流体多孔介质的一些物理属性给出：

$M_c\left(K_f\right)=K_G\left(K_f\right)+\frac{4}{3}\mathrm{\μ}---\left(5\right)$

式中，K_m、K_s、K_f分别为干岩石、固体颗粒和孔隙流体的体积模量，μ为含流体多孔介质的剪切模量，根据Hill理论，假设其余各混合组分的剪切模量都相同，φ为孔隙度，K_G(K_f)则根据对应流体体积模量和Gassmann方程算出：

当地震波频率足够低时，流体之间有足够的时间来达到孔隙压力平衡，基于Wood定律可以得到孔隙流体的有效体积模量：

$K_R={\left(\frac{f_i}{K_i}\right)}^{-1}---\left(7\right)$

其中f_i(i＝g,w)为介质孔隙中流体的体积分量，这里以气和水的饱和度代替，K_i(i＝g,w)为对应气和水的体积模量；于是在低频背景下，根据Gassmann方程计算出岩石的等效体积模量：

相对的，当地震波频率很高时，流体之间没有时间达到孔隙压力平衡，此时孔隙中产生的压力是不均匀的，若假设孔隙压力不同但均为常数，可根据Gassmann方程首先计算不同流体的体积模量，再根据Hill理论计算高频下的岩石等效体积模量：

$K_H=\underset{i=g,w}{\Σ}\frac{M_c\left(K_i\right)}{S_i}-\frac{4}{3}\mathrm{\μ}---\left(9\right)$

根据高频和低频下的岩石等效体积模量，求得对应中间频率时含气、水岩石的等效体积模量：

$K_D=K_H-\frac{K_H-K_L}{1-\mathrm{\θ}+\mathrm{\θ}\sqrt{1+100i\mathrm{\ω}\mathrm{\σ}/{\mathrm{\θ}}^2}}---\left(10\right)$

式中，i为虚数单位，θ和σ可分别根据干岩石、流体的物理性质等求出，这里给出相应的计算公式：

$\mathrm{\θ}=\frac{K_H-K_L}{2K_L}\·\frac{\mathrm{\σ}}{T}---\left(11\right)$

$\mathrm{\σ}=\frac{{(K_H-K_L)}^2}{K_H^2{G}^2}---\left(12\right)$

$G={\[\frac{\underset{i=g,w}{\Σ}(Z_i+Q_i)M_c\left(K_i\right)}{\underset{i=g,w}{\Σ}{\mathrm{\φS}}_i{K}_G\left(K_i\right)M_c\left(K_i\right)}\]}^2\·\frac{S}{V}\sqrt{F}---\left(13\right)$

Z_i＝φ²K_av(K_i) (14)

$F={\left(\frac{{\mathrm{\κK}}_H}{\underset{i=g,w}{\Σ}{\mathrm{\η}}_i\sqrt{D_i}}\right)}^2---\left(16\right)$

上式中，D_i(i＝g,w)即为对应流体的扩散系数，可由之前提到的Biot理论计算得到，在空间周期排列斑块饱和模型中，S/V、T可表示为：

$\frac{S}{V}=3\frac{R_g^2}{R_w^2}---\left(17\right)$

$T=\frac{1}{\mathrm{\κ}}\frac{K_L{\mathrm{\φ}}^2}{30R_w^3}\left\{\begin{array}{c}\[3{\mathrm{\η}}_w{g}_w^2+5({\mathrm{\η}}_g-{\mathrm{\η}}_w)g_g{g}_w-3{\mathrm{\η}}_g{g}_g^2\]R_g^5\\ +5g_w\[3{\mathrm{\η}}_w{g}_w-(2{\mathrm{\η}}_w+{\mathrm{\η}}_g)g_g\]R_g^2{R}_w^3\\ -15{\mathrm{\η}}_w{g}_w(g_w-g_g)R_g^3{R}_w^2-3{\mathrm{\η}}_w{g}_w^2{R}_w^5\end{array}\right\}---\left(18\right)$

其中：

通过上述公式计算得到地震纵波在含气、含水储层中传播的复速度：

$V_R=\sqrt{\frac{K_D+\frac{4}{3}\mathrm{\μ}}{{\mathrm{\ρ}}_b}}---\left(20\right)$

其中，ρ_b为介质的体密度，式中ρ_s、ρ_i分别为固体颗粒和流体的密度，S_i为对应的孔隙流体的饱和度；根据Carcoine衰减理论，利用复速度给出地震波在衰减介质中的等效相速度和品质因子：

$V_P={\[\mathrm{Re}\left(\frac{1}{V_R}\right)\]}^{-1}---\left(21\right)$

$Q=\frac{\mathrm{Re}\left(V_R^2\right)}{\mathrm{Im}\left(V_R^2\right)}---\left(22\right)$

依此得到斑块饱和介质中地震波相速度和品质因子的频散关系，利用上述频散关系，结合二维黏滞-弥散波动方程，进行数值模拟。

3.根据权利要求2所述的基于岩石物理模型的含气砂岩储层地震响应数值模拟方法，其特征在于：黏滞-弥散波动方程正演的具体实现步骤如下：

首先根据二维黏滞-弥散波动方程：

$\frac{{\∂}^{2}u}{\∂t^2}+\mathrm{\γ}\frac{\∂u}{\∂t}-\mathrm{\η}(\frac{{\∂}^{3}u}{\∂x^2\∂t}+\frac{{\∂}^{3}u}{\∂z^2\∂t})-v^2(\frac{{\∂}^{2}u}{\∂x^2}+\frac{{\∂}^{2}u}{\∂z^2})=0---\left(23\right)$

上式中，u为质点位移，γ为含流体介质的弥散系数，η为含流体介质的黏滞系数，v为介质中地震波速度，给出方程的简谐波形式的解：

$u=e^{i(k_{z}z+k_{x}x)}{e}^{i\mathrm{\ω}t}---\left(24\right)$

式中，ω为圆频率，k_x、k_z分别为x、z方向的波数，单位为1/m,i为虚数单位。将上式代入波动方程有：

于是有：

$k_z^2=\frac{{\mathrm{\ω}}^2(v^2-\mathrm{\γ}\mathrm{\η})}{v^4+{\mathrm{\ω}}^2{\mathrm{\η}}^2}-k_x^2+i\frac{\mathrm{\ω}({\mathrm{\γv}}^2+{\mathrm{\η\ω}}^2)}{v^4+{\mathrm{\ω}}^2{\mathrm{\η}}^2}---\left(26\right)$

令k_z＝k+iα(27)

等式两边平方得到得k和α关于黏滞弥散系数的表达式；

根据Gazdag等的频率-波数域波场延拓公式：

$P(z_0+dz)=P\left(z_0\right)e^{{\mathrm{ik}}_{z}dz}=P\left(z_0\right)e^{i(k+i\mathrm{\α})dz}---\left(28\right)$

其中，e^i(k+iα)dz为相移因子，利用上式进行叠前黏滞-弥散波动方程正演数值模拟；

根据斑块饱和模型中地震波相速度的频散关系式(21)，将其代入计算波数的公式式(26-28)求得相移因子，设计模型对含气砂岩的地震响应特征进行较准确的分析。