CN110609325A - 弹性波场数值模拟方法及系统 - Google Patents
弹性波场数值模拟方法及系统 Download PDFInfo
- Publication number
- CN110609325A CN110609325A CN201810620610.XA CN201810620610A CN110609325A CN 110609325 A CN110609325 A CN 110609325A CN 201810620610 A CN201810620610 A CN 201810620610A CN 110609325 A CN110609325 A CN 110609325A
- Authority
- CN
- China
- Prior art keywords
- wave field
- equation
- boundary
- wave
- velocity model
- 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.)
- Granted
Links
- 238000004088 simulation Methods 0.000 title claims abstract description 47
- 238000000034 method Methods 0.000 title claims abstract description 26
- 238000010521 absorption reaction Methods 0.000 claims abstract description 120
- 238000004364 calculation method Methods 0.000 claims abstract description 33
- 238000005070 sampling Methods 0.000 claims abstract description 21
- 238000006073 displacement reaction Methods 0.000 claims description 41
- 230000009466 transformation Effects 0.000 claims description 6
- 238000004590 computer program Methods 0.000 claims description 3
- 230000001902 propagating effect Effects 0.000 abstract description 15
- 238000010586 diagram Methods 0.000 description 15
- 238000005516 engineering process Methods 0.000 description 7
- 238000004422 calculation algorithm Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 230000005284 excitation Effects 0.000 description 3
- 239000011159 matrix material Substances 0.000 description 2
- 230000000644 propagated 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
- 230000002378 acidificating effect Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
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
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
公开了一种弹性波场数值模拟方法及系统。该方法可以包括：步骤1：建立速度模型，确定空间和时间的采样间隔，进行有限差分网格离散化；步骤2：根据弹性波动方程计算内部波场数据；步骤3：根据角点波场数值计算公式计算角点波场数据；步骤4：根据波场吸收边界方程和波场吸收系数计算边界波场数据；步骤5：获得速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，执行步骤6，若否，重新确定波场吸收系数，重复步骤4‑5；步骤6：输出最终波场数据。本发明的模型波场值能准确反映出弹性波在二维各向同性介质中传播的真实性，达到提高弹性波数值模拟结果精度的目的。
Description
技术领域
本发明涉及弹性波波场数值模拟技术领域，更具体地，涉及一种弹性波场数值模拟方法及系统。
背景技术
弹性波波场数值模拟技术是地震勘探中的重要技术，利用弹性波波场数值模拟技术中的正演技术，可以检验经过地震资料处理、地质解释最终得到的地下物性参数以及反演的速度模型是否能真实地反映地下的储油构造，如果离开了弹性波波场数值模拟技术，我们目前的许多地震资料的处理和解释工作将寸步难行，目前广泛应用于实际工作中的弹性波波场数值模拟技术主要是弹性波波场有限差分波场数值模拟技术，在弹性波有限差分波场数值模拟中，由于计算模型的限制，人们所碰到的一个不可避免的问题是如何解决由计算网格边界所引起的人为边界反射能量，这些来自边界的人为反射扭曲了波在无限介质中传播的真实性，因此有必要提出一种波在边界上传播的算法能使往外传播的波的能量透过边界，否则，我们将不得不扩大模型的边界，大大增加计算的工作量。为了解决边界反射问题，几十年来，人们进行了大量的研究，提出了许多不同有关边界波场的计算方法，但至今，仍没有一种十分完美的解决方法，而与声波波场数值模拟相比，弹性波有限差分波场数值模拟中的边界问题更复杂、也更难解决，目前，弹性波有限差分波场数值模拟中的边界问题研究仍是该领域里的一个重要课题，同时也相继诞生了一些针对不同边界问题的算法，其中通过波场吸收消除来自边界反射的方法比较通用且效果也较明显，国内国外都有不少这方面的研究成果出现和文献发表，目前勘探中广泛使用两类吸收边界条件：海绵吸收边界条件和傍轴近似吸收边界条件。海绵吸收边界条件利用黏滞边界或靠近边界的条带范围内对入射波进行衰减；傍轴近似吸收边界条件基于单程波传播的方法原理，将不同近似条件下的单程波方程作为不同边界处的吸收边界条件。在国内，董良国在文献中(董良国.弹性波数值模拟中的吸收边界条件.石油地球物理勘探,1999,34(1):45～56)利用特征分析方法将前人提出的一维情况下无边界反射概念推广至三维各向异性介质弹性波的数值模拟中，得到TI介质中的吸收边界条件，在国外，Clayton和Engquist等人在文献中(RobertClayton,Bjorn Engquist.absorbing boundary conditions for acoustic and elasticwave equations.Bull.Seis.Soc.Am.,1977,67(6):1529～1540)利用弹性波传播方程傍轴近似理论得到了各向同性介质中弹性波数值模拟的吸收边界条件，至今，他们所提出的吸收边界条件仍是在各向同性介质中进行弹性波有限差分数值模拟用来解决边界问题最常用的方法，此后，Clayton和Engquist在文献中(Robert Clayton,Bjorn Engquist,absorbing boundary conditions for wave equations migration.Geophysics,1980,45(3):895～904)又对弹性波有限差分波场数值模拟中的吸收边界条件进行了更深入的研究。在国内，也有不少研究人员对利用有限差分法进行弹性波场数值模拟情况下的吸收边界条件做过较详细的研究，得到了许多有价值的研究成果，他们的算法对某些模型得到了一些较好的结果，本发明第一发明人在文献中(李文杰等.一种新的弹性波数值模拟吸收边界条件.石油地球物理勘探,2009,44(4):501～507)利用二维各向同性介质的弹性波传播方程傍轴近似理论得到了另一种的吸收边界条件,该方法与Clayton和Engquist等人提出的吸收边界方法不一样，经过数值模拟结果的证明，该吸收边界条件具有较好的吸收效果。遍观已有的弹性波波场数值模拟算法，每种方法都有自己各自的特点，但所有这些吸收边界条件都不能完全吸收来自边界的反射波场能量，只是部分吸收来自边界的反射能量。因此，对于弹性波数值模拟来说，确实需要发明一种可以完全吸收来自边界反射能量的吸收边界条件，这样才能提高弹性波数值模拟结果的精度，使数值模拟结果更加符合实际波场的传播规律。因此，有必要开发一种弹性波场数值模拟方法及系统。
公开于本发明背景技术部分的信息仅仅旨在加深对本发明的一般背景技术的理解，而不应当被视为承认或以任何形式暗示该信息构成已为本领域技术人员所公知的现有技术。
发明内容
本发明提出了一种弹性波场数值模拟方法及系统，其能够通过不断调整波场吸收系数，有效地消除来自有限的速度模型边界所产生的反射能量，当吸收系数调整到最佳状态时，则可以完全吸收来自模型边界的反射能量，从而使计算得到的模型波场值能准确地反映出弹性波在二维各向同性介质中传播的真实性，进而达到提高弹性波数值模拟结果精度的目的。
根据本发明的一方面，提出了一种弹性波场数值模拟方法。所述方法可以包括：步骤1：建立速度模型，确定所述速度模型的空间采样间隔和时间采样间隔，对所述速度模型进行有限差分网格离散化；步骤2：根据弹性波动方程计算所述速度模型的内部各个网格点的各个时刻的波场值，获得内部波场数据；步骤3：根据角点波场数值计算公式计算所述速度模型的四个角点所对应的各个时刻的波场值，获得角点波场数据；步骤4：根据波场吸收边界方程和波场吸收系数计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据；步骤5：根据所述边界波场数据、所述内部波场数据、所述角点波场数据，获得所述速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5；步骤6：输出所述速度模型的最终波场数据。
优选地，所述速度模型的内部各个网格点的各个时间的波场值的表达式为：
其中，
u、w分别为弹性波场所对应的水平位移和垂直位移；
α,β分别为纵、横波速度，t为时间，x和z是二维介质情况下沿X轴和Z轴的坐标。
优选地，所述波场吸收边界方程通过以下方式获得：对所述弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程；将所述傅里叶变换后的弹性波动方程进行因式分解与化简，获得化简方程；将所述化简方程转换至时间域，获得波场吸收边界方程。
优选地，傅里叶变换后的弹性波动方程为：
优选地，所述化简方程为：
优选地，通过方程(4a)-(4d)计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，
速度模型底边所对应的吸收边界方程为：
速度模型顶边所对应的吸收边界方程为：
速度模型左边所对应的吸收边界方程为：
速度模型右边所对应的吸收边界方程为：
其中，
u、w分别为弹性波场所对应的水平位移和垂直位移；
r1，r2、r3，r4为波场吸收系数，1≥r1≥0；1≥r2≥0，1≥r3≥0；1≥r4≥0，α、β分别为纵、横波速度。
优选地，通过方程(5a)-(5d)计算所述速度模型四个角点所对应的各个时刻的波场值，
右下角角点的波场计算方程为：
左上角角点的波场计算方程为：
右上角角点的波场计算方程为：
左下角角点的波场计算方程为：
其中：
u、w分别为弹性波场所对应的水平位移场和垂直位移场；
α、β分别为纵、横波速度，x和z是二维介质情况下沿X轴和Z轴的坐标；
其中：
θ为波的入射角。
根据本发明的另一方面，提出了一种弹性波场数值模拟系统，其上存储有计算机程序，其中，所述程序被处理器执行时实现以下步骤：步骤1：建立速度模型，确定所述速度模型的空间采样间隔和时间采样间隔，对所述速度模型进行有限差分网格离散化；步骤2：根据弹性波动方程计算所述速度模型的内部各个网格点的各个时刻的波场值，获得内部波场数据；步骤3：根据角点波场数值计算公式计算所述速度模型的四个角点所对应的各个时刻的波场值，获得角点波场数据；步骤4：根据波场吸收边界方程和波场吸收系数计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据；步骤5：根据所述边界波场数据、所述内部波场数据、所述角点波场数据，获得所述速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5；步骤6：输出所述速度模型的最终波场数据。
优选地，所述速度模型的内部各个网格点的各个时间的波场值的表达式为：
其中，u、w分别为弹性波场所对应的水平位移和垂直位移；
α,β分别为纵、横波速度，t为时间，x和z是二维介质情况下沿X轴和Z轴的坐标。
优选地，所述波场吸收边界方程通过以下方式获得：对所述弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程；将所述傅里叶变换后的弹性波动方程进行因式分解与化简，获得化简方程；将所述化简方程转换至时间域，获得波场吸收边界方程。
优选地，傅里叶变换后的弹性波动方程为：
优选地，所述化简方程为：
优选地，通过方程(4a)-(4d)计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，
速度模型底边所对应的吸收边界方程为：
速度模型顶边所对应的吸收边界方程为：
速度模型左边所对应的吸收边界方程为：
速度模型右边所对应的吸收边界方程为：
其中，u、w分别为弹性波场所对应的水平位移和垂直位移；
r1，r2、r3，r4为波场吸收系数，1≥r1≥0；1≥r2≥0，1≥r3≥0；1≥r4≥0，α、β分别为纵、横波速度。
优选地，通过方程(5a)-(5d)计算所述速度模型四个角点所对应的各个时刻的波场值，
右下角角点的波场计算方程为：
左上角角点的波场计算方程为：
右上角角点的波场计算方程为：
左下角角点的波场计算方程为：
其中：u、w分别为弹性波场所对应的水平位移场和垂直位移场；
α、β分别为纵、横波速度，x和z是二维介质情况下沿X轴和Z轴的坐标；
其中：
θ为波的入射角。
本发明具有其它的特性和优点，这些特性和优点从并入本文中的附图和随后的具体实施方式中将是显而易见的，或者将在并入本文中的附图和随后的具体实施方式中进行详细陈述，这些附图和具体实施方式共同用于解释本发明的特定原理。
附图说明
通过结合附图对本发明示例性实施例进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显，其中，在本发明示例性实施例中，相同的参考标号通常代表相同部件。
图1示出了根据本发明的弹性波场数值模拟方法的步骤的流程图。
图2示出了根据本发明的一个实施例的二维各向同性介质常速度模型的示意图。
图3a和图3b分别示出了根据本发明的一个实施例的刚性边界条件下传播到120毫秒时波场X分量与Z分量的示意图。
图4a和图4b分别示出了根据本发明的一个实施例的吸收边界条件下传播到120毫秒时波场X分量与Z分量的示意图。
图5示出了根据本发明的一个实施例的双层介质速度模型的示意图。
图6a和图6b分别示出了根据本发明的一个实施例的弹性波双重吸收边界条件下传播到260毫秒时水平位移与垂直位移的示意图。
图7a和图7b分别示出了根据本发明的一个实施例的Clayton吸收边界条件下传播到260毫秒时波场水平位移与垂直位移的示意图。
具体实施方式
下面将参照附图更详细地描述本发明。虽然附图中显示了本发明的优选实施例，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了使本发明更加透彻和完整，并且能够将本发明的范围完整地传达给本领域的技术人员。
图1示出了根据本发明的弹性波场数值模拟方法的步骤的流程图。
在该实施例中，根据本发明的弹性波场数值模拟方法可以包括：步骤1：建立速度模型，确定速度模型的空间采样间隔和时间采样间隔，对速度模型进行有限差分网格离散化；步骤2：根据弹性波动方程计算速度模型的内部各个网格点的各个时刻的波场值，获得内部波场数据；步骤3：根据角点波场数值计算公式计算速度模型的四个角点所对应的各个时刻的波场值，获得角点波场数据；步骤4：根据波场吸收边界方程和波场吸收系数计算速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据；步骤5：根据边界波场数据、内部波场数据、角点波场数据，获得速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5；步骤6：输出速度模型的最终波场数据。
在一个示例中，速度模型的内部各个网格点的各个时间的波场值的表达式为：
其中，u,w分别为弹性波场所对应的水平位移和垂直位移；
α,β分别为纵、横波速度，t为时间。
在一个示例中，波场吸收边界方程通过以下方式获得：对弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程；将傅里叶变换后的弹性波动方程进行因式分解与化简，获得化简方程；将化简方程转换至时间域，获得波场吸收边界方程。
在一个示例中，傅里叶变换后的弹性波动方程为：
在一个示例中，化简方程为：
在一个示例中，通过方程(4a)-(4d)计算速度模型的边界上的网格点所对应的各个时刻的波场值，
速度模型底边所对应的吸收边界方程为：
速度模型顶边所对应的吸收边界方程为：
速度模型左边所对应的吸收边界方程为：
速度模型右边所对应的吸收边界方程为：
其中，u、w分别为弹性波场所对应的水平位移和垂直位移；
r1，r2、r3，r4为波场吸收系数，1≥r1≥0；1≥r2≥0，1≥r3≥0；1≥r4≥0，α、β分别为纵、横波速度。
在一个示例中，通过方程(5a)-(5d)计算速度模型四个角点所对应的各个时刻的波场值，
右下角角点的波场计算方程为：
左上角角点的波场计算方程为：
右上角角点的波场计算方程为：
左下角角点的波场计算方程为：
其中：u、w分别为弹性波场所对应的水平位移场和垂直位移场；
α、β分别为纵、横波速度，x和z是二维介质情况下沿X轴和Z轴的坐标；
其中：
θ为波的入射角。
具体地，根据本发明的弹性波场数值模拟方法可以包括：
步骤1：建立速度模型，确定炮点(即激发源)位置，选择地震子波，确定速度模型的空间采样间隔和时间采样间隔，对速度模型进行有限差分网格离散化。
步骤2：在各向同性介质条件下，假设x和z是二维介质情况下沿X轴和Z轴的坐标，X轴的方向向右，Z轴的方向朝下，我们可以利用两个共轭的二阶偏微分方程来描述介质中P波的运动以及垂向偏振的SV波的运动，在这里不考虑沿水平方向偏振的SH波，u、w分别为水平方向与垂直方向的位移，ρ为介质的密度，t为时间，λ和μ为具体介质中的拉梅系数，那么可以得到关于各向同性非均匀介质的弹性波动方程为公式(5)：
假定密度ρ为常数，这样就可以把公式(6)看成是随空间位置变化的P波和SV波速度的函数。其中λ和μ与纵、横波速度α(x,z)和β(x,z)的关系为公式(6)：
进而将公式(6)改写为矩阵形式，即为公式(1)，计算速度模型的内部各个网格点的各个时刻的波场值为公式(1)，获得内部波场数据。
步骤3：根据角点波场数值计算公式计算速度模型的四个角点的波场值分别为公式(5a)-(5d)，获得角点波场数据，根据炮点位置坐标以及四个角点所对应的位置坐标，计算出sinθ和cosθ，由此得到J矩阵向量中的各项值。
步骤4：根据弹性波动方程与波场吸收系数获得波场吸收边界方程，根据波场吸收边界方程计算对应波场吸收系数的速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据，包括：
对弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程为公式(2)；为了能对公式(2)的左边项进行分解，考虑将公式(2)的左边项变为公式(8)：
其中：
对(8)式进行整理，得到公式(9)：
由(9)式得到(10)式：
令：则有公式(11)：
进而获得公式(12)：
其中：1≥r1≥0；1≥r2≥0。
于是，(10)式可以写成下面的推导式子：
可以将(13)式可改写成下面的形式：
I+C(kz/ω)-A(kx/ω)-B(kz/ω)＝0 (14)
其中：
这样，通过因式分解，公式(2)便可以将改写成下面的形式：
令G＝I+C(kz/ω)+A(kx/ω)+B(kz/ω)，在(15)式的两边乘上G-1，则获得化简方程为公式(3)；将化简方程转换至时间域，获得波场吸收边界方程；根据波场吸收边界方程计算对应波场吸收系数的速度模型的边界上的网格点所对应的各个时刻的波场值为公式(4a)-(4d)，获得边界波场数据。
步骤5：根据边界波场数据、内部波场数据、角点波场数据，获得速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5，其中，本领域技术人员可以根据速度模型的不同来确定预定标准。
步骤6：输出速度模型的最终波场数据。
本方法通过不断调整波场吸收系数，有效地消除来自有限的速度模型边界所产生的反射能量，当吸收系数调整到最佳状态时，则可以完全吸收来自模型边界的反射能量，从而使计算得到的模型波场值能准确地反映出弹性波在二维各向同性介质中传播的真实性，进而达到提高弹性波数值模拟结果精度的目的。
为便于理解本发明实施例的方案及其效果，以下给出两个具体应用示例。本领域技术人员应理解，该示例仅为了便于理解本发明，其任何具体细节并非意在以任何方式限制本发明。
应用示例1
图2示出了根据本发明的一个实施例的二维各向同性介质常速度模型的示意图。
步骤1：建立二维各向同性介质常速度模型，如图2所示，纵波速度为3000米/秒，横波速度为2000米/秒，激发源为纵波震源，位于速度模型的中央。选取主频为50Hz的雷克子波作为我们的地震子波，空间采样间隔为1米，时间采样间隔为0.25ms，然后根据所选的空间采样对所选的速度模型进行网格离散化。
步骤2：根据公式(1)计算速度模型的内部各个网格点的各个时刻的波场值，获得内部波场数据。
步骤3：根据角点波场数值计算公式计算速度模型的角点的波场值为公式(5a)-(5d)，获得角点波场数据。
步骤4：根据弹性波动方程与波场吸收系数获得波场吸收边界方程，根据波场吸收边界方程计算对应波场吸收系数的速度模型的边界上的网格点所对应的各个时刻的波场值为公式(4a)-(4d)，获得边界波场数据。
步骤5：根据边界波场数据、内部波场数据、角点波场数据，获得速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5。
步骤6：输出速度模型的最终波场数据。
图3a和图3b分别示出了根据本发明的一个实施例的刚性边界条件(即波场传播到边界时出现全反射的情况)下传播到120毫秒时波场X分量与Z分量的示意图。图4a和图4b分别示出了根据本发明的一个实施例的吸收边界条件下传播到120毫秒时波场X分量与Z分量的示意图。对比图3和图4，很容易看出，本发明所使用的吸收边界条件可以对传播到边界的波场起到很好的吸收作用。
应用示例2
图5示出了根据本发明的一个实施例的双层介质速度模型的示意图。
步骤1：建立二维各向同性层状介质速度模型，如图5所示，第一层介质的纵波速度为2000米/秒，横波速度为1000米/秒，第二层介质的纵波速度为4000米/秒，横波速度为2000米/秒，激发源位于模型顶部以下300米，距左右两边各300米。选取主频为50Hz的雷克子波作为我们的地震子波，空间采样间隔为1米，时间采样间隔为0.25ms，然后根据所选的空间采样对所选的速度模型进行网格离散化。
步骤2：根据公式(1)计算速度模型的内部各个网格点的各个时刻的波场值，获得内部波场数据。
步骤3：根据角点波场数值计算公式计算速度模型的角点的波场值分别为公式(5a)-(5d)，获得角点波场数据。
步骤4：根据弹性波动方程与波场吸收系数获得波场吸收边界方程，根据波场吸收边界方程计算对应波场吸收系数的速度模型的边界上的网格点所对应的各个时刻的波场值为公式(4a)-(4d)，获得边界波场数据。
步骤5：根据边界波场数据、内部波场数据、角点波场数据，获得速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5。
步骤6：输出速度模型的最终波场数据。
图6a和图6b分别示出了根据本发明的一个实施例的弹性波双重吸收边界条件下传播到260毫秒时水平位移与垂直位移的示意图。图7a和图7b分别示出了根据本发明的一个实施例的Clayton吸收边界条件下传播到260毫秒时波场水平位移与垂直位移的示意图。在上层介质，这时直达纵波已经传出了模型外，来自两层介质分界面上的反射波也快到达顶部边界；在下层介质，透射波也已传到底部，此时模型的波场分布图上只有透射波和反射波，从图上对比Clayton吸收边界条件和弹性波双重吸收边界条件两种情况下的波场分布状况，可以看到：利用Clayton吸收边界条件得到的波场分布图上出现了来自边界的反射，而利用本发明提出的弹性波双重吸收边界条件却能收到很好的吸收效果。
综上所述，本发明通过不断调整波场吸收系数，有效地消除来自有限的速度模型边界所产生的反射能量，当吸收系数调整到最佳状态时，则可以完全吸收来自模型边界的反射能量，从而使计算得到的模型波场值能准确地反映出弹性波在二维各向同性介质中传播的真实性，进而达到提高弹性波数值模拟结果精度的目的。
本领域技术人员应理解，上面对本发明的实施例的描述的目的仅为了示例性地说明本发明的实施例的有益效果，并不意在将本发明的实施例限制于所给出的任何示例。
根据本发明的另一方面，提出了一种弹性波场数值模拟系统，其上存储有计算机程序，其中，程序被处理器执行时实现以下步骤：步骤1：建立速度模型，确定速度模型的空间采样间隔和时间采样间隔，对速度模型进行有限差分网格离散化；步骤2：根据弹性波动方程计算速度模型的内部各个网格点的各个时刻的波场值，获得内部波场数据；步骤3：根据角点波场数值计算公式计算速度模型的四个角点所对应的各个时刻的波场值，获得角点波场数据；步骤4：根据波场吸收边界方程和波场吸收系数计算速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据；步骤5：根据边界波场数据、内部波场数据、角点波场数据，获得速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5；步骤6：输出速度模型的最终波场数据。
在一个示例中，速度模型的内部各个网格点的各个时间的波场值的表达式为：
其中，u、w分别为弹性波场所对应的水平位移和垂直位移；
α,β分别为纵、横波速度，t为时间。
在一个示例中，波场吸收边界方程通过以下方式获得：对弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程；将傅里叶变换后的弹性波动方程进行因式分解与化简，获得化简方程；将化简方程转换至时间域，获得波场吸收边界方程。
在一个示例中，傅里叶变换后的弹性波动方程为：
在一个示例中，化简方程为：
在一个示例中，通过方程(4a)-(4d)计算速度模型的边界上的网格点所对应的各个时刻的波场值，
速度模型底边所对应的吸收边界方程为：
速度模型顶边所对应的吸收边界方程为：
速度模型左边所对应的吸收边界方程为：
速度模型右边所对应的吸收边界方程为：
其中，u、w分别为弹性波场所对应的水平位移和垂直位移； r1，r2、r3，r4为波场吸收系数，1≥r1≥0；1≥r2≥0，1≥r3≥0；1≥r4≥0，α、β分别为纵、横波速度。
在一个示例中，通过方程(5a)-(5d)计算速度模型四个角点所对应的各个时刻的波场值，
右下角角点的波场计算方程为：
左上角角点的波场计算方程为：
右上角角点的波场计算方程为：
左下角角点的波场计算方程为：
其中：u、w分别为弹性波场所对应的水平位移场和垂直位移场；
α、β分别为纵、横波速度，x和z是二维介质情况下沿X轴和Z轴的坐标；
其中：
θ为波的入射角。
本系统通过不断调整波场吸收系数，有效地消除来自有限的速度模型边界所产生的反射能量，当吸收系数调整到最佳状态时，则可以完全吸收来自模型边界的反射能量，从而使计算得到的模型波场值能准确地反映出弹性波在二维各向同性介质中传播的真实性，进而达到提高弹性波数值模拟结果精度的目的。
以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。
Claims (10)
1.一种弹性波场数值模拟方法，包括：
步骤1：建立速度模型，确定所述速度模型的空间采样间隔和时间采样间隔，对所述速度模型进行有限差分网格离散化；
步骤2：根据弹性波动方程计算所述速度模型的内部各个网格点的各个时刻的波场值，获得内部波场数据；
步骤3：根据角点波场数值计算公式计算所述速度模型的四个角点所对应的各个时刻的波场值，获得角点波场数据；
步骤4：根据波场吸收边界方程和波场吸收系数计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据；
步骤5：根据所述边界波场数据、所述内部波场数据、所述角点波场数据，获得所述速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5；
步骤6：输出所述速度模型的最终波场数据。
2.根据权利要求1所述的弹性波场数值模拟方法，其中，所述速度模型的内部各个网格点的各个时间的波场值的表达式为：
其中，
u、w分别为弹性波场所对应的水平位移和垂直位移；
α,β分别为纵、横波速度，t为时间，x和z是二维介质情况下沿X轴和Z轴的坐标。
3.根据权利要求1所述的弹性波场数值模拟方法，其中，所述波场吸收边界方程通过以下方式获得：
对所述弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程；
将所述傅里叶变换后的弹性波动方程进行因式分解与化简，获得化简方程；
将所述化简方程转换至时间域，获得波场吸收边界方程。
4.根据权利要求3所述的弹性波场数值模拟方法，其中，傅里叶变换后的弹性波动方程为：
5.根据权利要求3所述的弹性波场数值模拟方法，其中，所述化简方程为：
6.根据权利要求3所述的弹性波场数值模拟方法，其中，通过方程(4a)-(4d)计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，
速度模型底边所对应的吸收边界方程为：
速度模型顶边所对应的吸收边界方程为：
速度模型左边所对应的吸收边界方程为：
速度模型右边所对应的吸收边界方程为：
其中，u、w分别为弹性波场所对应的水平位移和垂直位移；
r1，r2、r3，r4为波场吸收系数，1≥r1≥0；1≥r2≥0，1≥r3≥0；1≥r4≥0，α、β分别为纵、横波速度。
7.根据权利要求1所述的弹性波场数值模拟方法，其中，通过方程(5a)-(5d)计算所述速度模型四个角点所对应的各个时刻的波场值：
右下角角点的波场计算方程为：
左上角角点的波场计算方程为：
右上角角点的波场计算方程为：
左下角角点的波场计算方程为：
其中：u、w分别为弹性波场所对应的水平位移场和垂直位移场；
α、β分别为纵、横波速度，x和z是二维介质情况下沿X轴和Z轴的坐标；
其中：
θ为波的入射角。
8.一种弹性波场数值模拟系统，其上存储有计算机程序，其中，所述程序被处理器执行时实现以下步骤：
步骤1：建立速度模型，确定所述速度模型的空间采样间隔和时间采样间隔，对所述速度模型进行有限差分网格离散化；
步骤2：根据弹性波动方程计算所述速度模型的内部各个网格点的各个时刻的波场值，获得内部波场数据；
步骤3：根据角点波场数值计算公式计算所述速度模型的四个角点所对应的各个时刻的波场值，获得角点波场数据；
步骤4：根据波场吸收边界方程和波场吸收系数计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据；
步骤5：根据所述边界波场数据、所述内部波场数据、所述角点波场数据，获得所述速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5；
步骤6：输出所述速度模型的最终波场数据。
9.根据权利要求8所述的弹性波场数值模拟系统，其中，根据所述弹性波动方程获得波场吸收边界方程，根据所述波场吸收边界方程计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据包括：
对所述弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程；
将所述傅里叶变换后的弹性波动方程进行因式分解与化简，获得化简方程；
将所述化简方程转换至时间域，获得波场吸收边界方程；
根据所述波场吸收边界方程计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据。
10.根据权利要求9所述的弹性波场数值模拟系统，其中，通过方程(4a)-(4d)计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，
速度模型底边所对应的吸收边界方程为：
速度模型顶边所对应的吸收边界方程为：
速度模型左边所对应的吸收边界方程为：
速度模型右边所对应的吸收边界方程为：
其中，u、w分别为弹性波场所对应的水平位移和垂直位移；
r1，r2、r3，r4为波场吸收系数，1≥r1≥0；1≥r2≥0，1≥r3≥0；1≥r4≥0，α、β分别为纵、横波速度。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810620610.XA CN110609325B (zh) | 2018-06-15 | 2018-06-15 | 弹性波场数值模拟方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810620610.XA CN110609325B (zh) | 2018-06-15 | 2018-06-15 | 弹性波场数值模拟方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110609325A true CN110609325A (zh) | 2019-12-24 |
| CN110609325B CN110609325B (zh) | 2021-02-09 |
Family
ID=68888238
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810620610.XA Active CN110609325B (zh) | 2018-06-15 | 2018-06-15 | 弹性波场数值模拟方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110609325B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111257930A (zh) * | 2020-02-18 | 2020-06-09 | 中国石油大学(华东) | 一种黏弹各向异性双相介质区域变网格求解算子 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6687659B1 (en) * | 2000-03-24 | 2004-02-03 | Conocophillips Company | Method and apparatus for absorbing boundary conditions in numerical finite-difference acoustic applications |
| CN105807317A (zh) * | 2016-05-06 | 2016-07-27 | 中国地质大学（北京） | 基于切比雪夫伪谱法的各向异性衰减面波模拟方法 |
| WO2016130208A1 (en) * | 2015-02-13 | 2016-08-18 | Exxonmobil Upstream Research Company | Efficient and stable absorbing boundary condition in finite-difference calculations |
| CN107315192A (zh) * | 2016-04-26 | 2017-11-03 | 中国石油化工股份有限公司 | 基于二维各向同性介质的弹性波波场数值的模拟方法 |
| CN107942375A (zh) * | 2017-11-17 | 2018-04-20 | 河海大学 | 一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法 |
-
2018
- 2018-06-15 CN CN201810620610.XA patent/CN110609325B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6687659B1 (en) * | 2000-03-24 | 2004-02-03 | Conocophillips Company | Method and apparatus for absorbing boundary conditions in numerical finite-difference acoustic applications |
| WO2016130208A1 (en) * | 2015-02-13 | 2016-08-18 | Exxonmobil Upstream Research Company | Efficient and stable absorbing boundary condition in finite-difference calculations |
| CN107315192A (zh) * | 2016-04-26 | 2017-11-03 | 中国石油化工股份有限公司 | 基于二维各向同性介质的弹性波波场数值的模拟方法 |
| CN105807317A (zh) * | 2016-05-06 | 2016-07-27 | 中国地质大学（北京） | 基于切比雪夫伪谱法的各向异性衰减面波模拟方法 |
| CN107942375A (zh) * | 2017-11-17 | 2018-04-20 | 河海大学 | 一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法 |
Non-Patent Citations (3)
| Title |
|---|
| J. P. NARAYAN: "Absorbing boundary conditions in a fourth-order accurate SH-wave staggered grid finite difference algorithm", 《ACTA GEOPHYSICA》 * |
| ROBERT CLAYTON，等: "ABSORBING BOUNDARY CONDITIONS FOR ACOUSTIC AND ELASTIC WAVE EQUATIONS", 《BULLETIN OF THE SEISMOLOGICAL SOCIETY OF AMERICA》 * |
| 李文杰，等: "一种新的弹性波数值模拟吸收边界条件", 《石油地球物理勘探》 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111257930A (zh) * | 2020-02-18 | 2020-06-09 | 中国石油大学(华东) | 一种黏弹各向异性双相介质区域变网格求解算子 |
| CN111257930B (zh) * | 2020-02-18 | 2022-03-29 | 中国石油大学(华东) | 一种黏弹各向异性双相介质区域变网格求解算子 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110609325B (zh) | 2021-02-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Yao et al. | A review on reflection-waveform inversion | |
| Ren et al. | A hierarchical elastic full-waveform inversion scheme based on wavefield separation and the multistep-length approach | |
| US20170299745A1 (en) | Prestack egs migration method for seismic wave multi-component data | |
| Gazdag et al. | Migration of seismic data by phase shift plus interpolation | |
| US20210103065A1 (en) | Determining properties of a subterranean formation using an acoustic wave equation with a reflectivity parameterization | |
| AU2013205827A1 (en) | Methods and systems for imaging subterranean formations with primary and multiple reflections | |
| US20180210101A1 (en) | System and method for seismic inversion | |
| CN108051855B (zh) | 一种基于拟空间域声波方程的有限差分计算方法 | |
| CN108108331A (zh) | 一种基于拟空间域弹性波方程的有限差分计算方法 | |
| KR20190075176A (ko) | 다중반사파 없는 데이터 세트를 생성하는 다단식 전 파동장 역산 프로세스 | |
| US11105945B2 (en) | Processes and systems that attenuate source signatures and free-surface effects in recorded seismic data | |
| Keers et al. | Viscoacoustic crosswell imaging using asymptotic waveforms | |
| Pica | Fast and accurate finite-difference solutions of the 3D eikonal equation parametrized in celerity. | |
| Diouane et al. | A parallel evolution strategy for an earth imaging problem in geophysics | |
| Fu et al. | Infinite boundary element absorbing boundary for wave propagation simulations | |
| CN105447225B (zh) | 一种应用于声波有限差分数值模拟的组合吸收边界条件 | |
| CN110609325B (zh) | 弹性波场数值模拟方法及系统 | |
| Sun et al. | Multiple attenuation using λ-f domain high-order and high-resolution Radon transform based on SL0 norm | |
| US10379245B2 (en) | Method and system for efficient extrapolation of a combined source-and-receiver wavefield | |
| Hu et al. | Numerical modeling of seismic waves using frequency-adaptive meshes | |
| US9791580B2 (en) | Methods and systems to separate wavefields using pressure wavefield data | |
| CN104732093A (zh) | 一种基于弥散黏滞性波动方程的fct-fdm正演模拟方法 | |
| Bader | Seismic and Newtonian noise modeling for advanced Virgo and Einstein telescope | |
| CN107315192B (zh) | 基于二维各向同性介质的弹性波波场数值的模拟方法 | |
| Zhao et al. | Frequency-domain finite-difference elastic wave modeling in the presence of surface topography |
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 |