技术领域
[0001] 本
发明是关于地球物理勘探中的地震
岩石物理模型及
地震波数值模拟技术，特别是关于一种非饱和双重孔隙介质地震波频散衰减分析方法及装置。
背景技术
[0002] 上
地壳沉积岩在形成过程中普遍具有一定的孔隙度，随着距地表深度的增加，在
压实、固结与充填作用下，岩石的孔隙度会显著降低，但是，局部高孔隙度发育的砂层仍然会为气态
烃的储集创造条件。在含
水层以下，岩石孔隙中同时饱和气、水两相
流体的情况非常普遍，因此，有必要深入探索非饱和岩石中的弹性、
电阻率、
磁性等物理性质，为实际的地球探测与矿产勘查工作提供依据。近年来，中国
天然气供应的需求急剧增加，为满足人民生产生活的需要，生产部
门为天然气田的勘测与开发提出了更高的技术要求，为达到提高气层检测
精度、逐步实现定量识别的目的，开展深入研究以探索含气、水两相流体非饱和岩石中的地震波传播与响应规律的工作势在必行。
[0003] 早期针对含流体多孔介质中的弹性波传播规律的研究，在理论上主要基于宏观均匀性假设，近似认为流体与固体均匀分布于岩石内部的每一个
力学微元中，每一个力学单元内部仅含有一个流体相与一个固体相，而介质内部的所有力学微元都是相同的、均一的。这样的做法明显忽略了岩石内部的局部非均匀性，在实际的应用中也被发现明显低估了地震频段内(几十到几百Hz)岩石内部的
波速频散与
能量衰减现象。
[0004] 基于宏观均一的Biot理论，White等人在1975年首次引入了岩石内部小尺度非均匀性的概念，考虑了含水岩石内部局部分布的气泡对地震波传播的影响。在实际岩石中，局部气体集中形成的非浸润状态的气泡在地震波的
挤压作用下会发生耗散性的振动，即局部流体流动效应。在低频端，局部流体流动的充分发生会造成地震波能量的大量损失，并使得岩石变得松弛，而在高频地震波的激励下，由于流体在单个振荡周期内不能完成完整的局部荡动，岩石会呈现出比较“硬”的等效状态。Dutta和Sheriff在1979年对White理论进行了改进，改善了White理论的低频预测结果与Gassmann理论零频极限的吻合状况。Shapiro和Muller在1999年调查了局部流体流动中的流体压力与岩石参数的关系，得出中观尺度非均匀性诱发的局部流体流动是导致低频地震波能量强衰减的主因结论，2006年，Carcione等的研究显示，若渗透率下降、流体黏性上升或
地层厚度增加，则地震波的衰减与频散曲线有向低频段移动的趋势。
[0005] 聂建新等人在2004年基于同时包含Biot摩擦机制与喷射流机制的非饱和多孔隙BISQ模型，利用小生境遗传
算法实现了储层参数(孔隙度、渗透率、含气
饱和度等)的反演。刘炯等人在2009年用孔隙介质力学的方法，研究了嵌入有球状气泡的含水岩石模型中的地震波传播规律。聂建新、杨顶辉与巴晶在2010年借鉴等效介质的思想，将含水饱和度引入
波动力学控制方程，并考虑了不同波
频率下孔隙流体分布模式对其等效体积模量的影响，给出了能处理含非饱和流体孔隙介质中波传播问题的黏弹性BISQ模型。目前，研究者们在相关问题上已基本达到一致的认识，即局部流体流动是导致孔隙岩石中地震波速度频散与能力衰减的决定性因素，然而，至今为止，文献中还没有出现一个能够普遍被各种研究与产业部门接受的描述局域流的
基础理论与定量技术。
[0006] 综上所述，
现有技术的研究中主要存在以下问题：
[0007] (1)基于黏弹性理论与BISQ理论的相关研究中，由于引入了一些并不具有明确物理意义或并不易于直接物理实现的参数或系数(如各类松弛时间、黏弹性系数或特征喷射流长度等)，使得这些理论虽然在进行岩石波速现象的描述与预测中非常有效，但相关理论与核心参数的数学基础与物理内涵却难以得到实现或验证，因此即使抛开实验验证，这类方法在实际工程中很难得到推广应用。
[0008] (2)在基于精细模型的数值模拟与地震响应分析的相关研究中，由于需要建立过于详细的岩石模型并设置边界条件，使得这类方法虽然能够较准确的给出实际岩石中的地震波响应，但在实际工程中却往往由于缺乏足够的先验信息，并且由于相关方法本身建模的高复杂性与算法的高消耗性，使得此类方法也无法在实际工程应用中得到很好的实现。
[0009] (3)某些地震岩石模型在动力学方程组的推导过程中采用了统计
热力学关系，使得方程的一些核心参数在计算过程中仅能通过频率域的关系式进行上、下限的估算，而无法做到更准确的地震响应的定量预测。
发明内容
[0010] 本发明提供一种非饱和双重孔隙介质地震波频散衰减分析方法及系统，以分析非饱和岩石中的纵横波频散与衰减随频率、孔隙度、含气饱和度的变化规律。
[0011] 为了实现上述目的，本发明提供一种非饱和双重孔隙介质地震波频散衰减分析方法，该方法包括：步骤1：获取包括渗透率、孔隙率、泥质含量及矿物成分的地质数据，并根据所述地质数据生成储层岩石干骨架模型；步骤2：获取包括钻井数据、录井数据、
测井信息及流体实验测量数据的测量数据，根据所述测量数据及Van der Waals方程生成储
层流体模型；步骤3：根据所述储层岩石干骨架模型、储层流体模型及描述非饱和双重孔隙介质地震波传播的Biot-Rayleigh方程求解平面波解，生成纵波、横波的相速度及逆品质因子；步骤4：根据所述纵波、横波的相速度及逆品质因子生成速度、衰减、波阻抗、AVO响应特征等与频率之间的关系。
[0012] 进一步地，在所述步骤1中，生成储层岩石干骨架模型的步骤包括采用Voigt-Reuss-Hill平均模型计算骨架基质的等效
弹性模量，采用 包体模型计算干燥岩石的等效体积模量和
剪切模量：
[0013]
[0014]
[0015] —— 理论求得的等效体积模量和剪切模量；
[0016] Km，μm——骨架基质的体积模量和剪切模量；
[0017] Ki，μi——第i个包体的体积模量和剪切模量；
[0018] xi——每种包体的体积比率；
[0019] N——包体的个数；
[0020] ζm＝μm(9Km+8μm)/6(Km+2μm)；
[0021] Pmi，Qmi——骨架基质m中加入包体材料i后的弹性模量。
[0022] 进一步地，在所述步骤2中，生成储层流体模型的步骤包括采用所述Van der Waals方程计算超
临界状态下甲烷气体的
密度和体积模量：
[0023]
[0024]
[0027] a——气体有关的参数，对于甲烷a＝0.225Pa(m3/mol)2＝879.9Pa(m3/kg)2；
[0028] b——气体有关的参数，对于甲烷b＝42.7cm3/mol＝2.675×10-3m3/kg；
[0029] ρg——气体密度；
[0030] Kg——气体体积模量；
[0031] R——气体常数，R＝8.31J/(mol°K)＝519.4J/(kg°K)。
[0032] 进一步地，在所述步骤3中，描述非饱和双重孔隙介质地震波传播的Biot-Rayleigh方程如下：
[0033]
[0034]
[0035]
[0036]
[0037]
[0038]
[0039] ui、Ui——分别为固相和液相的位移；
[0040] e、ξ——分别为固相、液相的体应变；
[0041] ζ——地震波激励过程中产生的一种局域性的流体
变形增量；
[0042] ρf1——背景相流体的密度；
[0043] R0——嵌入体(气泡)半径；
[0044] bi——(i＝1，2)分别表示两种流体饱和区域中各自的耗散系数；
[0045] κ10——表示主骨架(背景相)中的渗透率；
[0046] η1——表示背景相中的流体黏度；
[0047] φ——岩石中总的孔隙度，φ＝φ1+φ2＝v1φ10+v2φ20，v1与v2分别表示两种流体所占据的骨架占总骨架的体积比率，φ10与φ20分别表示两个区域内部的局部孔隙度；
[0048] A、N、Q1、R1、Q2、R2——双孔介质中的六个Biot弹性参数；
[0049]
[0050]
[0051]
[0052]
[0053]
[0054]
[0055] N＝μb
[0056] ks，kb，μs，μb——分别为岩石基质和骨架的体积模量、剪切模量；
[0057] ρ11、ρ12、ρ13、ρ22、ρ33——表示双孔介质中的五个密度参数，满足下式：
[0058] (1-φ1-φ2)ρs＝ρ11+ρ12+ρ13
[0059]
[0060]
[0061]
[0062]
[0063] 进一步地，在所述步骤3中，所述生成纵波、横波的相速度及逆品质因子，包括：以位移场ui和Ui表示在非饱和双重孔隙介质中传播的弹性平面波的通解形式，将位移场平面波解 j＝x，y，z.代入到所述Biot-Rayleigh方程，得到Christoffel方程，在A和 有非零解的情况下，根据Christoffel方程的行列式等于零得到纵波和横波的频散关系，最后由复
波数 可求得弹性波的相速度Vi和逆品质因子
[0064]
[0065]
[0066] 为实现上述目的，本发明还提供一种非饱和双重孔隙介质地震波频散衰减分析装置，该装置包括：岩石骨架弹性模量生成模
块，用于根据包括渗透率、孔隙率、泥质含量及矿物成分的地质数据，采用Voigt-Reuss-Hill平均模型生成骨架基质的等效弹性模量，并应用等效介质理论或者
接触理论计算干燥岩石的等效体积模量和剪切模量；储层流体模型建立模块，用于根据包括钻井数据、录井数据、测井信息及流体实验测量数据的测量数据，根据流体相态分析，应用Van der Waals方程或者流体试验测量结果建立储层流体模型；非饱和双重孔隙介质地震岩石物理建模模块，用于输入岩石骨架模型和流体模型，应用双重孔隙介质Biot-Rayleigh理论建立非饱和储层岩石物理模型；频散特性分析模块，用于分析非饱和岩石中的纵横波频散与衰减随频率、孔隙度及含气饱和度的变化规律。
[0067] 本发明
实施例的有益效果在于，本发明基于Rayleigh理论描述纵波激励下非饱和岩石中气泡的局域流体流动，从经典力学的哈密顿原理导出双重孔隙介质中的波传播方程，即Biot-Rayleigh方程，该方程格式简洁，所有相关参数物理可测，具有较好的物理可实现性，可用于分析非饱和岩石中的纵横波频散与衰减随频率、孔隙度、含气饱和度的变化规律。
附图说明
[0068] 为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。在附图中：
[0069] 图1为本发明实施例非饱和双重孔隙介质地震波频散衰减分析方法
流程图；
[0070] 图2为本发明实施例岩石干骨架建模流程图；
[0071] 图3为本发明实施例储层流体建模流程图；
[0072] 图4为本发明实施例非饱和双重孔隙介质地震波频散衰减计算流程图；
[0073] 图5为本发明实施例非饱和双重孔隙介质地震波频散衰减分析装置结构
框图；
[0074] 图6及图7分别为三种理论预测纵波速度、衰减的结果对比示意图；
[0075] 图8及图9分别为频率为0.1～10MHz，饱和度为0～100％时纵波速度和衰减变化示意图；
[0076] 图10为实验室观测数据和Biot-Rayleigh理论预测数据在地震频率，纵波速度随饱和度变化趋势图；
[0077] 图11为Biot-Rayleigh理论预测的含70％甲烷水
砂岩纵波速度示意图；
[0078] 图12为预测的横波速度示意图；
[0079] 图13为泥岩+70％甲烷水砂岩的AVO特征曲线示意图。
具体实施方式
[0080] 为使本发明实施例的目的、技术方案和优点更加清楚明白，下面结合附图对本发明实施例做进一步详细说明。在此，本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。
[0081] 如图1所示，本实施例提供一种非饱和双重孔隙介质地震波频散衰减分析方法，该方法包括：
[0082] 步骤S101：获取包括渗透率、孔隙率、泥质含量及矿物成分的地质数据，并根据所述地质数据生成储层岩石干骨架模型。
[0083] 图2为本发明岩石干骨架建模(建立储层岩石干骨架模型)流程图，其具体实现原理如下：
[0084] 已知岩石中各种矿物成分与其弹性模量，而矿物间的组合方式未知的情况下，Voigt界限和Reuss界限给出了岩石弹性模量的上下限。
[0085] Voigt界限又被称为等应变平均，其计算式如下：
[0086]
[0087] Reuss界限又被称为等
应力平均，其计算式如下：
[0088]
[0089] (1)、(2)两个公式中，N表示构成成分的个数，fi为第i个介质的体积分量，Mi为第i个介质的弹性模量，MV表示Voigt上限，MR表示Reuss下限。M可以代表任何一种模量：K，μ，E等。
[0090] Voigt-Reuss-Hill平均模型是Voigt上限和Reuss下限二者的算术平均，其表达式如下：
[0091]
[0092] (3)式中，MVRH即为Voigt-Reuss-Hill平均模量，可以表示任何一种弹性模量，也是其后计算岩石干骨架模量时所需的岩石骨架基质模量。
[0093] 然后，采用 (KT)包体模型计算干燥岩石的等效体积模量和剪切模量：
[0094]
[0095]
[0096] 其中， —— 理论求得的等效体积模量和剪切模量；
[0097] Km，μm——骨架基质的体积模量和剪切模量；
[0098] Ki，μi——第i个包体的体积模量和剪切模量；
[0099] xi——每种包体的体积比率；
[0100] N——包体的个数；
[0101] ζm＝μm(9Km+8μm)/6(Km+2μm)；mi mi
[0102] P ，Q ——骨架基质m中加入包体材料i后的弹性模量。
[0103] 步骤S102：获取包括钻井数据、录井数据、测井信息及流体实验测量数据的测量数据，根据所述测量数据及Van der Waals方程生成储层流体模型。
[0104] 图3为本发明实施例储层流体建模流程图，在步骤S102中，采用Van der Waals方程计算目的层温度、压力下的甲烷气体的密度和体积模量：
[0105]
[0106]
[0107] 其中，公式(6)、(7)P——地层压力；
[0108] T——地层温度；
[0109] a——气体有关的参数，对于甲烷a＝0.225Pa(m3/mol)2＝879.9Pa(m3/kg)2；
[0110] b——气体有关的参数，对于甲烷b＝42.7cm3/mol＝2.675×10-3m3/kg；
[0111] ρg——气体密度；
[0112] Kg——气体体积模量；
[0113] R——气体常数，R＝8.31J/(mol°K)＝519.4J/(kg°K)。
[0114] 步骤S103：根据所述储层岩石干骨架模型、储层流体模型及描述非饱和双重孔隙介质地震波传播的Biot-Rayleigh方程求解平面波解，生成纵波、横波的相速度及逆品质因子。
[0115] 图4为本发明实施例非饱和双重孔隙介质地震波频散衰减计算流程图，描述非饱和双重孔隙介质地震波传播的Biot-Rayleigh方程组如下：
[0116]
[0117]
[0118]
[0119]
[0120]
[0121]
[0122] (8)式中，ui、Ui——分别为固相和液相的位移；
[0123] e、ξ——分别为固相、液相的体应变；
[0124] ζ——地震波激励过程中产生的一种局域性的流体变形增量；
[0125] ρf1——背景相流体的密度；
[0126] R0——嵌入体(气泡)半径；
[0127] bi——(i＝1，2)分别表示两种流体饱和区域中各自的耗散系数；
[0128] κ10——表示主骨架(背景相)中的渗透率；
[0129] η1——表示背景相中的流体黏度；
[0130] φ——岩石中总的孔隙度，φ＝φ1+φ2＝v1φ10+v2φ20，v1与v2分别表示两种流体所占据的骨架占总骨架的体积比率，φ10与φ20分别表示两个区域内部的局部孔隙度；
[0131] A、N、Q1、R1、Q2、R2——双孔介质中的六个Biot弹性参数；
[0132] (9)
[0133]
[0134]
[0135]
[0136]
[0137]
[0138] N＝μb
[0139] ks，kb，μs，μb——分别为岩石基质和骨架的体积模量、剪切模量；
[0140] ρ11、ρ12、ρ13、ρ22、ρ23——表示双孔介质中的五个密度参数，满足下式：
[0141] (1-φ1-φ2)ρs＝ρ11+ρ12+ρ13 (14)
[0142]
[0143]
[0144]
[0145]
[0146]
[0147] 接下来求解Biot-Rayleigh方程组的平面波解，得到纵波和横波的相速度Vi和逆品质因子
[0148]
[0149]
[0150] 步骤S104：根据所述纵波、横波的相速度及逆品质因子生成速度、衰减、波阻抗、AVO响应特征等与频率之间的关系，根据上式可以分析得到非饱和岩石中的纵波频散与衰减随频率、孔隙度、含气饱和度的变化规律。
[0151] 如图5所示，本发明提供一种非饱和双重孔隙介质地震波频散衰减分析装置，该非饱和双重孔隙介质地震波频散衰减分析装置包括：岩石骨架弹性模量生成模块501，储层流体模型建立模块502，非饱和双重孔隙介质地震岩石物理建模模块503及频散特性分析模块504。
[0152] 岩石骨架弹性模量生成模块501用于根据包括渗透率、孔隙率、泥质含量及矿物成分的地质数据，采用Voigt-Reuss-Hill平均模型生成骨架基质的等效弹性模量，并应用等效介质理论或者接触理论计算干燥岩石的等效体积模量和剪切模量。
[0153] 根据目的区的文献调研、地质报告、测井资料、岩芯照片与电镜分析得到矿物组分、矿物体积比率、渗透率、孔隙率、泥质含量等，采用Voigt-Reuss-Hill平均模型计算骨架基质的等效弹性模量，应用等效介质理论或者接触理论计算干燥岩石的等效体积模量和剪切模量，如上述公式(1)至(5)所示。
[0154] 储层流体模型建立模块502用于根据包括钻井数据、录井数据、测井信息及流体实验测量数据的测量数据，根据流体相态分析，应用Van der Waals方程或者流体试验测量结果建立储层流体模型。
[0155] 根据钻井、录井、测井资料和文献调研等得到目的区的温度、压力信息，根据流体相态分析，应用Van der Waals方程或者流体试验测量结果建立储层流体模型，如上述公式(6)、(7)所示.
[0156] 非饱和双重孔隙介质地震岩石物理建模模块503用于输入岩石骨架模型和流体模型，应用双重孔隙介质Biot-Rayleigh理论建立非饱和储层岩石物理模型，如公式(8)至(21)所示。
[0157] 频散特性分析模块504用于分析非饱和岩石中的纵横波频散与衰减随频率、孔隙度及含气饱和度的变化规律。
[0158] 下面以具体的实例详细说明本发明：
[0159] 实例1
[0160] Biot-Rayleigh理论、Johnson理论和White理论预测纵波速度、衰减的结果对比：
[0161] 为检查拓展到非饱和岩石状况的Biot-Rayleigh理论的合理性，设计一个非饱和岩石物理模型，应用Biot-Rayleigh理论、Johnson理论和White三种理论进行了纵波速度和能量衰减的预测，对预测结果进行了对比分析。岩石参数为：基质体积模量35GPa，骨架体积模量7Gpa，骨架剪切模量9Gpa，
水体积模量2.25GPa，空气体积模量0.1MPa，水黏度3 3
0.001Pa*s，空气黏度0.00001Pa*s，基质平均密度2.65g/cm，水密度0.99g/cm，空气密度
3
0.1g/cm，孔隙度为0.15，渗透率0.1D，平均气泡尺寸0.25m，空气饱和度为0.95。图6及图7分别为三种理论预测纵波速度、衰减的结果对比示意图。
[0162] 实例2
[0163] 全频带非饱和岩石纵波的速度与衰减的定量预测：
[0164] 设计一个非饱和岩石物理模型，基本的岩石物理参数为：
石英体积模量39GPa，砂9 5
岩固结系数15，水体积模量2.51×10Pa，气体积模量1.44×10Pa，水黏度0.001Pa*s，气黏
3 3 3
度0.000022Pa*s，砂岩颗粒平均密度2.69g/cm，水密度1.04g/cm，气密度0.01g/cm，平均气泡尺寸lmm，孔隙度为0.11，渗透率0.5mD。图8及图9分别为频率为0.1～10MHz，饱和度为0～100％时纵波速度和衰减变化示意图。
[0165] 实例3
[0166] 北海高孔砂岩宽频带非饱和岩石波速观测数据对比：
[0167] 为验证本文理论成果的适用性，采用了Adam与Batzle于2006年发表的宽频带非饱和岩石波速观测数据与本发明Biot-Rayleigh理论预测的结果进行对比分析。高孔砂岩来自北海，主要成分为石英，孔隙度为0.35，渗透率为8.7D，气泡平均尺寸为5mm，图10为实验室观测数据和Biot-Rayleigh理论预测数据在地震频率(＜100Hz)、测井频率(50KHz)、超声频率下(500kHz)，纵波速度随饱和度变化趋势图。
[0168] 实例4
[0169] 含非饱和甲烷砂岩纵波速度、衰减及AVO特征：
[0170] 以南海渐新统陵水组
长石岩屑质石英砂岩为例设计一个非饱和甲烷岩石物理3
模型，基本的岩石物理参数为：石英体积模量37GPa，剪切模量44Gpa，密度2.65g/cm ；长
3
石体积模量37.5GPa，剪切模量15Gpa，密度2.62g/cm ；粘土体积模量21GPa，剪切模量
3 9 3
7Gpa，密度2.6g/cm ；水体积模量2.5×10Pa，水密度1.04g/cm，埋深3.5Km下甲烷体积模
3
量0.081GPa，密度0.17g/cm ；孔隙度为0.15，甲烷饱和度为70％，盖层为泥岩，纵波速度
3400m/s，横波速度1600m/s。图11为Biot-Rayleigh理论预测的含70％甲烷水砂岩纵波速度示意图，图12为预测的横波速度示意图，图13为泥岩+70％甲烷水砂岩的AVO特征曲线示意图。
[0171] 本发明实施例的有益效果在于，本发明基于Rayleigh理论描述纵波激励下非饱和岩石中气泡的局域流体流动，从经典力学的哈密顿原理导出双重孔隙介质中的波传播方程，即Biot-Rayleigh方程，该方程格式简洁，所有相关参数物理可测，具有较好的物理可实现性，可用于分析非饱和岩石中的纵横波频散与衰减随频率、孔隙度、含气饱和度的变化规律。
[0172] 以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何
修改、等同替换、改进等，均应包含在本发明的保护范围之内。
