CN112526605B - 一种采用地震数值模拟勘探天然气水合物的方法 - Google Patents

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

Info

Publication number
CN112526605B
CN112526605B CN202011545296.7A CN202011545296A CN112526605B CN 112526605 B CN112526605 B CN 112526605B CN 202011545296 A CN202011545296 A CN 202011545296A CN 112526605 B CN112526605 B CN 112526605B
Authority
CN
China
Prior art keywords
data
formula
optimized
excitation
difference
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202011545296.7A
Other languages
English (en)
Other versions
CN112526605A (zh
Inventor
马超
高妍
陆程
沙志斌
耿澜涛
边航
齐荣荣
邢东辉
孟凡乐
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Guangzhou Marine Geological Survey
Original Assignee
Guangzhou Marine Geological Survey
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Guangzhou Marine Geological Survey filed Critical Guangzhou Marine Geological Survey
Priority to CN202011545296.7A priority Critical patent/CN112526605B/zh
Publication of CN112526605A publication Critical patent/CN112526605A/zh
Application granted granted Critical
Publication of CN112526605B publication Critical patent/CN112526605B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30Analysis
    • G01V1/306Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/62Physical property of subsurface
    • G01V2210/624Reservoir parameters
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明实施例公开了一种采用地震数值模拟勘探天然气水合物的方法,包括:测定模型参数,根据参数设定时间域弹性波动方程;将时间域弹性波动方程变换为频率域波动方程;基于交错网格,消除应力分量,得到四阶精度条件下的差分格式;推导得到频散关系;对频散关系最优化,得到最优化差分系数,并带入差分格式中,形成最优化交错网格有限差分格式;将模型参数带入最优化交错网格有限差分格式,构建大型稀疏方程组,求解,得到天然气水合物的高精度弹性波响应。本发明基于变网格有限差分离散方法,推导出差分格式以及频散关系,并对差分格式进行高斯牛顿优化,通过将优化的差分系数重新代入并构建大型稀疏方程组,实现弹性波响应的高效精确求解。

Description

一种采用地震数值模拟勘探天然气水合物的方法
技术领域
本发明实施例涉及天然气水合物资源勘查技术领域,具体涉及一种采用地震数值模拟勘探天然气水合物的方法。
背景技术
天然气水合物为一种似冰状混合物,由水分子笼型结构和甲烷或其他低分子量气体组成。水合物中存储的大量天然气是重要的潜在能源,海域水合物向海水和大气中大量释放天然气,在过去的气候变化中也扮演着重要的角色。对天然气水合物资源的精确勘查是未来进行水合物规模化开采和温室气候研究等领域的基础。在海域水合物系统的地震成像中常可以看到似海底反射信号(BSR)。BSR为水合物-水层和下伏气-水层(足够高的气体浓度)的强烈波阻抗差异所致,通常与水合物稳定域的底界面一致。水合物在海底的埋藏具有多种分布模式,常具有强烈的非均质性。地震全波形反演(FWI)可直接得到地层地震参数的高分辨率成像,是天然气水合物资源勘查的重要方法。
然而,地震FWI需求大量的计算资源和时间,减少计算量的习惯做法是采用声波近似,即只反演介质的纵波速度。然而,这只在考虑小入射角和非转换纵波的情况下,且密度和横波速度剖面足够平滑时才有效。对海洋天然气水合物地震勘探来说,受地震反射振幅与偏移距关系影响,且考虑含气量对纵波速度的影响关系等因素,用声波近似反演地震数据的应用效果并不够理想。弹性波FWI算法,特别是在频率域,可以考虑到纵波、横波、转换波、首波和瑞雷波等,只用有限频点数的数据就能很好的恢复地下的模型参数结构,具有天然的多尺度性。在宽方位、多分量地震数据采集处理和解释技术的发展带动下,反演中可用的剖面越来越多,信噪比也逐渐提高,地震波对地下的照明角度也逐渐增加,而照明角度增加意味着可以得到更加完整的水合物地层信息。
频率域弹性波FWI强烈依赖于频率域正演,有限差分方法是最重要的正演方法之一。需要对频率域弹性波方程进行有限差分离散,根据模型大小和计算频率来综合考虑频散关系,然后选取合适的差分阶数并推导差分格式,最后带入模型参数构建大型稀疏方程组进行求解。传统的频率域有限差分算法精度低,容易发生网格频散,严重影响计算的效率。
发明内容
为此,本发明实施例提供一种采用地震数值模拟勘探天然气水合物的方法,将得到的参数变换为频率域波动方程,并基于此进行最优化变网格有限差分格式的建立,进一步通过该最优化差分格式构建大型稀疏方程组以实现海洋天然气水合物储层模型弹性波响应的高效求解,从而有效驱动频率域全波形反演,为天然气水合物资源的勘查的效率提高提供有效的技术保障。
为了实现上述目的,本发明的实施方式提供如下技术方案:
在本发明实施例的一个方面,提供了一种采用地震数值模拟勘探天然气水合物的方法,其特征在于,包括:
S100、在二维情形下,测定天然气水合物地震勘探模型的参数,并根据得到的参数设定原始的时间域弹性波动方程;
S200、将得到的时间域弹性波动方程变换为频率域波动方程;
S300、基于交错网格,消除应力分量,对频率域波动方程进行空间离散,得到四阶精度条件下的差分格式;
S400、对得到的差分格式进行推导,得到频散关系;
S500、对得到的频散关系采用高斯牛顿方法进行最优化,得到一系列最优化差分系数,并将得到的最优化差分系数带入步骤S300中的差分格式中,重新构建形成最优化交错网格有限差分格式;
S600、将模型参数带入步骤S500中得到的最优化交错网格有限差分格式,构建大型稀疏方程组,并进行求解,得到天然气水合物的高精度弹性波响应。
作为本发明的一种优选方案,步骤S100中,所述时间域波动方程如式I所示,且所述式I具体包括:
式IA:
Figure BDA0002855794480000031
式IB:
Figure BDA0002855794480000032
式IC:
Figure BDA0002855794480000033
式ID:
Figure BDA0002855794480000034
式IE:
Figure BDA0002855794480000035
其中,ρ为密度,x表示空间上的水平方向、z表示空间上的竖直方向,vx和vz分别表示水平方向和竖直方向上的速度分量,σxx和σzz为正应力分量,σxz为剪应力分量,c11、c33和c13为弹性模量;其中,在各向同性介质中,
Figure BDA0002855794480000036
c13=c31=c11-2c55,VP和VS分别代表介质的纵波速度和横波速度。
作为本发明的一种优选方案,步骤S200中,将时间域弹性波动方程变换为频率域波动方程为采用傅里叶变换得到,且得到的频率域波动方程如式II所示,所述式II具体包括:
式IIA:
Figure BDA0002855794480000037
式IIB:
Figure BDA0002855794480000038
式IIC:
Figure BDA0002855794480000039
式IID:
Figure BDA00028557944800000310
式IIE:
Figure BDA0002855794480000041
其中,i2=-1,ω=2πf,f为频率。
作为本发明的一种优选方案,步骤S300中得到的差分格式如式III所示,且所述式III具体包括:
式IIIA:
Figure BDA0002855794480000042
式IIIB:
Figure BDA0002855794480000043
其中,
Figure BDA0002855794480000044
Figure BDA0002855794480000045
Figure BDA0002855794480000046
Figure BDA0002855794480000047
Figure BDA0002855794480000048
Figure BDA0002855794480000049
Figure BDA00028557944800000410
Figure BDA00028557944800000411
Figure BDA00028557944800000412
Figure BDA0002855794480000051
Figure BDA0002855794480000052
Figure BDA0002855794480000053
Figure BDA0002855794480000054
Figure BDA0002855794480000055
Figure BDA0002855794480000056
Figure BDA0002855794480000057
Figure BDA0002855794480000058
Figure BDA0002855794480000059
Figure BDA00028557944800000510
Figure BDA00028557944800000511
Figure BDA00028557944800000512
Figure BDA00028557944800000513
Figure BDA00028557944800000514
Figure BDA0002855794480000061
Figure BDA0002855794480000062
Figure BDA0002855794480000063
Figure BDA0002855794480000064
Figure BDA0002855794480000065
Figure BDA0002855794480000066
Figure BDA0002855794480000067
Figure BDA0002855794480000068
Figure BDA0002855794480000069
且,sx和sz分别代表x方向和z方向上与吸收边界条件相关的表达式,c1、c2、c3和c4为4个待优化的差分系数,上标a和b分别表示变网格剖分引起的两种差分表示类型。
作为本发明的一种优选方案,步骤S400中的频散关系如式IV所示,且式IV具体包括:
式IVA:
Figure BDA00028557944800000610
式IVB:
Figure BDA00028557944800000611
其中,αph和βph分别表示纵横波相速度,K为波数矢量,R为纵横波速度比,A、B和C分别为与入射角、差分系数c、质量加权系数a以及网格间距相关的表达式。
作为本发明的一种优选方案,步骤S500中最优化差分系数的计算具体包括:
S501、对得到的频散关系建立泛函,得到如式V所示的公式:
式V:
Figure BDA0002855794480000071
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:
Figure BDA0002855794480000101
式IB:
Figure BDA0002855794480000102
式IC:
Figure BDA0002855794480000103
式ID:
Figure BDA0002855794480000104
式IE:
Figure BDA0002855794480000105
其中,ρ为介质的密度,x、y和z分别表示空间的三个方向,对应地,x表示空间上的水平方向、z表示空间上的竖直方向,vx和vz分别表示水平方向和竖直方向上的速度分量,σxx和σzz为正应力分量,σxz为剪应力分量,c11、c33和c13为弹性模量;其中,在各向同性介质中,
Figure BDA0002855794480000106
Figure BDA0002855794480000107
c13=c31=c11-2c55,VP和VS分别代表介质的纵波速度和横波速度。
步骤2、采用傅里叶变换,将时间域弹性波动方程变换到频率域,对应得到频率域波动方程:
式IIA:
Figure BDA0002855794480000108
式IIB:
Figure BDA0002855794480000111
式IIC:
Figure BDA0002855794480000112
式IID:
Figure BDA0002855794480000113
式IIE:
Figure BDA0002855794480000114
其中,i2=-1,ω=2πf,f为频率。
步骤3、对式IIA-式IIE的求解可以在经典的交错网格中进行,具体如图2所示,主应力分配在网格的顶点中,剪应力分配在网格中心,而速度场则分配在网格边界的中点。交错网格有限差分法对介质泊松比变化剧烈的地层有很好的适应性,因此是时间域模拟复杂介质地震响应的重要方法,但在频率域的弹性波应用上还不是很广泛。在此基础上,由于式II中既有速度场也有应力场,这样将导致自由度的成倍增长,如果以这样的形式直接求解会引起巨大的不必要浪费,因此需要将其改造成只包含速度场分量的表达式,从而进行求解。为了更好的示意推导过程,以二阶精度为例推导非均匀各向同性介质情况下二维弹性波方程,四阶精度可按照同样的过程进行推导。具体的推导过程如下:
考虑到各向同性介质,采用吸收边界条件来吸收来自人工边界的虚假反射,对x和z方向,其表达式分别记为Sx和Sz,令
Figure BDA0002855794480000115
代入式IIA可得:
Figure BDA0002855794480000116
由二阶中心差分格式得到应力和应变的偏导数近似为:
Figure BDA0002855794480000117
Figure BDA0002855794480000118
同样地,对式IIC和式IIE进行推导可得:
Figure BDA0002855794480000121
Figure BDA0002855794480000122
速度场偏导数的二阶差分近似为:
Figure BDA0002855794480000123
Figure BDA0002855794480000124
Figure BDA0002855794480000125
Figure BDA0002855794480000126
对上述表达式进行层层带入,即首先写出不同位置应力场的表达式,再将这些表达式带入到应力场的偏导数表达式中,最后将结果带入到速度表达式中即可消除所有的应力项,形成的有限差分格式中只含有速度分量。
本发明为了提高模拟的精度和效率,采用四阶精度进行离散推导,结合变网格离散技术来提高建模的灵活性。由变网格引起的两种有限差分不同表示方法如图3所示。以vx为例,四阶精度有限差分格式可以写成如下的形式:
Figure BDA0002855794480000127
经过推导和整理,步骤3中四阶精度条件下的差分格式如下所示:
式IIIA:
Figure BDA0002855794480000128
式IIIB:
Figure BDA0002855794480000129
其中,
Figure BDA0002855794480000131
Figure BDA0002855794480000132
Figure BDA0002855794480000133
Figure BDA0002855794480000134
Figure BDA0002855794480000135
Figure BDA0002855794480000136
Figure BDA0002855794480000137
Figure BDA0002855794480000138
Figure BDA0002855794480000139
Figure BDA00028557944800001310
Figure BDA00028557944800001311
Figure BDA00028557944800001312
Figure BDA00028557944800001313
Figure BDA00028557944800001314
Figure BDA0002855794480000141
Figure BDA0002855794480000142
Figure BDA0002855794480000143
Figure BDA0002855794480000144
Figure BDA0002855794480000145
Figure BDA0002855794480000146
Figure BDA0002855794480000147
Figure BDA0002855794480000148
Figure BDA0002855794480000149
Figure BDA00028557944800001410
Figure BDA00028557944800001411
Figure BDA00028557944800001412
Figure BDA00028557944800001413
Figure BDA00028557944800001414
Figure BDA0002855794480000151
Figure BDA0002855794480000152
Figure BDA0002855794480000153
Figure BDA0002855794480000154
且,sx和sz分别代表x方向和z方向上与吸收边界条件相关的表达式,c1、c2、c3和c4为4个待优化的差分系数,上标a和b分别表示变网格剖分引起的两种差分表示类型。
步骤4、设在全空间均匀介质中的平面谐波为vω=Vωeik·r,其中V=[Vx Vz]T为水平速度分量和垂直速度分量的幅度,k=[kx kz]T为波数向量,r=[x z]T为位置向量,将vω代入步骤3中得到的差分格式中,得到:
Figure BDA0002855794480000155
其中,M表示质量加权平均算子,α和β为纵横波速度,FXX、FXZ、FZX、FZZ以及FX和FZ分别表示空间导数的有限差分算子。在均匀各向同性介质情况下FX=FZ、FXX=FZZ和FXZ=FZX。考虑均匀网格条件下
Figure BDA0002855794480000156
Figure BDA0002855794480000157
得到如下表达式:
Figure BDA0002855794480000158
Figure BDA0002855794480000159
Figure BDA00028557944800001510
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,群速度
Figure BDA0002855794480000161
得到频散关系为:
式IVA:
Figure BDA0002855794480000162
式IVB:
Figure BDA0002855794480000163
其中,A=M2,B=M(1+R2)(FXX+FXZ),
Figure BDA0002855794480000164
和C=(FXX+R2FXZ)(R2FXX+FXZ)-M2(1-R2)2
步骤5、对得到的频散关系建立泛函,得到:
式V:
Figure BDA0002855794480000165
设定初始值: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/m3、1000kg/m3。图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作为天然气水合物地震勘探模型的参数。这样的方式能够将有效数据和相对无效数据均进行分类收集,不仅降低了对数据采集设备数量等要求,降低了采集成本,且对数据信息的收集更为全面,能够更好地提高数据收集的有效性和覆盖范围的广泛性。
虽然,上文中已经用一般性说明及具体实施例对本发明作了详尽的描述,但在本发明基础上,可以对之作一些修改或改进,这对本领域技术人员而言是显而易见的。因此,在不偏离本发明精神的基础上所做的这些修改或改进,均属于本发明要求保护的范围。

Claims (5)

1.一种采用地震数值模拟勘探天然气水合物的方法,其特征在于,包括:
S100、在二维情形下,测定天然气水合物地震勘探模型的参数,并根据得到的参数设定原始的时间域弹性波动方程;
S200、将得到的时间域弹性波动方程变换为频率域波动方程;
S300、基于交错网格,消除应力分量,对频率域波动方程进行空间离散,得到四阶精度条件下的差分格式;
S400、对得到的差分格式进行推导,得到频散关系;
S500、对得到的频散关系采用高斯牛顿方法进行最优化,得到一系列最优化差分系数,并将得到的最优化差分系数带入步骤S300中的差分格式中,重新构建形成最优化交错网格有限差分格式;
S600、将模型参数带入步骤S500中得到的最优化交错网格有限差分格式,构建大型稀疏方程组,并进行求解,得到天然气水合物的高精度弹性波响应;
其中:
步骤S100中,所述时间域弹性波动方程如式I所示,且所述式I具体包括:
式IA:
Figure FDA0003728352440000011
式IB:
Figure FDA0003728352440000012
式IC:
Figure FDA0003728352440000013
式ID:
Figure FDA0003728352440000014
式IE:
Figure FDA0003728352440000015
其中,ρ为密度,x表示空间上的水平方向、z表示空间上的竖直方向,vx和vz分别表示水平方向和竖直方向上的速度分量,σxx和σzz为正应力分量,σxz为剪应力分量,c11、c33和c13为弹性模量;其中,在各向同性介质中,
Figure FDA0003728352440000021
c13=c31=c11-2c55,VP和VS分别代表介质的纵波速度和横波速度;
步骤S200中,将时间域弹性波动方程变换为频率域波动方程为采用傅里叶变换得到,且得到的频率域波动方程如式II所示,所述式II具体包括:
式IIA:
Figure FDA0003728352440000022
式IIB:
Figure FDA0003728352440000023
式IIC:
Figure FDA0003728352440000024
式IID:
Figure FDA0003728352440000025
式IIE:
Figure FDA0003728352440000026
其中,i2=-1,ω=2πf,f为频率;
步骤S300中得到的差分格式如式III所示,且所述式III具体包括:
式IIIA:
Figure FDA0003728352440000027
式IIIB:
Figure FDA0003728352440000028
其中,
Figure FDA0003728352440000029
Figure FDA0003728352440000031
Figure FDA0003728352440000032
Figure FDA0003728352440000033
Figure FDA0003728352440000034
Figure FDA0003728352440000035
Figure FDA0003728352440000036
Figure FDA0003728352440000037
Figure FDA0003728352440000038
Figure FDA0003728352440000039
Figure FDA00037283524400000310
Figure FDA00037283524400000311
Figure FDA00037283524400000312
Figure FDA00037283524400000313
Figure FDA00037283524400000314
Figure FDA0003728352440000041
Figure FDA0003728352440000042
Figure FDA0003728352440000043
Figure FDA0003728352440000044
Figure FDA0003728352440000045
Figure FDA0003728352440000046
Figure FDA0003728352440000047
Figure FDA0003728352440000048
Figure FDA0003728352440000049
Figure FDA00037283524400000410
Figure FDA00037283524400000411
Figure FDA00037283524400000412
Figure FDA00037283524400000413
Figure FDA00037283524400000414
Figure FDA0003728352440000051
Figure FDA0003728352440000052
Figure FDA0003728352440000053
且,sx和sz分别代表x方向和z方向上与吸收边界条件相关的表达式,c1、c2、c3和c4为4个待优化的差分系数,上标a和b分别表示变网格剖分引起的两种差分表示类型;
步骤S400中的频散关系如式IV所示,且式IV具体包括:
式IVA:
Figure FDA0003728352440000054
式IVB:
Figure FDA0003728352440000055
其中,αph和βph分别表示纵横波相速度,K为波数矢量,R为纵横波速度比,A、B和C分别为与入射角、差分系数c、质量加权系数a以及网格间距相关的表达式;
步骤S500中最优化差分系数的计算具体包括:
S501、对得到的频散关系建立泛函,得到如式V所示的公式:
式V:
Figure FDA0003728352440000056
S502、选择a1=1.0、a2=0、a3=0、a4=0、c1=1.125和c2=-0.042作为初始值,采用高斯牛顿法求最小值,得到一系列最优化系数。
2.根据权利要求1所述的方法,其特征在于,单位横波波长内,交错网格的网格数量为3个。
3.根据权利要求1所述的方法,其特征在于,天然气水合物地震勘探模型的参数的测定包括:
S101、选定地震勘探模型的待测定区域;
S102、将所述待测定区域中均匀划分为多个激发区域,并将每个激发区域的中点设置为激发点,并在每个激发点上设置炮点;
S103、在每个激发区域中根据设定的检测间距,设置多个检波点;
S104、顺次激发炮点,并收集每个检波点数据;
S105、将检波点的数据进行收集,并根据对应的激发区域将每个激发区域中的数据分类收集;
S106、对收集的数据进行筛分过滤,滤除无效信息,整合有效信息,得到参数。
4.根据权利要求3所述的方法,其特征在于,步骤S106具体包括:
S1061、对每个激发区域中的数据分别根据对应激发区域中的炮点激发时间筛选出对应时间段内的有效数据,并记为数据A;
S1062、对得到的数据A中以对应的炮点激发时间为起始点向后计算,总时长满足预设时长的数据进行收集,并记为数据B;
S1063、将步骤S1062中筛除的数据记为数据C,将4-8个相邻的检波点作为一个收集区域进行划分,将每个激发区域分别划分为多个收集区域,对每个不含有数据B的收集区域中的数据C进行计数,当某个所述收集区域中的数据C的数量大于2个,则将该收集区域作为整体计数单元,并对其中的多个数据C的数据以时间为轴进行整合,整合后的数据记为数据D;
S1064、当数据D的总时长满足预设时长时,则对该数据D进行收集;
S1065、将收集到的数据B和数据D作为天然气水合物地震勘探模型的参数。
5.根据权利要求4所述的方法,其特征在于,预设时长为1-2s,且数据采样过程中时间间隔为0.5-1.5ms。
CN202011545296.7A 2020-12-24 2020-12-24 一种采用地震数值模拟勘探天然气水合物的方法 Active CN112526605B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011545296.7A CN112526605B (zh) 2020-12-24 2020-12-24 一种采用地震数值模拟勘探天然气水合物的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011545296.7A CN112526605B (zh) 2020-12-24 2020-12-24 一种采用地震数值模拟勘探天然气水合物的方法

Publications (2)

Publication Number Publication Date
CN112526605A CN112526605A (zh) 2021-03-19
CN112526605B true CN112526605B (zh) 2022-09-02

Family

ID=74976529

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011545296.7A Active CN112526605B (zh) 2020-12-24 2020-12-24 一种采用地震数值模拟勘探天然气水合物的方法

Country Status (1)

Country Link
CN (1) CN112526605B (zh)

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103149585B (zh) * 2013-01-30 2016-02-17 中国石油天然气集团公司 一种弹性偏移地震波场构建方法及装置
WO2016108896A1 (en) * 2014-12-31 2016-07-07 Landmark Graphics Corporation Seismic elastic wave simulation for tilted transversely isotropic media using adaptive lebedev staggered grid
CN105158797B (zh) * 2015-10-16 2018-02-09 成都理工大学 一种基于实际地震资料的交错网格波动方程正演的方法
CN108108331B (zh) * 2017-12-13 2018-11-02 国家深海基地管理中心 一种基于拟空间域弹性波方程的有限差分计算方法
CN109490956B (zh) * 2018-11-14 2020-12-08 深圳市勘察研究院有限公司 一种基于交错网格的声波波动方程正演模拟方法及装置
CN109490955B (zh) * 2018-11-14 2021-07-20 深圳市勘察研究院有限公司 一种基于规则网格的声波波动方程正演模拟方法及装置
CN109725346B (zh) * 2018-12-17 2021-06-18 中国石油天然气集团有限公司 一种时间-空间高阶vti矩形网格有限差分方法及装置
CN110109177B (zh) * 2019-06-05 2020-07-28 吉林大学 基于旋转时空双变网格有限差分法的地震波正演模拟方法

Also Published As

Publication number Publication date
CN112526605A (zh) 2021-03-19

Similar Documents

Publication Publication Date Title
CN111239802B (zh) 基于地震反射波形和速度谱的深度学习速度建模方法
US11609352B2 (en) Machine learning-augmented geophysical inversion
US8406081B2 (en) Seismic imaging systems and methods employing tomographic migration-velocity analysis using common angle image gathers
US8296069B2 (en) Pseudo-analytical method for the solution of wave equations
CN109711110B (zh) 任一方向入射平面波激振边坡地震动响应模拟方法
EP2409176B1 (en) Seismic imaging systems and methods employing a fast target-oriented illumination calculation
EP2497043B1 (en) Seismic imaging systems and methods employing a 3d reverse time migration with tilted transverse isotropy
AU2019340410B2 (en) Reservoir characterization utilizing inversion of resampled seismic data
CN111025387B (zh) 一种页岩储层的叠前地震多参数反演方法
CN112883564A (zh) 一种基于随机森林的水体温度预测方法及预测系统
CN110780351B (zh) 纵波和转换波叠前联合反演方法及系统
CN110007340A (zh) 基于角度域直接包络反演的盐丘速度密度估计方法
CN112526605B (zh) 一种采用地震数值模拟勘探天然气水合物的方法
US20180299573A9 (en) Method and system for efficient extrapolation of a combined source-and-receiver wavefield
CN110703318B (zh) 叠前地震数据的直接反演方法
ZHANG et al. Numerical simulation of the six-component elastic-wave field
NL2031350B1 (en) Method and system for determining geosteering irregular observation system
CN108732619B (zh) 一种海底地球物理数据采集方法
CN115407394A (zh) 一种基于波动方程的地震记录合成方法
Sun et al. Model-Data-Driven Seismic AVO Inversion Method for Viscoelastic Media based on Improved CGANs
Khoshkholgh Uncertainty Quantification in Seismic Subsurface Modelling by Informed Proposal Monte Carlo
Guo et al. Numerical Modeling of Underwater Acoustic Propagation using the Pseudo-spectral Method
Yang et al. The Fluid Identification Forward Modeling Study of Carbonate Cave Reservoir
CN115963542A (zh) 各向异性参数敏感度矩阵的确定方法、装置、设备及介质
HE et al. Finite Difference Modeling of the Acoustic Field by Sidewall Logging Devices

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant