技术领域
[0001] 本
发明涉及地质技术领域，尤其涉及一种基于岩石物理模型的含气砂岩储层地震响应数值模拟方法。
背景技术
[0002] 目前基于岩石物理模型的
地震波衰减特征分析多停留在对反射系数讨论，而传统
波动方程正演模拟中模型的建立多是基于对实际
地层中地震
波速度的估计，因此由于增加了人为干预的主观性变得不准确，基于岩石物理模型建立的地质模型更加符合实际，因此岩石物理模型与波动方程正演的结合是数值模拟方法必然的发展趋势。
[0003] 斑
块饱和模型是一种饱含
流体，且部分含气部分含
水的岩石物理模型，它能够用来模拟常规含气砂岩的岩石物理特征；黏滞弥散波动方程主要用来模拟含流体介质的地震响应特征。
发明内容
[0004] 本发明的目的在于解决上述
现有技术存在的
缺陷，提供一种基于岩石物理模型的含气砂岩储层地震响应数值模拟方法，能够有效地将岩石物理模型的地球物理特征表现出来，更加直观且准确地给出常规含气砂岩的地震响应特征。
[0005] 未达到上述目的，本发明采用如下方案：
[0006] 基于岩石物理模型的含气砂岩储层地震响应数值模拟方法，岩石物理模型为空间周期排列球状斑块饱和模型；该方法包括以下步骤：
[0007] 根据Wood定律计算低频情况下孔隙流体的有效体积模量，从而利用Gassmann方程得到岩石等效体积模量；高频下首先根据Gassmann方程计算不同流体的体积模量，再利用Hill理论计算岩石体积模量；利用Johnson理论给出动态体积模量，进而计算地震纵波复速度；根据Carcoine衰减理论得到地震P波相速度和品质因子的频散关系，并将其引入
频率波数域得到复波数，结合利用频变反射系数公式与子波褶积得到的初始波场，代入二维黏滞弥散波动方程进行数值模拟。
[0008] 进一步地，具体实现过程如下：
[0009] 基于空间周期排列球状斑块饱和模型，球体内部饱含气，外部饱含水，当地震波传播到斑块饱和介质时，由于气和水的存在导致产生了不同的孔隙压
力，利用Biot理论解释相应的孔隙压力动态平衡过程，并给出相应的地震波扩散长度：
[0010]
[0011] 其中ω为
角频率，Fi为扩散系数，i＝g,w；其中g代表气，w代表水：
[0012]
[0013] 其中，κ为斑块饱和模型的渗透率，η为孔隙流体的黏滞系数， Mc、Kav根据含流体多孔介质的一些物理属性给出：
[0014]
[0015]
[0016]
[0017] 式中，Km、Ks、Kf分别为干岩石、固体颗粒和孔隙流体的体积模量，μ为含流体多孔介质的
剪切模量，根据Hill理论，假设其余各混合组分的剪切模量都相同，φ为孔隙度，KG(Kf)则根据对应流体体积模量和Gassmann方程算出：
[0018]
[0019] 当地震波频率足够低时，流体之间有足够的时间来达到孔隙压力平衡，基于Wood定律可以得到孔隙流体的有效体积模量：
[0020]
[0021] 其中fi(i＝g,w)为介质孔隙中流体的体积分量，这里以气和水的
饱和度代替，Ki(i＝g,w)为对应气和水的体积模量；于是在低频背景下，根据Gassmann方程计算出岩石的等效体积模量：
[0022]
[0023] 相对的，当地震波频率很高时，流体之间没有时间达到孔隙压力平衡，此时孔隙中产生的压力是不均匀的，若假设孔隙压力不同但均为常数，可根据Gassmann方程首先计算不同流体的体积模量，再根据Hill理论计算高频下的岩石等效体积模量：
[0024]
[0025] 根据高频和低频下的岩石等效体积模量，求得对应中间频率时含气、水岩石的等效体积模量：
[0026]
[0027] 式中，i为虚数单位，θ和σ可分别根据干岩石、流体的物理性质等求出，这里给出相应的计算公式：
[0028]
[0029]
[0030]
[0031] Zi＝φ2Kav(Ki) (14)
[0032]
[0033]
[0034] 上式中，Di(i＝g,w)即为对应流体的扩散系数，可由之前提到的Biot理论计算得到，在空间周期排列斑块饱和模型中，S/V、T可表示为：
[0035]
[0036]
[0037] 其中：
[0038]
[0039] 通过上述公式计算得到地震纵波在含气、含水储层中传播的复速度：
[0040]
[0041] 其中，ρb为介质的体
密度， 式中ρs、ρi分别为固体颗粒和流体的密度，Si为对应的孔隙流体的饱和度；根据Carcoine衰减理论，利用复速度给出地震波在衰减介质中的等效相速度和品质因子：
[0042]
[0043]
[0044] 依此得到斑块饱和介质中地震波相速度和品质因子的频散关系，利用上述频散关系，结合二维黏滞-弥散波动方程，进行数值模拟。
[0045] 进一步地，黏滞-弥散波动方程正演的具体实现步骤如下：
[0046] 首先根据二维黏滞-弥散波动方程：
[0047]
[0048] 上式中，u为质点位移，γ为含流体介质的弥散系数，η为含流体介质的黏滞系数，v为介质中地震波速度，给出方程的简谐
波形式的解：
[0049]
[0050] 式中，ω为圆频率，kx、kz分别为x、z方向的波数，单位为1/m,i为虚数单位。将上式代入波动方程有：
[0051]
[0052] 于是有：
[0053]
[0054] 令kz＝k+iα(27)
[0055] 等式两边平方得到 得k和α关于黏滞弥散系数的表达式；
[0056] 根据Gazdag等的频率-波数域波场延拓公式：
[0057]
[0058] 其中，ei(k+iα)dz为
相移因子，利用上式进行叠前黏滞-弥散波动方程正演数值模拟；
[0059] 根据斑块饱和模型中地震波相速度的频散关系式(21)，将其代入计算波数的公式式(26-28)求得相移因子，设计模型对含气砂岩的地震响应特征进行较准确的分析。
[0060] 本发明基于球状斑块饱和模型，能够较好地描述常规含气砂岩的岩石物理和地震动力学特征，应用黏滞弥散波动方程理论，引入地震波在传播时因流体的存在而造成的衰减和
能量损失。利用频率-波数域黏滞弥散波动方程进行正演模拟，能够有效地将岩石物理模型的地球物理特征表现出来，结合黏滞弥散理论带来的含流体介质的黏滞特征和弥散特征，能够更加直观且准确地给出常规含气砂岩的地震响应特征。
附图说明
[0061] 图1地质模型
[0062] 图2不同含气饱和度下地震波相速度与逆品质因子的频散关系曲线
[0063] 图3不同含气饱和度下三类砂岩储层地震响应
[0064] 图4不同含气饱和度时下界面反射波
频谱具体实施方式
[0065] 为使本发明的目的、技术方案和优点更加清楚，下面本发明中的技术方案进行清楚、完整地描述，显然，所描述的
实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0066] 本发明提供了一种基于岩石物理模型的含气砂岩储层地震响应数值模拟方法，该模型为空间周期排列球状斑块饱和模型，该方法包括以下步骤：根据Wood定律计算低频情况下孔隙流体的有效体积模量，从而利用Gassmann方程得到岩石等效体积模量；高频下首先根据Gassmann方程计算不同流体的体积模量，再利用Hill理论计算岩石体积模量；利用Johnson理论给出动态体积模量，进而计算地震纵波复速度；根据Carcoine衰减理论得到地震P波相速度和品质因子的频散关系，并将其引入频率波数域得到复波数，结合利用频变反射系数公式与子波褶积得到的初始波场，代入二维黏滞弥散波动方程即可进行数值模拟。本发明给出了一种更加贴合实际的含气砂岩储层数值模拟方法。以下为具体实现过程：
[0067] 基于空间周期排列球状斑块饱和模型，球体内部饱含气，外部饱含水。当地震波传播到斑块饱和介质时，由于气和水的存在导致产生了不同的孔隙压力，利用Biot理论解释相应的孔隙压力动态平衡过程，并给出相应的地震波扩散长度：
[0068]
[0069] 其中ω为角频率，Fi(i＝g,w；其中g代表气，w代表水)为扩散系数：
[0070]
[0071] 其中，κ为斑块饱和模型的渗透率，η为孔隙流体的黏滞系数， Mc、Kav则可根据含流体多孔介质的一些物理属性给出：
[0072]
[0073]
[0074]
[0075] 式中，Km、Ks、Kf分别为干岩石、固体颗粒和孔隙流体的体积模量，μ为含流体多孔介质的剪切模量，根据Hill理论，我们假设其余各混合组分的剪切模量都相同，φ为孔隙度，KG(Kf)则可根据对应流体体积模量和Gassmann方程算出：
[0076]
[0077] 当地震波频率足够低时，流体之间有足够的时间来达到孔隙压力平衡，基于Wood定律可以得到孔隙流体的有效体积模量：
[0078]
[0079] 其中fi(i＝g,w)为介质孔隙中流体的体积分量，这里以气和水的饱和度代替，Ki(i＝g,w)为对应气和水的体积模量。于是在低频背景下，根据Gassmann方程可以计算出岩石的等效体积模量：
[0080]
[0081] 相对的，当地震波频率很高时，流体之间没有时间达到孔隙压力平衡，此时孔隙中产生的压力是不均匀的，若假设孔隙压力不同但均为常数，可根据Gassmann方程首先计算不同流体的体积模量，再根据Hill理论计算高频下的岩石等效体积模量：
[0082]
[0083] 根据高频和低频下的岩石等效体积模量，可以求得对应中间频率时含气、水岩石的等效体积模量：
[0084]
[0085] 式中，i为虚数单位，θ和σ可分别根据干岩石、流体的物理性质等求出，这里给出相应的计算公式：
[0086]
[0087]
[0088]
[0089] Zi＝φ2Kav(Ki) (14)
[0090]
[0091]
[0092] 上式中，Di(i＝g,w)即为对应流体的扩散系数，可由之前提到的Biot理论计算得到。在空间周期排列斑块饱和模型中，S/V、T可表示为：
[0093]
[0094]
[0095] 其中：
[0096]
[0097] 通过上述公式计算得到地震纵波在含气、含水储层中传播的复速度：
[0098]
[0099] 其中，ρb为介质的体密度， 式中ρs、ρi分别为固体颗粒和流体的密度，Si为对应的孔隙流体的饱和度。根据Carcoine衰减理论，利用复速度给出地震波在衰减介质中的等效相速度和品质因子：
[0100]
[0101]
[0102] 依此即可得到斑块饱和介质中地震波相速度和品质因子的频散关系。
[0103] 利用上述频散关系，结合二维黏滞-弥散波动方程，即可进行数值模拟。下面将介绍黏滞-弥散波动方程正演的具体实现过程。
[0104] 首先根据二维黏滞-弥散波动方程：
[0105]
[0106] 上式中，u为质点位移，γ为含流体介质的弥散系数，η为含流体介质的黏滞系数，v为介质中地震波速度。给出方程的简谐波形式的解：
[0107]
[0108] 式中，ω为圆频率，kx、kz分别为x、z方向的波数，单位为1/m,i为虚数单位。将上式代入波动方程有：
[0109]
[0110] 于是有：
[0111]
[0112] 令kz＝k+iα(27)
[0113] 等式两边平方得到 根据前人研究，可省去高阶α2项。于是可得k和α关于黏滞弥散系数的表达式。
[0114] 根据Gazdag等的频率-波数域波场延拓公式可有：
[0115]
[0116] 其中，ei(k+iα)dz为相移因子，利用上式即可进行叠前黏滞-弥散波动方程正演数值模拟。
[0117] 根据斑块饱和模型中地震波相速度的频散关系(式21)，将其代入上述计算波数的公式(式26-28)求得相移因子，设计模型即可对含气砂岩的地震响应特征进行较准确的分析。
[0118] 附图1为设计的三层地质模型，其中上层和底层为
页岩，中间为含气砂岩储层。设计一个欠
压实的砂岩，给定其渗透率、孔隙度等物性参数，其中岩石颗粒的体积模量均为38GPa，密度均为2.65g/cm3，干岩石的体积模量和剪切模量分别为1.56GPa和1.10GPa，孔隙度为33％，孔隙中所含气的体积模量为0.03GPa，密度为0.15g/cm3，水的体积模量为
2.42GPa，密度为1.0g/cm3。
[0119] 利用所述技术路线，首先可得到如图2所示基于斑块饱和模型的欠压实砂岩储层中地震P波相速度和逆品质因子的频散关系曲线。从图中可以看出，当含气饱和度增大，介质中流体对地震波的衰减产生了巨大的影响，导致其速度发生了数百米每秒的降低，且在低频段发生了较明显的频散现象。而逆品质因子的频散曲线则表明了，含气饱和度增大能够使逆品质因子的峰值向高频方向移动，峰值频率时地震波速度衰减最强。
[0120] 图2不同含气饱和度下地震波相速度与逆品质因子的频散关系曲线
[0121] 进而根据二维黏滞-弥散波动方程，代入斑块饱和模型的速度频散关系，并基于图1所示的地质模型，可以得到如图3所示的第128道叠前数值模拟结果，从图中可以明显地看到由于含气量的增加，流体对地震波的衰减增强，导致储层下界面反射出现明显的时间延迟和振幅变化，其中振幅受储层波阻抗变化的影响而增大。
[0122] 提取储层下反射界面地震波的频谱，如附图4所示，可以看到地震波的主频随着储层中含流体量增加逐渐向低频移动，反映出地震波在高频段的衰减变强。
[0123] 根据上述方法，基于球状斑块饱和模型，利用黏滞-弥散波动方程正演，能够有效地反映经典含气砂岩中地震波因流体的存在而发生的衰减和损失，岩石物理模型与地震波动方程正演的结合能够提供更加准确的地震模型正演方法，这能够为进一步理解含流体层的地震响应提供一定理论支持。
[0124] 最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行
修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。
