CN112764105B - Hti介质准纵波正演模拟方法、装置、存储介质及处理器 - Google Patents
Hti介质准纵波正演模拟方法、装置、存储介质及处理器 Download PDFInfo
- Publication number
- CN112764105B CN112764105B CN202011578765.5A CN202011578765A CN112764105B CN 112764105 B CN112764105 B CN 112764105B CN 202011578765 A CN202011578765 A CN 202011578765A CN 112764105 B CN112764105 B CN 112764105B
- Authority
- CN
- China
- Prior art keywords
- equation
- wave
- quasi
- hti medium
- longitudinal wave
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 64
- 238000004088 simulation Methods 0.000 title claims abstract description 56
- 238000003860 storage Methods 0.000 title claims abstract description 20
- 238000010521 absorption reaction Methods 0.000 claims abstract description 35
- 238000009795 derivation Methods 0.000 claims description 19
- 230000010354 integration Effects 0.000 claims description 12
- 230000001413 cellular effect Effects 0.000 claims description 8
- 230000015654 memory Effects 0.000 description 27
- 230000008569 process Effects 0.000 description 23
- 238000010586 diagram Methods 0.000 description 20
- 238000004590 computer program Methods 0.000 description 10
- 238000004364 calculation method Methods 0.000 description 8
- 230000008878 coupling Effects 0.000 description 7
- 238000010168 coupling process Methods 0.000 description 7
- 238000005859 coupling reaction Methods 0.000 description 7
- 239000011435 rock Substances 0.000 description 5
- 230000000694 effects Effects 0.000 description 4
- 230000006870 function Effects 0.000 description 4
- 239000011159 matrix material Substances 0.000 description 4
- 230000010287 polarization Effects 0.000 description 4
- 238000012512 characterization method Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 230000005540 biological transmission Effects 0.000 description 2
- 239000012530 fluid Substances 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 239000011148 porous material Substances 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 230000002123 temporal effect Effects 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- QNRATNLHPGXHMA-XZHTYLCXSA-N (r)-(6-ethoxyquinolin-4-yl)-[(2s,4s,5r)-5-ethyl-1-azabicyclo[2.2.2]octan-2-yl]methanol;hydrochloride Chemical compound Cl.C([C@H]([C@H](C1)CC)C2)CN1[C@@H]2[C@H](O)C1=CC=NC2=CC=C(OCC)C=C21 QNRATNLHPGXHMA-XZHTYLCXSA-N 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000001902 propagating effect Effects 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/63—Seismic attributes, e.g. amplitude, polarity, instant phase
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/67—Wave propagation modeling
- G01V2210/673—Finite-element; Finite-difference
-
- Y—GENERAL 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
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
- Y02E60/00—Enabling technologies; Technologies with a potential or indirect contribution to GHG emissions mitigation
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)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明实施例提供了一种HTI介质准纵波正演模拟方法、装置、存储介质及处理器，属于地球物理勘查技术领域，所述方法包括：推导HTI介质准纵波积分方程；采用分裂完全匹配层吸收边界，为所述HTI介质准纵波积分方程设置边界条件；利用高阶有限差分对设置边界条件后的HTI介质准纵波积分方程进行正演模拟。
Description
技术领域
本发明涉及地球物理勘查技术领域，尤其涉及一种HTI介质准纵波正演模拟方法、装置、存储介质及处理器。
背景技术
地震波正演模拟主要是求取地震波在已知的地下地质模型中的传播规律,包括传播时间、路径、能量等。通过正演模拟，可以正确认识地震波在复杂介质中传播的运动学和动力学特征,以及准确地分析地下地质构造所产生的反射地震波场特征。
相关技术中，常采用波动方程对地震波进行正演模拟，波动方程法是应对日益复杂的地质构造地震波传播数值模拟的最重要方法。
现实中，由于地下介质存在裂隙、裂缝、砂泥岩薄互层等，是各向异性的，采用各向同性波动方程进行正演模拟，模拟波在地下传播时会有较大的误差，因此通常采用各向异性波动方程来模拟波在地下传播过程。
另外，由于受到应力场的作用，岩石中形成择优取向排列的裂缝、裂隙和孔隙，这些裂缝、裂隙或孔隙可能充满气体或流体等充填物，地震波在这些裂隙岩石中的传播相当于在均匀弹性各向异性固体中的传播，因此可称此裂隙岩石具有等效各向异性，也可以把这种具有气体或流体等充填物的择优取向裂隙称为广泛扩容各向异性。通常采用具有水平对称轴的横向各向同性 (Horizontal Transverse Isotropy，HTI)介质来等效这种裂隙岩石进行正演模拟。
然而，相关技术中，基于HTI介质的地震波正演模拟过程尚需优化。
发明内容
本发明实施例的目的是提供一种HTI介质准纵波正演模拟方法、装置、存储介质及处理器。
为了实现上述目的，本发明第一方面提供一种HTI介质准纵波正演模拟方法，包括：
推导HTI介质准纵波积分方程；
采用分裂完全匹配层吸收边界，为HTI介质准纵波积分方程设置边界条件；
利用高阶有限差分对设置边界条件后的HTI介质准纵波积分方程进行正演模拟。
在本发明实施例中，推导HTI介质准纵波积分方程，包括：
基于柯西方程、几何方程和纳维尔方程得到HTI介质弹性波方程，令弹性波方程体力项为零，并将平面波解带入弹性波方程得到开尔文-克里斯托弗方程；
对所述开尔文-克里斯托弗方程进行求解，并利用声学假设令横波速度为零，得到HTI介质准纵波近似方程；
对HTI介质准纵波近似方程进行双重时间积分，以推导出HTI介质准纵波积分方程。
在本发明实施例中，基于柯西方程、几何方程和纳维尔方程得到开尔文 -克里斯托弗方程，包括：
基于柯西方程、几何方程和纳维尔方程得到二阶弹性波波动方程；
利用所述二阶弹性波波动方程与平面波的解得到所述开尔文-克里斯托弗方程。
在本发明实施例中，HTI介质准纵波近似方程被定义为：
其中，ε、δ表示各向异性参数、vp0表示纵波速度，P表示压强，x、y、 z为三个相互垂直的空间方向，t表示时间，
表示对x方向取空间偏导数，
表示对y方向取空间偏导数，
表示对z方向取空间偏导数，
表示对时间的偏导数。
在本发明实施例中，HTI介质准纵波积分方程被定义为：
其中，ε、δ表示各向异性参数、vp0表示纵波速度，P表示压强，x、y、 z为三个相互垂直的空间方向，t表示时间，
表示对x方向取空间偏导数，
表示对y方向取空间偏导数，
表示对z方向取空间偏导数，
表示对时间的偏导数。
在本发明实施例中，采用分裂完全匹配层吸收边界，为HTI介质准纵波积分方程设置边界条件，包括：
确定分裂完全匹配层吸收边界的复拉伸变量；
基于复拉伸变量进行一阶偏导，得到复数坐标系下的一阶空间导数；
基于一阶空间导数得到二阶空间偏导和四阶空间偏导；
基于二阶空间偏导和四阶空间偏导，将HTI介质准纵波积分方程分解为多个波场，完成边界条件设置。
在本发明实施例中，利用高阶有限差分对设置边界条件后的HTI介质准纵波积分方程进行正演模拟，包括：
利用时间二阶差分格式、空间十阶差分格式对HTI介质准纵波积分方程进行有限差分处理，得到离散方程；
基于离散方程进行HTI介质准纵波正演模拟。
本发明第二方面提供一种HTI介质准纵波正演模拟装置，包括：
推导模块，用于确定HTI介质准纵波积分方程；
设置模块，用于采用分裂完全匹配层吸收边界，为HTI介质准纵波积分方程设置边界条件；
正演模拟模块，用于利用高阶有限差分对设置边界条件后的HTI介质准纵波积分方程进行正演模拟。
本发明第三方面提供一种机器可读存储介质，机器可读存储介质上存储有指令，指令用于使得机器执行本申请上述任一项的HTI介质准纵波正演模拟方法。
本发明第四方面提供一种处理器，程序被处理器运行时用于执行本申请上述任一项的HTI介质准纵波正演模拟方法。
通过上述技术方案，由于qP波积分方程为时间的二阶导数，只有一组共轭解，而qP波近似方程为时间的四阶导数，其有两组解，会存在慢纵波干扰，因此，相比于qP波近似方程，采用qP波积分方程进行正演模拟，消除了慢纵波干扰，同时由于采用qP波积分方程进行正演模拟过程中，无需引入中间变量，因此占用内存低，计算效率高。
本发明实施例的其它特征和优点将在随后的具体实施方式部分予以详细说明。
附图说明
附图是用来提供对本发明实施例的进一步理解，并且构成说明书的一部分，与下面的具体实施方式一起用于解释本发明实施例，但并不构成对本发明实施例的限制。在附图中：
图1是本发明实施例HTI介质准纵波正演模拟方法的流程示意图；
图2是本发明应用实施例数值模拟过程示意图；
图3a是本发明应用实施例各向同性介质波动方程XOZ平面的波场快照示意图；
图3b是本发明应用实施例HTI介质近似方程XOZ平面的波场快照示意图；
图3c是本发明应用实施例HTI介质积分方程XOZ平面的波场快照示意图；
图4是本发明应用实施例添加分裂完全匹配层吸收边界下近似方程 XOZ平面的波场快照示意图；
图5是本发明应用实施例XOZ平面的Marmousi纵波速度模型示意图；
图6是本发明应用实施例XOZ平面的Marmousi纵波速度模型的地震记录示意图；
图7是本发明实施例HTI介质地震波正演模拟装置的结构框图；
图8是本发明实施例计算机设备的内部结构图。
具体实施方式
以下结合附图对本发明实施例的具体实施方式进行详细说明。应当理解的是，此处所描述的具体实施方式仅用于说明和解释本发明实施例，并不用于限制本发明实施例。
相关技术中，对于各向异性介质波动方程的正演模拟，常用弹性波动方程进行模拟。然而采用弹性波动方程进行模拟虽然能得到丰富的地震信息，但是计算量大。
因此，Alkhalifah(1998)提出，通过将沿对称轴的剪切速度设置为零来引入横向各向同性(TI)介质的伪声学近似，也就是准纵波(也称为qP波) 近似方程模拟。相对于弹性波方程的正演模拟，qP波方程的正演模拟具有占用资源小、计算速度快、频散小于弹性波的优点。
但由于qP波近似方程为时间的四阶导数，其有两组解，因此正演模拟过程会存在慢纵波干扰，同时，由于该方法在正演模拟过程中需要引入中间变量，因此需要较多的内存、计算效率较低。
基于此，本发明实施例通过获取HTI介质准纵波积分方程；采用分裂完全匹配层吸收边界，为HTI介质准纵波积分方程设置边界条件；利用高阶有限差分数值对设置边界条件后的HTI介质准纵波积分方程进行正演模拟。由于在正演模拟过程中，是采用qP波积分方程，而非采用qP波近似方程，qP 波积分方程为时间的二阶导数，只有一组共轭解，因此在正演模拟过程中，消除了慢纵波干扰，同时，由于正演模拟过程中无需引入中间变量，因此，占用内存低，计算效率高。
本发明实施例提供了一种准纵波正演模拟方法，如图1所示，该方法包括：
步骤101：推导HTI介质准纵波积分方程；
步骤102：采用分裂完全匹配层吸收边界，为HTI介质准纵波积分方程设置边界条件；
步骤103：利用高阶有限差分对设置边界条件后的HTI介质准纵波积分方程进行正演模拟。
实际应用时，HTI介质表示各向同性介质中分布着一组平行的定向排列的垂直裂隙所构成的各向异性介质模型。
实际应用时，准纵波是一种标量波，不具有方向性。
在一实施例中，推导HTI介质准纵波积分方程，包括：
基于柯西方程、几何方程和纳维尔方程得到HTI介质弹性波方程，令弹性波方程体力项为零，并将平面波解带入弹性波方程得到开尔文-克里斯托弗方程；
对所述开尔文-克里斯托弗方程进行求解，并利用声学假设令横波速度为零，得到HTI介质准纵波近似方程；
对所述HTI介质准纵波近似方程进行双重时间积分，以推导出所述HTI 介质准纵波积分方程。
对HTI介质准纵波近似方程进行双重时间积分，以推导出HTI介质准纵波积分方程。
具体地，基于柯西方程、几何方程和纳维尔方程得到Kelvin-Christoffel 方程，包括：
基于柯西方程、几何方程和纳维尔方程得到二阶弹性波波动方程；
利用所述二阶弹性波波动方程与平面波的解得到所述开尔文-克里斯托弗方程。
实际应用时，二阶弹性波波动方程(体力项为零)被定义为：
其中，U表示波场；ρ为密度；C为刚度张量矩阵；L为偏导算子矩阵；LT为偏导算子矩阵L的转置矩阵。
弹性波方程(1)的平面波解为：
U＝Pexp[ik(n*x-vt)] (3)
其中U＝(ux,uy,uz)T为位移矢量，x＝(x,y,z)T为位置矢量，n＝(nx,ny,nz)T为波的传播方向，v为平面波传播速度，k表示波数，P(px,py,pz)T表示波的偏振方向。
实际应用时，Kelvin-Christoffel方程被定义为：
其中，Γ为多项式，v是速度；ρ为密度；px、py、pz分别为x、y、z方向的偏振方向。
实际应用时，由于偏振方向不为零，左侧矩阵的模为零，将HTI介质刚度张量参数代入Kelvin-Christoffel方程中，对Kelvin-Christoffel方程进行求解，获得HTI介质标量波耦合方程。
实际应用时，通过对HTI介质标量波耦合方程进行解耦，并采用Thomesn 表征参数进行表征，将横波速度设置为零，可以获得HTI介质qP波近似方程。
实际应用时，HTI介质准纵波近似方程被定义为：
其中，ε、δ表示各向异性参数、vp0表示纵波速度，P表示压强，x、y、 z为三个相互垂直的空间方向，t表示时间，
表示对x方向取空间偏导数，
表示对y方向取空间偏导数，
表示对z方向取空间偏导数，
表示对时间的偏导数。
实际应用时，由于HTI介质qP波近似方程为时间的四阶导数，具有两组解，一组为正常波动解，另一组为慢纵波干扰解。而慢纵波会对正演模拟过程造成干扰，因此，本发明实施例采用对近似方程做双重时间积分得到 HTI介质qP波积分方程，qP波积分方程为时间的二阶导数，只有一组共轭解，可直接消除HTI介质中慢纵波干扰。
实际应用时，HTI介质准纵波积分方程被定义为：
其中，ε、δ表示各向异性参数、vp0表示纵波速度，P表示压强，x、y、 z为三个相互垂直的空间方向，t表示时间，
表示对x方向取空间偏导数，
表示对y方向取空间偏导数，
表示对z方向取空间偏导数，
表示对时间的偏导数。
实际应用时，由于正演模拟过程中是通过有限的区域模拟无限区域的波的传播过程，因此，在正演模拟过程中，需要添加人工边界。
实际应用时，可以采用分裂完全匹配层吸收边界消除人工边界反射的影响。分裂完全匹配层吸收边界能在模型区域外人工增加了一层吸收介质，当波传播到该介质中时，波的能量将被介质以指数衰减的速度迅速吸收，直至消失，而不会反弹回模型区域，无论什么频率的波以任意角度从模型区域进入吸收介质层时，都不产生反射，因此具有良好的吸收效果。
在一实施例中，采用分裂完全匹配层吸收边界，为HTI介质准纵波积分方程设置边界条件，包括：
确定分裂完全匹配层吸收边界的复拉伸变量；
基于复拉伸变量进行一阶偏导，得到复数坐标系下的一阶空间导数；
基于一阶空间导数得到二阶空间偏导和四阶空间偏导；
基于二阶空间偏导和四阶空间偏导，将HTI介质准纵波积分方程分解为多个波场，完成边界条件设置。
实际应用时，以x方向为例，分裂完全匹配层吸收边界的复拉伸变量被定义为：
其中，sx(x)表示复频域拉伸变量，
dx(x)表示衰减因子，在 PML吸收区域dx(x)>0，数值表示衰减的强弱；在地震正演模拟区域 dx(x)＝0。
复数坐标系下变量的一阶导数与原坐标系下变量的一阶偏导数的对应关系：
实际应用时，二阶空间偏导被定义为：
实际应用时，四阶空间偏导被定义为：
实际应用时，采用如下公式将积分方程P分解为5个波场。
P＝P1+P2+P3+P4+P5 (11)
在一实施例中，利用高阶有限差分对设置边界条件后的HTI介质准纵波积分方程进行正演模拟，包括：
利用时间二阶差分格式、空间十阶差分格式对HTI介质准纵波积分方程进行有限差分处理，得到离散方程；
基于离散方程进行HTI介质准纵波正演模拟。
实际应用时，时间二阶差分格式被定义为：
其中，Δt表时间采样间隔，u表示波场，t表示时间。
实际应用时，空间高阶差分格式被定义为：
其中，Δx、Δy、Δz分别表示x、y、z方向空间步长，w表示常数。
获得离散方程后，可利用离散方程对qP波在HTI介质中的传播过程进行模拟，完成qP波的正演模拟过程，为后续的波形反演和偏移成像提供基础。
本发明实施例的方案，由于qP波积分方程为时间的二阶导数，只有一组共轭解，而qP波近似方程为时间的四阶导数，其有两组解，会存在慢纵波干扰，因此，相比于qP波近似方程，采用qP波积分方程进行正演模拟，消除了慢纵波干扰，同时由于采用qP波积分方程进行正演模拟过程中，无需引入中间变量，因此占用内存低，计算效率高。
下面结合应用实施例对本发明再作进一步详细的描述。
本应用实施例提供了一种基于HTI介质qP波积分方程有限差分方法，参见图2，方法包括如下步骤：
步骤201：建立HTI介质qP波积分方程；
步骤202：添加分裂完全匹配层吸收边界，消除人工边界反射的影响；
步骤203：高阶有限差分数值模拟，展现HTI介质qP波积分方程的正确性与优点。
具体地，步骤201的实现过程如下：
通过柯西方程、几何方程以及纳维尔方程推导出二阶弹性波波动方程；二阶弹性波波动方程可以采用上述公式(1)表示；
将上述公式(3)平面波解带入上述弹性波公式(1)可得Kelvin-Christoffel 方程，Kelvin-Christoffel方程可以采用上述公式(4)表示；
由于偏振方向不为零，所以上述公式(4)中前项的模为零。带入HTI 介质刚度张量参数并求解Kelvin-Christoffel方程可得到HTI介质耦合波动方程。对HTI介质耦合波动方程进行解耦。解耦后的方程不易直接看出各参数的物理意义，因此利用Thomesn表征参数表征方程，并令横波速度为零可得到HTI介质qP波近似方程，HTI介质qP波近似方程可以采用上述公式(5) 表示；
得到的HTI介质qP波近似方程为时间的四阶导数，具有两组解，一组为正常波动解，另一组为慢纵波干扰解。在做数值模拟时，利用qP波近似方程进行正演模拟需要引入中间变量，会增大计算机内存需求、降低计算速度。因此对近似方程进行双重时间积分，得到积分方程，HTI介质qP波积分方程可以采用上述公式(6)表示；
积分方程为时间的二阶导数，只有一组共轭解，可直接消除HTI介质中慢纵波干扰。同时，在做数值模拟时不用引入中间变量，占用内存低，计算效率高。
具体地，步骤202的实现过程如下：
相较于其它吸收边界，完全匹配层吸收边界具有良好的吸收效果。因此在做数值模拟时，使用分裂完全匹配吸收边界吸收边界反射。
以x方向为例，完全匹配层吸收边界的复拉伸变量可以采用上述公式(7) 表示；
对公式(7)求一阶偏导可得到一阶空间导数，一阶空间导数可以采用上述公式(8)表示；
根据得到的复数坐标系下的一阶空间导数，对应求出二阶空间偏导和四阶空间偏导，二阶空间偏导可以采用上述公式(9)表示，四阶空间偏导可以采用上述公式(10)表示，而后为方便计算，将积分方程P分解为5个波场，具体地，按照上述公式(11)进行分解。
具体地，步骤203的实现过程如下：
为实现正演模拟，对连续方程进行时间二阶、空间十阶差分。时间二阶差分格式可以采用上述公式(12)表示，空间十阶差分格式可以采用上述公式(13)表示。
本应用实施例将几何方程、本构方程带入到运动微分方程中，并去掉体力项，得到各向异性介质弹性波的波动方程；而后通过求解Kelvin-Christoffel 方程将HTI弹性波方程转变为标量耦合方程，通过解耦、Thomsen表征参数并假设横波速度为零，得到HTI介质qP波近似方程，对近似方程积分得到 qP波积分方程。并采用常规网格高阶有限差分对HTI介质qP波积分方程进行数值模拟，可在正演模拟过程中，消除各向异性介质中慢纵波干扰和提高计算效率。
同时，本应用实施例还采用均匀模型与Marmous纵波速度模型并添加分裂完全匹配层吸收边界来进行数值模拟。通过数值模拟对比HTI介质qP 波近似方程与积分方程的运算速度、内存占用，以及分裂完全匹配层吸收边界的吸收效果，并利用marmousi纵波速度模型来验证HTI介质qP波积分方程对复杂模型的适应性。本次正演模拟过程中的主频为20Hz；时间采样间隔为1ms，总采样点为1000，均匀模型震源位于(1000m，1000m)处。
模拟结果参见图3-图6。图3a、图3b和图3c分别为各向同性介质波动方程波场快照、HTI介质近似方程波场快照和HTI介质积分方程的波场快照示意图；从图中可看出，各向同性介质声波方程波场快照为圆形，HTI介质声波方程正演模拟波场为椭圆形，符合各向异性介质波传播特征。HTI介质近似方程中除含有正常波场外还含有菱形慢纵波干扰，与其方程有两组共轭解相对应。HTI介质积分方程中只含有规则的椭圆波场，不含有慢纵波干扰，所以积分方程在做数值模拟时不含有慢纵波干扰。图4为添加分裂完全匹配层吸收边界后的HTI介质qP波波场快照，由图可以看出，波场快照内没有边界反射，所添加的分裂完全匹配层吸收边界具有良好的吸收效果。图5为 marmous纵波速度模型，图6为积分方程marmous纵波速度模型的地震记录，其震源位于(1500m,0m)处，由图5和图6可以看出，得到的地震记录都没有明显的虚假反射波，证明了本实施例方法对于复杂模型的数值模拟稳定性良好。
另外，参见表一，通过表一可以看出积分方程相较于近似方程需要内存更少、计算效率更高。
表一
为了实现本发明实施例的方法，本发明实施例还提供了一种HTI介质准纵波正演模拟装置，设置在电子设备上，如图7所示，HTI介质准纵波正演模拟装置700包括：推导模块701、设置模块702和正演模拟模块703；其中，
推导模块701，用于推导HTI介质准纵波积分方程；
设置模块702，用于采用分裂完全匹配层吸收边界，为所述HTI介质准纵波积分方程设置边界条件；
正演模拟模块703，用于利用高阶有限差分对设置边界条件后的HTI介质准纵波积分方程进行正演模拟。
在一实施例中，推导模块701，还用于：
基于柯西方程、几何方程和纳维尔方程得到Kelvin-Christoffel方程；
将Kelvin-Christoffel方程的体力项设置为零，以得到HTI介质标量波耦合方程；
对HTI介质标量波耦合方程进行解耦，以得到HTI介质准纵波近似方程；
对HTI介质准纵波近似方程进行双重时间积分，以推导出HTI介质准纵波积分方程。
在一实施例中，推导模块701，还用于：
基于柯西方程、几何方程和纳维尔方程得到二阶弹性波波动方程；
利用二阶弹性波波动方程得到Kelvin-Christoffel方程。
在一实施例中，推导模块701，还用于：
HTI介质准纵波近似方程被定义为：
其中，ε、δ表示各向异性参数、vp0表示纵波速度，P表示压强，x、y、 z为三个相互垂直的空间方向，t表示时间，
表示对x方向取空间偏导数，
表示对y方向取空间偏导数，
表示对z方向取空间偏导数，
表示对时间的偏导数。
在一实施例中，推导模块701，还用于：
HTI介质准纵波积分方程被定义为：
其中，ε、δ表示各向异性参数、vp0表示纵波速度，P表示压强，x、y、z为三个相互垂直的空间方向，t表示时间，
表示对x方向取空间偏导数，
表示对y方向取空间偏导数，
表示对z方向取空间偏导数，
表示对时间的偏导数。
在一实施例中，设置模块702，还用于：
确定分裂完全匹配层吸收边界的复拉伸变量；
基于复拉伸变量进行一阶偏导，得到复数坐标系下的一阶空间导数；
基于一阶空间导数得到二阶空间偏导和四阶空间偏导；
基于二阶空间偏导和四阶空间偏导，将HTI介质准纵波积分方程分解为多个波场，完成边界条件设置。
在一实施例中，正演模拟模块703，还用于：
利用时间二阶差分格式、空间十阶差分格式对HTI介质准纵波积分方程进行有限差分处理，得到离散方程；
基于离散方程进行HTI介质准纵波正演模拟。
实际应用时，推导模块701、设置模块702和正演模拟模块703可由HTI 介质准纵波正演模拟装置中的处理器实现。
需要说明的是：上述实施例提供的HTI介质准纵波正演模拟装置在对准纵波进行正演模拟时，仅以上述各程序模块的划分进行举例说明，实际应用时，可以根据需要而将上述处理分配由不同的程序模块完成，即将装置的内部结构划分成不同的程序模块，以完成以上描述的全部或者部分处理。另外，上述实施例提供的HTI介质准纵波正演模拟装置与HTI介质准纵波正演模拟方法实施例属于同一构思，其具体实现过程详见方法实施例，这里不再赘述。
本发明实施例提供了一种存储介质，其上存储有程序，该程序被处理器执行时实现上述HTI介质准纵波正演模拟方法。
本发明实施例提供了一种处理器，处理器用于运行程序，其中，程序运行时执行上述HTI介质准纵波正演模拟方法。
在一个实施例中，提供了一种计算机设备，该计算机设备可以是终端，其内部结构图可以如图8所示。该计算机设备包括通过系统总线连接的处理器A01、网络接口A02、显示屏A04、输入装置A05和存储器(图中未示出)。其中，该计算机设备的处理器A01用于提供计算和控制能力。该计算机设备的存储器包括内存储器A03和非易失性存储介质A06。该非易失性存储介质 A06存储有操作系统B01和计算机程序B02。该内存储器A03为非易失性存储介质A06中的操作系统B01和计算机程序B02的运行提供环境。该计算机设备的网络接口A02用于与外部的终端通过网络连接通信。该计算机程序被处理器A01执行时以实现一种HTI介质准纵波正演模拟方法。该计算机设备的显示屏A04可以是液晶显示屏或者电子墨水显示屏，该计算机设备的输入装置A05可以是显示屏上覆盖的触摸层，也可以是计算机设备外壳上设置的按键、轨迹球或触控板，还可以是外接的键盘、触控板或鼠标等。
本领域技术人员可以理解，图8中示出的结构，仅仅是与本申请方案相关的部分结构的框图，并不构成对本申请方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
本发明实施例提供了一种设备，设备包括处理器、存储器及存储在存储器上并可在处理器上运行的程序，处理器执行程序时实现上述HTI介质准纵波正演模拟方法：
本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/ 或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
在一个典型的配置中，计算设备包括一个或多个处理器(CPU)、输入/输出接口、网络接口和内存。
存储器可能包括计算机可读介质中的非永久性存储器，随机存取存储器 (RAM)和/或非易失性内存等形式，如只读存储器(ROM)或闪存(flashRAM)。存储器是计算机可读介质的示例。
计算机可读介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机的存储介质的例子包括，但不限于相变内存 (PRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器(DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体或其他内存技术、只读光盘只读存储器 (CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带，磁带磁磁盘存储或其他磁性存储设备或任何其他非传输介质，可用于存储可以被计算设备访问的信息。按照本文中的界定，计算机可读介质不包括暂存电脑可读媒体(transitorymedia)，如调制的数据信号和载波。
还需要说明的是，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、商品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、商品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括要素的过程、方法、商品或者设备中还存在另外的相同要素。
以上仅为本申请的实施例而已，并不用于限制本申请。对于本领域技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本申请的权利要求范围之内。
Claims (9)
1.一种HTI介质准纵波正演模拟方法，其特征在于，包括：
推导HTI介质准纵波积分方程；
采用分裂完全匹配层吸收边界，为所述HTI介质准纵波积分方程设置边界条件；
利用高阶有限差分对设置边界条件后的HTI介质准纵波积分方程进行正演模拟；
所述推导HTI介质准纵波积分方程，包括：
基于柯西方程、几何方程和纳维尔方程得到HTI介质弹性波方程，令弹性波方程体力项为零，并将平面波解带入弹性波方程得到开尔文-克里斯托弗方程；
对所述开尔文-克里斯托弗方程进行求解，并利用声学假设令横波速度为零，得到HTI介质准纵波近似方程；
对所述HTI介质准纵波近似方程进行双重时间积分，以推导出所述HTI介质准纵波积分方程。
2.根据权利要求1所述的HTI介质准纵波正演模拟方法，其特征在于，基于柯西方程、几何方程、纳维尔方程与平面波的解得到开尔文-克里斯托弗方程，包括：
基于柯西方程、几何方程和纳维尔方程得到二阶弹性波方程；
利用所述二阶弹性波方程与平面波的解得到所述开尔文-克里斯托弗方程。
3.根据权利要求1所述的HTI介质准纵波正演模拟方法，其特征在于，所述HTI介质准纵波近似方程被定义为：
其中，ε、δ表示各向异性参数、vp0表示纵波速度，P表示压强，x、y、z为三个相互垂直的空间方向，t表示时间，
表示对x方向取空间偏导数，
表示对y方向取空间偏导数，
表示对z方向取空间偏导数，
表示对时间的偏导数。
4.根据权利要求1所述的HTI介质准纵波正演模拟方法，其特征在于，所述HTI介质准纵波积分方程被定义为：
其中，ε、δ表示各向异性参数、vp0表示纵波速度，P表示压强，x、y、z为三个相互垂直的空间方向，t表示时间，
表示对x方向取空间偏导数，
表示对y方向取空间偏导数，
表示对z方向取空间偏导数，
表示对时间的偏导数。
5.根据权利要求1所述的HTI介质准纵波正演模拟方法，其特征在于，所述采用分裂完全匹配层吸收边界，为所述HTI介质准纵波积分方程设置边界条件，包括：
确定分裂完全匹配层吸收边界的复拉伸变量；
基于所述复拉伸变量进行一阶偏导，得到复数坐标系下的一阶空间导数；
基于所述一阶空间导数得到二阶空间偏导和四阶空间偏导；
基于所述二阶空间偏导和四阶空间偏导，将所述HTI介质准纵波积分方程分解为多个波场，完成边界条件设置。
6.根据权利要求1所述的HTI介质准纵波正演模拟方法，其特征在于，所述利用高阶有限差分对设置边界条件后的HTI介质准纵波积分方程进行正演模拟，包括：
利用时间二阶差分格式、空间十阶差分格式对所述HTI介质准纵波积分方程进行有限差分处理，得到离散方程；
基于所述离散方程进行HTI介质准纵波正演模拟。
7.一种HTI介质准纵波正演模拟装置，其特征在于，包括：
推导模块，用于推导HTI介质准纵波积分方程；
设置模块，用于采用分裂完全匹配层吸收边界，为所述HTI介质准纵波积分方程设置边界条件；
正演模拟模块，用于利用高阶有限差分对设置边界条件后的HTI介质准纵波积分方程进行正演模拟；
所述推导模块包括：
第一推导子模块，用于基于柯西方程、几何方程和纳维尔方程得到HTI介质弹性波方程，令弹性波方程体力项为零，并将平面波解带入弹性波方程得到开尔文-克里斯托弗方程；
第二推导子模块，用于对所述开尔文-克里斯托弗方程进行求解，并利用声学假设令横波速度为零，得到HTI介质准纵波近似方程；
第三推导子模块，用于对所述HTI介质准纵波近似方程进行双重时间积分，以推导出所述HTI介质准纵波积分方程。
8.一种存储介质，所述存储介质上存储有指令，其特征在于，所述指令用于使得机器执行根据权利要求1至6任一项所述的HTI介质准纵波正演模拟方法。
9.一种处理器，其特征在于，用于运行程序，其中，所述程序被所述处理器运行时用于执行根据权利要求1至6任意一项所述的HTI介质准纵波正演模拟方法。
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011109524 | 2020-10-16 | ||
| CN2020111095246 | 2020-10-16 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112764105A CN112764105A (zh) | 2021-05-07 |
| CN112764105B true CN112764105B (zh) | 2022-07-12 |
Family
ID=75696169
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011578765.5A Active CN112764105B (zh) | 2020-10-16 | 2021-02-22 | Hti介质准纵波正演模拟方法、装置、存储介质及处理器 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112764105B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115685337A (zh) * | 2022-10-24 | 2023-02-03 | 中国石油大学(华东) | 一种各向异性弹性波解耦方法、装置及计算机设备 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109946742A (zh) * | 2019-03-29 | 2019-06-28 | 中国石油大学(华东) | 一种TTI介质中纯qP波地震数据模拟方法 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8332156B2 (en) * | 2009-07-10 | 2012-12-11 | Chevron U.S.A. Inc. | Method for propagating pseudo acoustic quasi-P waves in anisotropic media |
| WO2011141440A1 (en) * | 2010-05-12 | 2011-11-17 | Shell Internationale Research Maatschappij B.V. | Seismic p-wave modelling in an inhomogeneous transversely isotropic medium with a tilted symmetry axis |
| CN105467443B (zh) * | 2015-12-09 | 2017-09-19 | 中国科学院地质与地球物理研究所 | 一种三维各向异性弹性波数值模拟方法及系统 |
| CN110261896B (zh) * | 2019-04-26 | 2021-07-20 | 中国石油化工股份有限公司 | 稳定的各向异性ti介质正演模拟方法 |
-
2021
- 2021-02-22 CN CN202011578765.5A patent/CN112764105B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109946742A (zh) * | 2019-03-29 | 2019-06-28 | 中国石油大学(华东) | 一种TTI介质中纯qP波地震数据模拟方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112764105A (zh) | 2021-05-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Yang et al. | A nearly analytic discrete method for acoustic and elastic wave equations in anisotropic media | |
| Roten et al. | High-frequency nonlinear earthquake simulations on petascale heterogeneous supercomputers | |
| de la Puente et al. | Discontinuous Galerkin methods for wave propagation in poroelastic media | |
| De Basabe et al. | Grid dispersion and stability criteria of some common finite-element methods for acoustic and elastic wave equations | |
| Maeda et al. | FDM simulation of seismic waves, ocean acoustic waves, and tsunamis based on tsunami-coupled equations of motion | |
| Yang et al. | An optimal nearly analytic discrete method for 2D acoustic and elastic wave equations | |
| CN108983285B (zh) | 一种基于矩张量的多种震源波场模拟方法及装置 | |
| Poul et al. | Efficient time-domain deconvolution of seismic ground motions using the equivalent-linear method for soil-structure interaction analyses | |
| Mercerat et al. | A nodal high-order discontinuous Galerkin method for elastic wave propagation in arbitrary heterogeneous media | |
| De Martin | Verification of a spectral-element method code for the Southern California Earthquake Center LOH. 3 viscoelastic case | |
| Yang et al. | Optimally accurate nearly analytic discrete scheme for wave-field simulation in 3D anisotropic media | |
| Duru et al. | Stable and high-order accurate boundary treatments for the elastic wave equation on second-order form | |
| Vamaraju et al. | Enriched Galerkin finite element approximation for elastic wave propagation in fractured media | |
| CN109143339A (zh) | 弹性逆时偏移成像方法及装置 | |
| Belonosov et al. | 3D numerical simulation of elastic waves with a frequency-domain iterative solver | |
| Alkhimenkov et al. | Resolving wave propagation in anisotropic poroelastic media using graphical processing units (GPUs) | |
| Zhang et al. | A procedure for 3D seismic simulation from rupture to structures by coupling SEM and FEM | |
| Martire et al. | SPECFEM2D-DG, an open-source software modelling mechanical waves in coupled solid–fluid systems: the linearized Navier–Stokes approach | |
| CN112764105B (zh) | Hti介质准纵波正演模拟方法、装置、存储介质及处理器 | |
| Meng et al. | Numerical dispersion analysis of discontinuous Galerkin method with different basis functions for acoustic and elastic wave equations | |
| Yang et al. | A strong stability-preserving predictor-corrector method for the simulation of elastic wave propagation in anisotropic media | |
| Huang et al. | A novel hybrid method based on discontinuous Galerkin method and staggered‐grid method for scalar wavefield modelling with rough topography | |
| Beresnev | Simulation of near-fault high-frequency ground motions from the representation theorem | |
| Zeng et al. | A multidomain PSTD method for 3D elastic wave equations | |
| CN109725346B (zh) | 一种时间-空间高阶vti矩形网格有限差分方法及装置 |
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 |