一种基于Norris-KG模型的地震纵波模拟分析方法与流程

drake，1972)，对计算机内存要求很高，计算量大。伪谱法在地震领域的应用则大大加快了计算速度，并节省了计算机存储空间(gazdag，1981；kosloff and baysal，1982)[39-40]，但该方法对边界的处理较难。有限差分法是目前声波或弹性波场数值模拟中应用最广的方法(alterman and karal，1968；boore，1972；alford et al.，1974；carcione，1996；牟永光，1996；孙卫涛，2009)，其实现简单，运行速度快，但难以克服空间离散带来的数值频散问题。高阶差分方法(virieux，1986)和一阶弹性波交错网格方法(dablain，1986)的提出与结合，形成了一阶交错网格高阶有限差分方法(董良国等，2000)，该方法有力地减弱了常规有限差分法的数值频散问题，因而在弹性介质(裴正林，2004)、粘弹性介质(李军等，2009；韩令贺等，2010)、非均匀介质(裴正林,牟永光，2003)、各向异性介质(裴正林，2006)、孔隙介质(王秀明等，2003)、裂缝-孔隙介质(杜启振等，2009；孔丽云等，2013)的正演模拟中取得了相对成功的应用。
[0005]
纵横波波场分离方面，目前常用的方法是亥姆霍兹分解(yan和sava，2008)。该方法所依据的基本原理是由于纵横波偏振方向的不同(纵波偏振方向平行于传播方向，横波偏振方向垂直于传播方向)，可以用波场的散度和旋度来分别对纵波和横波进行求取。这在各向同性的情况下是正确的，但是在各向异性介质中，地震波传播的方向发生了偏离，其振动方向并不再是纵波与其相同、横波与其垂直的关系，所以使用亥姆霍兹分解并不能很好地将各向异性介质中的耦合纵横波分离(yan和sava，2009)。另外，由于在分离过程中采用散度场代表纵波，旋度场代表横波，这就意味着，纵波仅代表整体位移的输入输出，而横波则代表旋转变换的角速度，其物理意义与原波场并不相同，且其振幅和相位信息都发生了改变。qp波近似方程方面，alkhalifah从vti介质的精确qp-qsv波频散关系出发，假设垂向剪切波速度v
s0
＝0，推导了vti介质四阶拟声波方程。为降低四阶偏微分方程计算的复杂性(alkhalifah，2008)，zhou等(2009)和du等(2008)从声学近似的频散关系出发，通过引入不同的辅助波场函数，分别推导了两种vti介质二阶耦合拟声波方程。为结合交错网格技术实现精确的有限差分数值模拟，hestholm推导了一阶偏微分形式的vti介质拟声波方程(hestholm，2009)。duveneck等(2008)从声学近似的胡克定律出发，也推导了二阶形式和一阶应力－速度形式的vti介质拟声波方程。与alkhalifah方程一样，几种拟声波方程均存在低速度、低振幅的qsv人为干扰波问题(grechka等，2004)，而这种干扰波的存在会影响正演模拟和偏移成像的最终结果。同时，由于无法考虑到转换波的信息，利用qp波近似方程得到的各向异性介质纵波数据并不完整。为得到振幅和相位信息相对准确且波场信息相对完整的各向异性介质纵波数据，zhang等提出“方程解耦”的分离方法(zhang等，2007)。在各向同性介质中，zhang等(2007)得到精确的纵横波分离的一阶速度应力方程；zhou等(2016)将这种方法进一步拓展到三维各向同性模型的情况。孔丽云等(2022)通过将弹性刚度系数分量中的耦合项c
13
近似为纵波模量和横波模量的线性函数，推导得到了各向异性介质的纵横波分离的一阶速度-应力方程，从而能够实现各向异性介质的纵横波波场分离。
[0006]
综上可知,目前的研究大多将等效介质建模、正演数值模拟、纵横波场分离这三方面的技术分开研究，无法直观体现裂缝性储层纵波地震响应特征与储层物性参数之间的联系。


技术实现要素：

[0007]
针对上述

背景技术：
中存在的问题，本发明提出了一种基于norris-kg模型的地震纵波模拟分析方法，其构思合理，以储层物性参数为输入，通过三方面理论的融合、计算，能够直接给出储层的纵波地震响应，并通过分析研究纵波地震响应特征随储层物性参数的变化规律，有助于揭示裂缝性储层地震响应特征的产生机理，为非常规储层的勘探开发提供可靠的理论依据。
[0008]
为解决上述技术问题，本发明提供的一种基于norris-kg模型的地震纵波模拟分析方法，基于norris-kg裂缝-孔隙等效介质模型，通过qp波相速度公式及纵、横波分离的一阶速度应力方程对介质的波场控制方程进行推导，并利用基于pml边界条件的交错网格高阶有限差分算法对地震波场进行正演数值模拟，最终，在给定参数下，对目的层的qp波相速度方位各向异性及其顶、底界面的反射波振幅方位各向异性分析，对目的层的各向异性及其影响参数进行优选分析,得到裂缝-孔隙介质地震纵波相速度和波场响应数据，最后进行纵波数值模拟分析。
[0009]
所述基于norris-kg模型的地震纵波模拟分析方法，其中，所述波场控制方程的推导过程为：
[0010]
1.1)构建qp波相速度方程
[0011]
假设裂缝倾角为θ0，地震波的传播方向为n＝(n
x
,ny,nz)
t
＝(sinθcosφ,sinθsinφ,cosθ)
t
，其中，θ为传播方向与观测系统z轴的夹角，φ为传播方向与x轴的夹角，任意入射角的纵波相速度表达式为：
[0012][0013]
其中，
[0014][0015][0016][0017][0018]
上式(2)建立了裂缝-多孔隙介质的纵波相速度vqp与储层弹性参数c之间的数学联系，结合norris-kg模型所给出的介质弹性模量表达式(1)，通过数值计算，得到介质的纵波相速度随储层物性参数的变化规律；
[0019]
1.2)构建纵横波分离的一阶速度-应力方程
[0020]
基于“方程分离”的思想，将弹性刚度系数分量中的耦合项c
13
近似为纵波模量和横波模量的线性函数，推导得到了各向异性介质的纵、横波分离的一阶速度-应力方程如下：
[0021][0022]
其中，上述方程(3)中下角标p和s分别代表纵波和横波，g为偏微分算子矩阵,c
p
为
与纵波相关的弹性系数矩阵,cs为与横波相关的弹性系数矩阵；且c
p
与cs表达式分别如下：
[0023][0024][0025]
上述的c
p
与cs表达式中,ε和γ为thomsen各向异性参数,v
p0
为qp波沿vti介质对称轴方向的相速度,v
s0
为qsv波沿vti介质对称轴方向的相速度；和b＝-2分别为纵波和横波模量的系数；
[0026]
上述方程(3)建立了裂缝-多孔隙介质的地震纵波波场v与储层物性参数c之间的数学联系，通过对波动方程的数值求解，可得到介质的纵波地震响应v
p
，进而对储层纵波地震响应特征随物性参数的变化规律进行分析。
[0027]
所述基于norris-kg模型的地震纵波模拟分析方法，其中，所述纵波数值模拟分析的具体过程为：
[0028]
2.1)qp波相速度分析
[0029]
2.1.1)基准模型结果分析
[0030]
基于任意入射角的纵波相速度公式(2)，根据给定井段的岩石物理测量数据和测井数据，对目的层qp波的相速度及储层物性参数的影响进行分析；
[0031]
2.1.2)储层参数的影响分析
[0032]
在给定参数的基础上，分别改变裂缝倾角、裂缝线密度、基质孔隙度、含气饱和度等参数，以分析不同储层物性参数对qp波相速度的影响；
[0033]
2.2)波场传播特征分析
[0034]
2.2.1)基准模型结果分析
[0035]
设计三层理论模型对目的层顶、底界面的地震波场传播特征进行分析，然后基于所述三层模型的参数，根据一阶速度-应力方程(3)，采用基于pml边界条件的交错网格高阶有限差分方法，对所述三层模型进行正演数值模拟；
[0036]
2.2.2)储层参数的影响分析
[0037]
在给定参数的基础上，保持其他参数不变，分别改变裂缝倾角、裂缝线密度、基质孔隙度和含气饱和度，以分析不同的参数对地震波场响应特征的影响。
[0038]
所述基于norris-kg模型的地震纵波模拟分析方法，其中，所述正演数值模拟是采
用交错网格高阶有限差分法对介质的一阶速度-应力方程进行数值求解；具体包括以下步骤：
[0039]
3.1)方程的离散化处理
[0040]
采用交错网格高阶有限差分算法，其时间、空间高阶差分近似公式分别为：
[0041][0042][0043]
其中，上式(4)-(5)中δt为时间步长，2m为时间差分精度，2l为空间差分精度；δx为差分网格间距，为不同阶的交错网格空间差分系数，可以通过如下方程组求取:
[0044][0045]
将该方程组带入到上述方程(3)中得到一阶速度应力方程的差分近似格式；
[0046]
3.2)pml边界条件处理
[0047]
采用pml边界条件处理方法，对上述方程(3)中的总体波场v进行变量分离，每一个波场变量可以分为纵波分量vp和横波分量vs两部分，t是总的应力场，上角标x、y、z代表三个坐标轴方向：
[0048][0049]
然后对于所给定的三维空间模型设置边界衰减系数。
[0050]
所述基于norris-kg模型的地震纵波模拟分析方法，其中，所述norris-kg裂缝-孔隙等效介质模型的具体构建过程为：基于norris周期成层模型和kg裂缝-孔隙等效介质模型建立norris-kg等效介质模型；基于norris-kg等效介质模型，假设裂缝为非常薄且松软的无限大裂缝平面，将norris-kg等效介质模型看成是周期成层两套裂缝-孔隙地层，根据norris理论及heuristic假设，norris-kg等效介质模型的弹性属性表达式为：
[0051][0052]
上式(1)中，为垂直裂缝面方向纵波模量，为模型高频极限弹性系数，为垂直裂缝面方向的高频极限弹性系数，为模型低频极限弹性系数；式(1)给出的是水平裂缝性储层的弹性系数矩阵，可计算任意入射角情况的频变纵波速度随流体饱和度的变化情况。
[0053]
采用上述技术方案，本发明具有如下有益效果：
[0054]
本发明基于norris-kg模型的地震纵波模拟分析方法是在已知储层物性参数(颗粒类型、孔隙形态、流体属性等)的情况下，给出地震纵波响应特征(速度、振幅等)与储层物性参数之间的对应关系，让我们有能力从地震纵波响应特征中直接对储层的物性参数进行识别。
[0055]
本发明将norris-kg等效介质模型、纵横波分离算法、数值模拟计算三者相融合，实现了基于norris-kg模型的地震纵波模拟分析。在给定参数下，通过对某研究区目的层的qp波相速度方位各向异性及其顶、底界面的反射波振幅方位各向异性分析，对目的层的各向异性及其影响参数进行优选。结果表明，在给定的岩石物理参数和测井参数条件下，该研究区目的层的观测各向异性较为明显，且主要受裂缝倾角和裂缝线密度的影响。
附图说明
[0056]
为了更清楚地说明本发明具体实施方式或现有技术中的技术方案下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0057]
图1为本发明基于norris-kg模型的地震纵波模拟分析方法中涉及的norris-kg模型示意图；
[0058]
图2为本发明基于norris-kg模型的地震纵波模拟分析方法中涉及的三维交错网格示意图；
[0059]
图3为本发明基于norris-kg模型的地震纵波模拟分析方法中涉及的三维完全匹配层吸收边界示意图；
[0060]
图4为本发明基于norris-kg模型的地震纵波模拟分析方法中涉及的qp波相速度随入射角和观测方位角的变化(红、紫、蓝线代表入射角分别为17
°
，32
°
和47
°
的情况)示意图；
[0061]
图5为本发明基于norris-kg模型的地震纵波模拟分析方法中涉及的不同裂缝倾角、地震波入射角下的qp波相速度随方位角的变化规律((a)-(i)分别对应裂缝倾角10
°‑
90
°
，间隔为10
°
；红、紫、蓝线代表入射角分别为17
°
，32
°
和47
°
的情况)示意图；
[0062]
图6为本发明基于norris-kg模型的地震纵波模拟分析方法中涉及的裂缝线密度分别为(a)1条/m、(b)15条/m和(c)30条/m时的qp波相速度随方位角的变化规律(红、紫、蓝线代表入射角分别为17
°
，32
°
和47
°
的情况)示意图；
[0063]
图7为本发明基于norris-kg模型的地震纵波模拟分析方法中涉及的基质孔隙度分别为(a)1.12％、(b)3.40％和(c)5.81％时的qp波相速度随方位角的变化规律(红、紫、蓝线代表入射角分别为17
°
，32
°
和47
°
)示意图；
[0064]
图8为本发明基于norris-kg模型的地震纵波模拟分析方法中涉及的含气饱和度分别为(a)0％、(b)25％、(c)50％、(d)75％和(e)100％时的qp波相速度(红、紫、蓝线代表入射角分别为17
°
，32
°
和47
°
的情况)示意图；
[0065]
图9为本发明基于norris-kg模型的地震纵波模拟分析方法中涉及的层顶界面(a)、底界面(b)的叠前全方位角道集模拟数据示意图；
[0066]
图10为本发明基于norris-kg模型的地震纵波模拟分析方法中涉及的目的层顶界面(a)、底界面(b)的振幅方位各向异性示意图；
[0067]
图11为本发明基于norris-kg模型的地震纵波模拟分析方法中涉及的裂缝倾角分别为1条/m(左)、15条/m(中)和30条/m(右)时，目的层顶界面蜗牛道集((a)-(c))、目的层底界面蜗牛道集((d)-(f))、目的层顶界面反射波振幅随方位角的变化((g)-(i))、目的层底界面反射波振幅随方位角的变化((j)-(l))示意图；
[0068]
图12为本发明基于norris-kg模型的地震纵波模拟分析方法中涉及的裂缝线密度分别为(a)1条/m、(b)15条/m和(c)30条/m时，目的层顶界面蜗牛道集((a)-(c))、目的层底界面蜗牛道集((d)-(f))、目的层顶界面反射波振幅随方位角的变化((g)-(i))、目的层底界面反射波振幅随方位角的变化((j)-(l))示意图；
[0069]
图13为本发明基于norris-kg模型的地震纵波模拟分析方法中涉及的基质孔隙度分别为(a)1.12％、(b)3.40％和(c)5.81％时，目的层顶界面蜗牛道集((a)-(c))、目的层底界面蜗牛道集((d)-(f))、目的层顶界面反射波振幅随方位角的变化((g)-(i))、目的层底界面反射波振幅随方位角的变化((j)-(l))示意图；
[0070]
图14为本发明基于norris-kg模型的地震纵波模拟分析方法中涉及的含气饱和度分别为(a)0％、(b)50％和(c)100％时，目的层顶界面蜗牛道集((a)-(c))、目的层底界面蜗牛道集((d)-(f))、目的层顶界面反射波振幅随方位角的变化((g)-(i))、目的层底界面反射波振幅随方位角的变化((j)-(l))示意图。
具体实施方式
[0071]
下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0072]
下面结合具体的实施方式对本发明做进一步的解释说明。
[0073]
本实施例提供的基于norris-kg模型的地震纵波模拟分析方法，基于norris-kg裂缝-孔隙等效介质模型，通过qp波相速度公式及纵、横波分离的一阶速度应力方程对介质的波场控制方程进行推导，并利用基于pml边界条件的交错网格高阶有限差分算法对地震波场进行正演数值模拟，最终，在给定参数下，对目的层的qp波相速度方位各向异性及其顶、底界面的反射波振幅方位各向异性分析，对目的层的各向异性及其影响参数进行优选分析,得到裂缝-孔隙介质地震纵波相速度和波场响应数据，最后进行纵波数值模拟分析。
[0074]
本实施例提供的基于norris-kg模型的地震纵波模拟分析方法，主要包括以下步骤：
[0075]
s100)norris-kg等效介质建模
[0076]
基于norris周期成层模型和kg裂缝-孔隙等效介质模型，建立图1所示的部分饱和裂缝-孔隙等效介质模型(以下称norris-kg模型)。图中，模型周期为h，两套地层均由干燥的kg模型所组成，蓝色代表的是水饱和kg模型(流体体积模量和粘度系数分别为k
water
、η
water
)，灰色代表的是气饱和kg模型(流体体积模量和粘度系数分别为k
gas
、η
gas
)。这样，含水饱和度sw和含气饱和度sg就分别代表水饱和和气饱和地层的空间占比。
[0077]
基于图1所示的模型，假设裂缝为非常薄且松软的无限大裂缝平面(符合kg模型假
设条件)，将norris-kg模型看成是周期成层两套裂缝-孔隙地层(单一流体饱和)，根据norris理论及heuristic假设(水平层状介质中，流体的能量传播方向始终垂直水平地层)，模型的弹性属性表达式为：
[0078][0079]
其中，为垂直裂缝面方向纵波模量，分别为模型高、低频极限弹性系数。
[0080]
式(1)给出的是水平裂缝性储层(即vti介质，vertical transverse isotropy)的弹性系数矩阵，能够计算任意入射角情况的频变纵波速度随流体饱和度的变化情况。
[0081]
s200)波场控制方程推导
[0082]
s210)qp波相速度方程
[0083]
假设裂缝倾角为θ0，地震波的传播方向为n＝(n
x
,ny,nz)
t
＝(sinθcosφ,sinθsinφ,cosθ)
t
，其中，θ为传播方向与观测系统z轴的夹角，φ为传播方向与x轴的夹角，任意入射角的纵波相速度表达式：
[0084][0085]
其中，
[0086][0087][0088][0089][0090]
式(2)建立了裂缝-多孔隙介质的纵波相速度vqp与储层弹性参数c之间的数学联系，结合norris-kg模型所给出的介质弹性模量表达式(1)，通过数值计算，可以得到介质的纵波相速度随储层物性参数的变化规律。
[0091]
s220)纵横波分离的一阶速度-应力方程
[0092]
基于“方程分离”的思想，将弹性刚度系数分量中的耦合项c
13
近似为纵波模量和横波模量的线性函数，推导得到了各向异性介质的纵、横波分离的一阶速度-应力方程：
[0093][0094]
其中，下角标p和s分别代表纵波和横波，g为偏微分算子矩阵；c
p
和cs分别为与纵波和横波相关的弹性系数矩阵，c
p
和cs的表达式如下：
[0095][0096][0097]
ε和γ为thomsen各向异性参数；v
p0
、v
s0
分别为qp波和qsv波沿vti介质对称轴方向的相速度；和b＝-2分别为纵波和横波模量的系数。
[0098]
表达式(3)建立了裂缝-多孔隙介质的地震纵波波场v与储层物性参数c之间的数学联系，通过对波动方程的数值求解，可以得到介质的纵波地震响应v
p
，进而对储层纵波地震响应特征随物性参数的变化规律进行分析。
[0099]
s300)正演数值模拟
[0100]
交错网格高阶有限差分法是通过将速度(应力)对时间的奇数阶高阶导数转化为应力(速度)对空间的导数，运用时间和空间差分精度均可达到任意阶的高阶差分法，采用交错网格技术，对介质的一阶速度-应力方程进行数值求解的过程。
[0101]
s310)方程的离散化处理
[0102]
采用交错网格高阶有限差分算法，其时间、空间高阶差分近似公式分别为：
[0103][0104][0105]
其中，上式(4)-(5)中δt为时间步长，2m为时间差分精度，2l为空间差分精度；δx为差分网格间距，为不同阶的交错网格空间差分系数，可以通过如下方程组求取:
[0106][0107]
将上述表达式带入到方程(3)中，就得到一阶速度应力方程的差分近似格式。
[0108]
进一步，时间上，将速度场放在主时间网格上，应力场放在次时间网格上；空间上，
速度场和应力场的位置如图2所示。这样，就可以得到差分近似格式的离散化形式。
[0109]
s320)pml边界条件处理
[0110]
采用pml边界条件处理方法，对上述方程(3)中的总体波场v进行变量分离，每一个波场变量可以分为纵波分量vp和横波分量vs两部分，t是总的应力场，上角标x、y、z代表三个坐标轴方向：
[0111][0112]
对于三维模型，按照图3设置边界衰减系数。边界面区：a，d
x
＝0，dy＞0，dz＝0；b，d
x
＝0，dy＝0，dz＞0；c，d
x
＞0，dy＝0，dz＝0。在边界棱区：d，d
x
＝0，dy＞0，dz＞0；e，d
x
＞0，dy＝0，dz＞0；f，d
x
＞0，dy＞0，dz＝0。在边界角点区：g，d
x
＞0，dy＞0，dz＞0。
[0113]
最终，通过步骤s100-s300，得到裂缝-孔隙介质地震纵波相速度和波场响应数据。
[0114]
s400)纵波数值模拟分析
[0115]
s410)qp波相速度分析
[0116]
s4101)基准模型结果分析
[0117]
基于任意入射角的纵波相速度公式(2)，根据给定井段的岩石物理测量数据和测井数据，对目的层qp波的相速度及储层物性参数的影响进行分析。
[0118]
岩石固体颗粒为白云石，体积模量、剪切模量和密度分别为94.9gpa，45gpa和2870kg/m3；储层内所含流体为气和水，其中，气的体积模量和密度分别为0.213gpa和258kg/m3，水的体积模量和密度分别为2.5gpa和1020kg/m3；储层中的基质孔隙度为5.81％；储层中裂缝的线密度为30条/m，裂缝平均宽度为0.0024m，裂缝倾角为90
°
。
[0119]
将以上参数带入到qp波相速度表达式(2)中，初始含气饱和度为100％，可以得到在给定参数下的qp波相速度随入射角和观测方位角的变化如图4所示。图中展示了入射角分别为17
°
(红线)，32
°
(紫线)和47
°
(蓝线)时qp波相速度随观测方位角从0
°‑
360
°
的变化规律。从图中可以看出，当裂缝倾角为90
°
时，随着观测方位角从0
°
到360
°
变化，qp波相速度呈现明显的方位各向异性特征，且入射角越大，方位各向异性越明显。
[0120]
s4102)储层参数的影响
[0121]
在给定参数的基础上，分别改变裂缝倾角、裂缝线密度、基质孔隙度、含气饱和度等参数，以分析不同储层物性参数对qp波相速度的影响。
[0122]
s41021)裂缝倾角的影响
[0123]
对不同裂缝倾角、地震波入射角以及观测方位角下的qp波相速度进行计算，结果如图5所示。其中，图5(a)-(i)分别对应裂缝倾角10
°‑
90
°
，间隔为10
°
；每幅图都展示了入射角分别为17
°
(红线)，32
°
(紫线)和47
°
(蓝线)时qp波相速度随观测方位角从0
°‑
360
°
的变化规律。
[0124]
从图5中可以看出，当入射角θ＝17
°
时(红线)，qp波相速度在所有裂缝倾角的情况下，几乎都关于观测方位角φ＝180
°
对称，这表明当入射角较小时，各向异性对qp波相速度的影响很小。当入射角θ＝47
°
时(紫线)，qp波相速度在大部分情况下关于φ＝180
°
是不对称的，且呈现以下规律：(1)当φ《180
°
时，裂缝倾角较小(《40
°
)时的qp波相速度随方位角出
现先增大后减小的趋势(图5(a)-(d))，且在φ＝90
°
时出现最大值，表明此时地震波传播方向与各向同性面(裂缝面)的平行方向最为接近；而当裂缝倾角较大时，地震波随方位角的增大逐渐透过裂缝面而后又返回，因此会出现qp波相速度先增大后减小然后又增大的趋势(图5(e)-(h)))；(2)当φ》180
°
时，qp波相速度随方位角的增大出现先增大后减小的趋势(图5(a)-(i))，且在φ＝270
°
时出现最小值，表明此时地震波传播方向与各向同性面(裂缝面)的垂直方向最为接近；(3)只有当裂缝倾角为90
°
时，qp波相速度是关于φ＝180
°
完全对称的(图5(i))。因此，对于各向异性较强的海相碳酸盐岩储层，其实际地震资料的采集和处理应尽可能考虑大角度(大于30
°
)入射时的全方位(0
°‑
360
°
)地震信息。
[0125]
s41022)裂缝线密度的影响
[0126]
图6给出了裂缝线密度分别为1条/m(图6(a))、15条/m(图6(b))和30条/m(图6(c))时的qp波相速度计算结果。从图中可以看到，随着裂缝线密度的增大，qp波相速度整体变小，且其方位各向异性逐渐增大。这是因为，随着裂缝线密度增大，等效介质的体积模量变小，会导致速度减小；而裂缝的增多会导致介质的等效各向异性增强。
[0127]
s41023)基质孔隙度的影响
[0128]
图7给出了基质孔隙度分别为1.12％(图7(a))、3.40％(图7(b))和5.81％(图7(c))时的qp波相速度计算结果。从图中可以看到，随着基质孔隙度的增大，qp波相速度整体变小，但其方位各向异性变化不大。这是由于，基质孔隙度增大，等效介质的弹性模量减小，会导致速度减小；但是各向异性不受基质孔隙度的影响。
[0129]
s41024)含气饱和度的影响
[0130]
图8给出了含气饱和度分别为0％(图8(a))、25％(图8(b))、50％(图8(c))、75％(图8(d))和100％(图8(e))时的qp波相速度计算结果。从图中可以看到，随着含气饱和度的增大，qp波相速度整体变小，但其方位各向异性变化不大。这是由于，随着当储层含气时，介质的等效体积模量减小，会导致速度减小；当储层含气后，混合流体的体积模量近似等于气体的体积模量，而等效介质的密度变小，从而导致整体的速度随着含气饱和度的增大而增大。但由于基质孔隙度较小(5.81％)，含气饱和度对速度的影响并不明显。
[0131]
s420)波场传播特征分析
[0132]
s4201)基准模型结果分析
[0133]
为了对目的层顶、底界面的地震波场传播特征进行分析，特设计三层理论模型。中间层是目的层，其参数与3.1中所给定的初始参数相同；模型顶层和底层均为各向同性弹性介质，通过测井资料给定其纵、横波速度和密度分别为：顶层，4729m/s，2606m/s，2600kg/m3；底层，5660m/s，3364m/s，2730kg/m3。
[0134]
基于上述模型参数，根据一阶速度-应力方程(3)，采用基于pml边界条件的交错网格高阶有限差分方法，对上述三层模型进行正演数值模拟。震源采用爆炸震源，主频25hz，时间间隔为0.4ms；模型大小为704m
×
704m
×
528m，网格间距为8m；交错网格差分法精度为o(δt2，δx
16
，δy
16
，δz
16
)。计算入射角分别为：0
°
、9
°
、17
°
、25
°
、32
°
、37
°
、43
°
、47
°
、51
°
、54
°
、57
°
、60
°
时的模型叠前全方位角道集数据，对于每一个入射角，方位角从0
°
到360
°
变化。图9给出了目的层顶界面(图9(a))、底界面(图9(b))的叠前全方位角道集数据，图10给出入射角分别为17
°
(红色)、32
°
(紫色)、47
°
(蓝色)时中间层顶界面(图10(a))、底界面(图10(b))的反射波振幅信息随观测方位角的变化曲线。
[0135]
从图9和10中可以看出，当入射角较小时(小于30
°
)，目的层顶、底界面反射波振幅的各向异性并不明显(图9，图10中红线)，随着入射角的增大，目的层顶界面反射波振幅呈现明显各向异性的特点(图9，图10紫线)，而底界面反射波的走时则呈现较强的各向异性特征(图9(b))。这表明，裂缝-多孔隙介质模型中纵波速度的各向异性会影响裂缝性储层的顶、底界面反射波走时和振幅的方位各向异性，可以利用qp波地震资料的这一特点对储层地震资料的各向异性时间、振幅属性进行分析，为储层的裂缝预测提供相应的理论依据。
[0136]
s4202)储层参数的影响
[0137]
在给定参数的基础上，保持其他参数不变，分别改变裂缝倾角、裂缝线密度、基质孔隙度和含气饱和度，以分析不同的参数对地震波场响应特征的影响。图11-14中，每组图的12幅图中，从上到下分别为目的层顶界面蜗牛道集、目的层底界面蜗牛道集、目的层顶界面反射波振幅随方位角的变化、目的层底界面反射波振幅随方位角的变化。从左向右，图11展示的是裂缝倾角分别为30
°
、60
°
、90
°
时的数值模拟结果，图12展示的是裂缝线密度分别为1条/m、15条/m、30条/m的数值模拟结果，图13展示的是基质孔隙度分别为1.12％、3.40％和5.81％的数值模拟结果，图14展示的是含气饱和度分别为0％、50％、100％时的数值模拟结果。
[0138]
s42021)裂缝倾角的影响
[0139]
从图11中可以看出，当裂缝倾角增大时，目的层顶、底界面的反射振幅的方位各向异性呈现增大的趋势(从左到右，反射波振幅随方位角变化的椭圆程度增强)，且不再关于180
°
方位角对称，这是由于qp波相速度的变化引起的。
[0140]
s42022)裂缝线密度的影响
[0141]
从图12中可以看出，当裂缝线密度变化时，目的层顶、底界面反射波振幅的方位各向异性随裂缝线密度的减小而减小(从右至左，反射波振幅随方位角变化的椭圆程度降低)。
[0142]
s42023)基质孔隙度的影响
[0143]
从图13中可以看出，当基质孔隙度变化时，目的层顶、底界面反射波振幅的方位各向异性基本不变，但是反射波振幅强度增强。
[0144]
s42024)含气饱和度的影响
[0145]
从图14中可以看出，当含气饱和度增大时，目的层顶、底界面反射波振幅的方位各向异性基本不变，但是反射波振幅强度增强。
[0146]
本发明将norris-kg等效介质模型、纵横波分离算法、数值模拟计算三者相融合，实现了基于norris-kg模型的地震纵波模拟分析。在给定参数下，通过对某研究区目的层的qp波相速度方位各向异性及其顶、底界面的反射波振幅方位各向异性分析，对目的层的各向异性及其影响参数进行优选。结果表明，在给定的岩石物理参数和测井参数条件下，该研究区目的层的观测各向异性较为明显，且主要受裂缝倾角和裂缝线密度的影响。
[0147]
最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

技术特征：
1.一种基于norris-kg模型的地震纵波模拟分析方法，其特征在于，基于norris-kg裂缝-孔隙等效介质模型，通过qp波相速度公式及纵、横波分离的一阶速度应力方程对介质的波场控制方程进行推导，并利用基于pml边界条件的交错网格高阶有限差分算法对地震波场进行正演数值模拟，最终，在给定参数下，对目的层的qp波相速度方位各向异性及其顶、底界面的反射波振幅方位各向异性分析，对目的层的各向异性及其影响参数进行优选分析,得到裂缝-孔隙介质地震纵波相速度和波场响应数据，最后进行纵波数值模拟分析。2.如权利要求1所述的基于norris-kg模型的地震纵波模拟分析方法，其特征在于，所述波场控制方程的推导过程为：1.1)构建qp波相速度方程假设裂缝倾角为θ0，地震波的传播方向为n＝(n
x
,n
y
,n
z
)
t
＝(sinθcosφ,sinθsinφ,cosθ)
t
，其中，θ为传播方向与观测系统z轴的夹角，φ为传播方向与x轴的夹角，任意入射角的纵波相速度表达式为：其中，其中，其中，其中，上式(2)建立了裂缝-多孔隙介质的纵波相速度vqp与储层弹性参数c之间的数学联系，结合norris-kg模型所给出的介质弹性模量表达式(1)，通过数值计算，得到介质的纵波相速度随储层物性参数的变化规律；1.2)构建纵横波分离的一阶速度-应力方程基于“方程分离”的思想，将弹性刚度系数分量中的耦合项c
13
近似为纵波模量和横波模量的线性函数，推导得到了各向异性介质的纵、横波分离的一阶速度-应力方程如下：其中，上述方程(3)中下角标p和s分别代表纵波和横波，g为偏微分算子矩阵,c
p
为与纵波相关的弹性系数矩阵,c
s
为与横波相关的弹性系数矩阵；且c
p
与c
s
表达式分别如下：
上述的c
p
与c
s
表达式中,ε和γ为thomsen各向异性参数,v
p0
为qp波沿vti介质对称轴方向的相速度,v
s0
为qsv波沿vti介质对称轴方向的相速度；和b＝-2分别为纵波和横波模量的系数；上述方程(3)建立了裂缝-多孔隙介质的地震纵波波场v与储层物性参数c之间的数学联系，通过对波动方程的数值求解，可得到介质的纵波地震响应v
p
，进而对储层纵波地震响应特征随物性参数的变化规律进行分析。3.如权利要求2所述的基于norris-kg模型的地震纵波模拟分析方法，其特征在于，所述纵波数值模拟分析的具体过程为：2.1)qp波相速度分析2.1.1)基准模型结果分析基于任意入射角的纵波相速度公式(2)，根据给定井段的岩石物理测量数据和测井数据，对目的层qp波的相速度及储层物性参数的影响进行分析；2.1.2)储层参数的影响分析在给定参数的基础上，分别改变裂缝倾角、裂缝线密度、基质孔隙度、含气饱和度等参数，以分析不同储层物性参数对qp波相速度的影响；2.2)波场传播特征分析2.2.1)基准模型结果分析设计三层理论模型对目的层顶、底界面的地震波场传播特征进行分析，然后基于所述三层模型的参数，根据一阶速度-应力方程(3)，采用基于pml边界条件的交错网格高阶有限差分方法，对所述三层模型进行正演数值模拟；2.2.2)储层参数的影响分析在给定参数的基础上，保持其他参数不变，分别改变裂缝倾角、裂缝线密度、基质孔隙度和含气饱和度，以分析不同的参数对地震波场响应特征的影响。4.如权利要求2所述的基于norris-kg模型的地震纵波模拟分析方法，其特征在于，所述正演数值模拟是采用交错网格高阶有限差分法对介质的一阶速度-应力方程进行数值求
解；具体包括以下步骤：3.1)方程的离散化处理采用交错网格高阶有限差分算法，其时间、空间高阶差分近似公式分别为：采用交错网格高阶有限差分算法，其时间、空间高阶差分近似公式分别为：其中，上式(4)-(5)中δt为时间步长，2m为时间差分精度，2l为空间差分精度；δx为差分网格间距，为不同阶的交错网格空间差分系数，可以通过如下方程组求取:将该方程组带入到上述方程(3)中得到一阶速度应力方程的差分近似格式；3.2)pml边界条件处理采用pml边界条件处理方法，对上述方程(3)中的总体波场v进行变量分离，每一个波场变量可以分为纵波分量vp和横波分量vs两部分，t是总的应力场，上角标x、y、z代表三个坐标轴方向：然后对于所给定的三维空间模型设置边界衰减系数。5.如权利要求1所述的基于norris-kg模型的地震纵波模拟分析方法，其特征在于，所述norris-kg裂缝-孔隙等效介质模型的具体构建过程为：基于norris周期成层模型和kg裂缝-孔隙等效介质模型建立norris-kg等效介质模型；基于norris-kg等效介质模型，假设裂缝为非常薄且松软的无限大裂缝平面，将norris-kg等效介质模型看成是周期成层两套裂缝-孔隙地层，根据norris理论及heuristic假设，norris-kg等效介质模型的弹性属性表达式为：上式(1)中，为垂直裂缝面方向纵波模量，为模型高频极限弹性系数，为垂直裂缝面方向的高频极限弹性系数，为模型低频极限弹性系数；式(1)给出的是水平裂缝性储层的弹性系数矩阵，可计算任意入射角情况的频变纵波速度随流体饱和度的变化情况。

技术总结
本发明提供了一种基于Norris-KG模型的地震纵波模拟分析方法，其基于Norris-KG裂缝-孔隙等效介质模型，通过qP波相速度公式及纵、横波分离的一阶速度应力方程对介质的波场控制方程进行推导，并利用基于PML边界条件的交错网格高阶有限差分算法对地震波场进行正演数值模拟，最终，在给定参数下，对目的层的qP波相速度方位各向异性及其顶、底界面的反射波振幅方位各向异性分析，对目的层的各向异性及其影响参数进行优选分析,得到裂缝-孔隙介质地震纵波相速度和波场响应数据，最后进行纵波数值模拟分析。本发明构思合理，构思合理，能从地震纵波响应特征中直接对储层的物性参数进行识别，分析效果好且精确。分析效果好且精确。分析效果好且精确。


技术研发人员：周惠 刘海浩 孔丽云 康海霞 李娟 罗卫锋
受保护的技术使用者：中国地质调查局油气资源调查中心
技术研发日：2023.04.07
技术公布日：2023/8/9