一种采用地震数值模拟勘探天然气水合物的方法与流程

本发明实施例涉及天然气水合物资源勘查技术领域，具体涉及一种采用地震数值模拟勘探天然气水合物的方法。

背景技术：

天然气水合物为一种似冰状混合物，由水分子笼型结构和甲烷或其他低分子量气体组成。水合物中存储的大量天然气是重要的潜在能源，海域水合物向海水和大气中大量释放天然气，在过去的气候变化中也扮演着重要的角色。对天然气水合物资源的精确勘查是未来进行水合物规模化开采和温室气候研究等领域的基础。在海域水合物系统的地震成像中常可以看到似海底反射信号(bsr)。bsr为水合物-水层和下伏气-水层(足够高的气体浓度)的强烈波阻抗差异所致，通常与水合物稳定域的底界面一致。水合物在海底的埋藏具有多种分布模式，常具有强烈的非均质性。地震全波形反演(fwi)可直接得到地层地震参数的高分辨率成像，是天然气水合物资源勘查的重要方法。

然而，地震fwi需求大量的计算资源和时间，减少计算量的习惯做法是采用声波近似，即只反演介质的纵波速度。然而，这只在考虑小入射角和非转换纵波的情况下，且密度和横波速度剖面足够平滑时才有效。对海洋天然气水合物地震勘探来说，受地震反射振幅与偏移距关系影响，且考虑含气量对纵波速度的影响关系等因素，用声波近似反演地震数据的应用效果并不够理想。弹性波fwi算法，特别是在频率域，可以考虑到纵波、横波、转换波、首波和瑞雷波等，只用有限频点数的数据就能很好的恢复地下的模型参数结构，具有天然的多尺度性。在宽方位、多分量地震数据采集处理和解释技术的发展带动下，反演中可用的剖面越来越多，信噪比也逐渐提高，地震波对地下的照明角度也逐渐增加，而照明角度增加意味着可以得到更加完整的水合物地层信息。

频率域弹性波fwi强烈依赖于频率域正演，有限差分方法是最重要的正演方法之一。需要对频率域弹性波方程进行有限差分离散，根据模型大小和计算频率来综合考虑频散关系，然后选取合适的差分阶数并推导差分格式，最后带入模型参数构建大型稀疏方程组进行求解。传统的频率域有限差分算法精度低，容易发生网格频散，严重影响计算的效率。

技术实现要素：

为此，本发明实施例提供一种采用地震数值模拟勘探天然气水合物的方法，将得到的参数变换为频率域波动方程，并基于此进行最优化变网格有限差分格式的建立，进一步通过该最优化差分格式构建大型稀疏方程组以实现海洋天然气水合物储层模型弹性波响应的高效求解，从而有效驱动频率域全波形反演，为天然气水合物资源的勘查的效率提高提供有效的技术保障。

为了实现上述目的，本发明的实施方式提供如下技术方案：

在本发明实施例的一个方面，提供了一种采用地震数值模拟勘探天然气水合物的方法，其特征在于，包括：

s100、在二维情形下，测定天然气水合物地震勘探模型的参数，并根据得到的参数设定原始的时间域弹性波动方程；

s200、将得到的时间域弹性波动方程变换为频率域波动方程；

s300、基于交错网格，消除应力分量，对频率域波动方程进行空间离散，得到四阶精度条件下的差分格式；

s400、对得到的差分格式进行推导，得到频散关系；

s500、对得到的频散关系采用高斯牛顿方法进行最优化，得到一系列最优化差分系数，并将得到的最优化差分系数带入步骤s300中的差分格式中，重新构建形成最优化交错网格有限差分格式；

s600、将模型参数带入步骤s500中得到的最优化交错网格有限差分格式，构建大型稀疏方程组，并进行求解，得到天然气水合物的高精度弹性波响应。

作为本发明的一种优选方案，步骤s100中，所述时间域波动方程如式i所示，且所述式i具体包括：

式ia:

式ib:

式ic:

式id:

式ie:

其中，ρ为密度，x表示空间上的水平方向、z表示空间上的竖直方向，vx和vz分别表示水平方向和竖直方向上的速度分量，σxx和σzz为正应力分量，σxz为剪应力分量，c11、c33和c13为弹性模量；其中，在各向同性介质中，c13＝c31＝c11-2c55，vp和vs分别代表介质的纵波速度和横波速度。

作为本发明的一种优选方案，步骤s200中，将时间域弹性波动方程变换为频率域波动方程为采用傅里叶变换得到，且得到的频率域波动方程如式ii所示，所述式ii具体包括：

式iia：

式iib：

式iic：

式iid：

式iie：

其中，i²＝-1，ω＝2πf，f为频率。

作为本发明的一种优选方案，步骤s300中得到的差分格式如式iii所示，且所述式iii具体包括：

式iiia：

式iiib：

其中，

且，sx和sz分别代表x方向和z方向上与吸收边界条件相关的表达式，c1、c2、c3和c4为4个待优化的差分系数，上标a和b分别表示变网格剖分引起的两种差分表示类型。

作为本发明的一种优选方案，步骤s400中的频散关系如式iv所示，且式iv具体包括：

式iva：

式ivb：

其中，αph和βph分别表示纵横波相速度，k为波数矢量，r为纵横波速度比，a、b和c分别为与入射角、差分系数c、质量加权系数a以及网格间距相关的表达式。

作为本发明的一种优选方案，步骤s500中最优化差分系数的计算具体包括：

s501、对得到的频散关系建立泛函，得到如式v所示的公式：

式v：

s502、选择a1＝1.0、a2＝0、a3＝0、a4＝0、c1＝1.125和c2＝-0.042作为初始值，采用高斯牛顿法求最小值，得到一系列最优化系数。

作为本发明的一种优选方案，单位横波波长内，交错网格的网格数量为3个。

作为本发明的一种优选方案，天然气水合物地震勘探模型的参数的测定包括：

s101、选定地震勘探模型的待测定区域；

s102、将所述待测定区域中均匀划分为多个激发区域，并将每个激发区域的中点设置为激发点，并在每个激发点上设置炮点；

s103、在每个激发区域中根据设定的检测间距，设置多个检波点；

s104、顺次激发炮点，并收集每个检波点数据；

s105、将检波点的数据进行收集，并根据对应的激发区域将每个激发区域中的数据分类收集；

s106、对收集的数据进行筛分过滤，滤除无效信息，整合有效信息，得到参数。

作为本发明的一种优选方案，步骤s106具体包括：

s1061、对每个激发区域中的数据分别根据对应激发区域中的炮点激发时间筛选出对应时间段内的有效数据，并记为数据a；

s1062、对得到的数据a中以对应的炮点激发时间为起始点向后计算，总时长满足预设时长的数据进行收集，并记为数据b；

s1063、将步骤s1062中筛除的数据记为数据c，将4-8个相邻的检波点作为一个收集区域进行划分，将每个激发区域分别划分为多个收集区域，对每个不含有数据b的收集区域中的数据c进行计数，当某个所述收集区域中的数据c的数量大于2个，则将该收集区域作为整体计数单元，并对其中的多个数据c的数据以时间为轴进行整合，整合后的数据记为数据d；

s1064、当数据d的总时长满足预设时长时，则对该数据d进行收集；

s1065、将收集到的数据b和数据d作为天然气水合物地震勘探模型的参数。

作为本发明的一种优选方案，预设时长为1-2s，且数据采样过程中时间间隔为0.5-1.5ms。

本发明的实施方式具有如下优点：

1、本发明基于交错网格的差分格式不受介质泊松比影响，能适应含流体介质模型的弹性波模拟条件，对海洋天然气水合物资源勘查具有明显的优势；

2、本发明通过推导出频率域弹性波交错网格有限差分格式的频散关系，并对查分系数进行最优化，使最终得到的差分算子模拟精度高，在单位横波博场内仅需3个网格，极大地提高了天然气水合物资源勘查中弹性波数值模拟的效率；

3、本发明推导得到基于变网格条件下的频率域弹性波四阶精度交错网格有限差分格式，有效提高了建模的灵活性，进一步提高了弹性波数值模拟的效率。

附图说明

为了更清楚地说明本发明的实施方式或现有技术中的技术方案，下面将对实施方式或现有技术描述中所需要使用的附图作简单地介绍。显而易见地，下面描述中的附图仅仅是示例性的，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图引伸获得其它的实施附图。

本说明书所绘示的结构、比例、大小等，均仅用以配合说明书所揭示的内容，以供熟悉此技术的人士了解与阅读，并非用以限定本发明可实施的限定条件，故不具技术上的实质意义，任何结构的修饰、比例关系的改变或大小的调整，在不影响本发明所能产生的功效及所能达成的目的下，均应仍落在本发明所揭示的技术内容得能涵盖的范围内。

图1为本发明实施例中一种采用地震数值模拟勘探天然气水合物的方法的流程图；

图2为本发明实施例中二维平面内弹性有限差分的网格配置图；

图3为本发明实施例中用于空间变网格的网格点计算(a)整节点和(b)相邻整节点之间的差分；

图4为本发明实施例中介质泊松比在0.48时最优化差分格式的频散曲线；

图5为本发明实施例中含流体介质的纵波速度模型；

图6为本发明实施例中常规有限差分方法与本发明有限差分方法计算含流体介质模型5hz的波场分量；

图7为本发明实施例中常规有限差分方法与本发明有限差分方法计算含流体介质模型30hz的波场分量。

具体实施方式

以下由特定的具体实施例说明本发明的实施方式，熟悉此技术的人士可由本说明书所揭露的内容轻易地了解本发明的其他优点及功效，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明中的最优化交错网格有限差分格式基于变网格有限差分离散方法，推导出频率域弹性波方程的差分格式以及频散关系，通过对得到的差分格式进行高斯牛顿优化，给出一系列最优化差分系数，通过将差分系数带入到差分格式并构建大型稀疏方程组，最终实现弹性波响应的高效精确求解。具体地，数学推演过程如下：

步骤1、对天然气水合物地震勘探模型，考虑二维情形，设置原始的时间域弹性波动方程：

式ia:

式ib:

式ic:

式id:

式ie:

其中，ρ为介质的密度，x、y和z分别表示空间的三个方向，对应地，x表示空间上的水平方向、z表示空间上的竖直方向，vx和vz分别表示水平方向和竖直方向上的速度分量，σxx和σzz为正应力分量，σxz为剪应力分量，c11、c33和c13为弹性模量；其中，在各向同性介质中，c13＝c31＝c11-2c55，vp和vs分别代表介质的纵波速度和横波速度。

步骤2、采用傅里叶变换，将时间域弹性波动方程变换到频率域，对应得到频率域波动方程：

式iia：

式iib：

式iic：

式iid：

式iie：

其中，i²＝-1，ω＝2πf，f为频率。

步骤3、对式iia-式iie的求解可以在经典的交错网格中进行，具体如图2所示，主应力分配在网格的顶点中，剪应力分配在网格中心，而速度场则分配在网格边界的中点。交错网格有限差分法对介质泊松比变化剧烈的地层有很好的适应性，因此是时间域模拟复杂介质地震响应的重要方法，但在频率域的弹性波应用上还不是很广泛。在此基础上，由于式ii中既有速度场也有应力场，这样将导致自由度的成倍增长，如果以这样的形式直接求解会引起巨大的不必要浪费，因此需要将其改造成只包含速度场分量的表达式，从而进行求解。为了更好的示意推导过程，以二阶精度为例推导非均匀各向同性介质情况下二维弹性波方程，四阶精度可按照同样的过程进行推导。具体的推导过程如下：

考虑到各向同性介质，采用吸收边界条件来吸收来自人工边界的虚假反射，对x和z方向，其表达式分别记为sx和sz，令代入式iia可得：

由二阶中心差分格式得到应力和应变的偏导数近似为：

同样地，对式iic和式iie进行推导可得：

速度场偏导数的二阶差分近似为：

对上述表达式进行层层带入，即首先写出不同位置应力场的表达式，再将这些表达式带入到应力场的偏导数表达式中，最后将结果带入到速度表达式中即可消除所有的应力项，形成的有限差分格式中只含有速度分量。

本发明为了提高模拟的精度和效率，采用四阶精度进行离散推导，结合变网格离散技术来提高建模的灵活性。由变网格引起的两种有限差分不同表示方法如图3所示。以vx为例，四阶精度有限差分格式可以写成如下的形式：

经过推导和整理，步骤3中四阶精度条件下的差分格式如下所示：

式iiia：

式iiib：

其中，

且，sx和sz分别代表x方向和z方向上与吸收边界条件相关的表达式，c1、c2、c3和c4为4个待优化的差分系数，上标a和b分别表示变网格剖分引起的两种差分表示类型。

步骤4、设在全空间均匀介质中的平面谐波为v^ω＝v^ωe^ik·r，其中v＝[vxvz]^t为水平速度分量和垂直速度分量的幅度，k＝[kxkz]^t为波数向量，r＝[xz]^t为位置向量，将v^ω代入步骤3中得到的差分格式中，得到：

其中，m表示质量加权平均算子，α和β为纵横波速度，fxx、fxz、fzx、fzz以及fx和fz分别表示空间导数的有限差分算子。在均匀各向同性介质情况下fx＝fz、fxx＝fzz和fxz＝fzx。考虑均匀网格条件下和得到如下表达式：

m＝a1+a2[cos(kx)+cos(kz)]+a3[cos(2kx)+cos(2kz)]+a4[cos(3kx)+cos(3kz)]+a5[cos(kx+kz)+cos(kx-kz)]+a6[cos(2kx+2kz)+cos(2kx-2kz)]+a7[cos(2kx+kz)+cos(2kx-kz)+cos(kx+2kz)+cos(kx-2kz)]

其中，kx＝2πkcos(θ)、kz＝2πksin(θ)分别为每个网格点的波数，这里为了分析方便，采用了正方形网格剖分。为了使频散方程有解，则必须满足左边的矩阵行列式为零，对纵、横波，通过相速度vph＝ω/k，群速度得到频散关系为：

式iva：

式ivb：

其中，a＝m²，b＝m(1+r²)(fxx+fxz)，和c＝(fxx+r²fxz)(r²fxx+fxz)-m²(1-r²)²。

步骤5、对得到的频散关系建立泛函，得到：

式v：

设定初始值：a1＝1.0、a2＝0、a3＝0、a4＝0、c1＝1.125和c2＝-0.042。利用高斯牛顿法求极小值，即可得到一系列最优化系数。具体地，本发明最终给出的参考最优化系数为：a1＝0.7527、a2＝0.2031、a3＝-0.112、a4＝0.0276、c1＝1.08和c2＝-0.049。图4为介质泊松比在0.48时最优化差分格式的频散曲线。可见当经过最优化处理后可以实现高精度的模拟，当单位波长内不小于2.9个网格的时候可将误差控制在2％以内。

步骤6、将模型参数带入步骤5得出的差分格式，构建大型稀疏方程组，并进行求解，最终得到天然气水合物的高精度弹性波响应。作为示例，设计了图5所示的模型，固体背景介质的纵、横波速度分别为3000m/s、1800m/s，流体块的纵横波速度分别为1500m/s和0，固体和流体介质的密度分别为2000kg/m³、1000kg/m³。图6和图7分别显示了采用爆炸震源激发得到的5hz和30hz弹性波场(实部)结果。图中可以看到，由于常规网格对于流固耦合界面的适应性差，即使在5hz低频情况下，在流体块内部也出现明显的频散。而采用交错网格模拟的结果在波场图上的流体块内基本没有表现频散特征。对30hz激发时的模拟结果，用常规网格模拟的结果出现严重的混乱，而交错网格的模拟结果则呈现出稳定的波场，在满足频散条件内给出了正确的模拟结果。其中，图6和图7中a和b为采用常规有限差分方法对应得到的波场分量，c和d为采用本发明的方法对应得到的波场分量。

进一步，为了更好地提高后期数据分析的精确度，对于前期的模型勘探的参数的测定可以进一步采用如下方法：

步骤101、根据需要测定的区域和区域范围实际大小以及测试区能源分布情况，选定地震勘探模型的待测定区域；

步骤102、将所述待测定区域中均匀划分为多个激发区域，并将每个激发区域的中点设置为激发点，并在每个激发点上设置炮点；

步骤103、在每个激发区域中根据设定的检测间距(具体地，相邻的两个检波点的距离一般在20-50m)，设置多个检波点；

步骤104、顺次激发炮点，并收集每个检波点数据；

步骤105、将检波点的数据进行收集，并根据对应的激发区域将每个激发区域中的数据分类收集；

步骤106、对每个激发区域中的数据分别根据对应激发区域中的炮点激发时间筛选出对应时间段内的有效数据，并记为数据a；对得到的数据a中以对应的炮点激发时间为起始点向后计算，总时长满足预设时长的数据进行收集，并记为数据b；将上述步骤中筛除的数据记为数据c，将4-8个相邻的检波点作为一个收集区域进行划分，将每个激发区域分别划分为多个收集区域，对每个不含有数据b的收集区域中的数据c进行计数，当某个所述收集区域中的数据c的数量大于2个，则将该收集区域作为整体计数单元，并对其中的多个数据c的数据以时间为轴进行整合，整合后的数据记为数据d；当数据d的总时长满足预设时长时，则对该数据d进行收集；将收集到的数据b和数据d作为天然气水合物地震勘探模型的参数。这样的方式能够将有效数据和相对无效数据均进行分类收集，不仅降低了对数据采集设备数量等要求，降低了采集成本，且对数据信息的收集更为全面，能够更好地提高数据收集的有效性和覆盖范围的广泛性。

虽然，上文中已经用一般性说明及具体实施例对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。