CN110609325B - 弹性波场数值模拟方法及系统 - Google Patents

弹性波场数值模拟方法及系统 Download PDF

Info

Publication number
CN110609325B
CN110609325B CN201810620610.XA CN201810620610A CN110609325B CN 110609325 B CN110609325 B CN 110609325B CN 201810620610 A CN201810620610 A CN 201810620610A CN 110609325 B CN110609325 B CN 110609325B
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.)
Active
Application number
CN201810620610.XA
Other languages
English (en)
Other versions
CN110609325A (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.)
China Petroleum and Chemical Corp
Sinopec Exploration and Production Research Institute
Original Assignee
China Petroleum and Chemical Corp
Sinopec Exploration and Production Research Institute
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 China Petroleum and Chemical Corp, Sinopec Exploration and Production Research Institute filed Critical China Petroleum and Chemical Corp
Priority to CN201810620610.XA priority Critical patent/CN110609325B/zh
Publication of CN110609325A publication Critical patent/CN110609325A/zh
Application granted granted Critical
Publication of CN110609325B publication Critical patent/CN110609325B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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

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：输出所述速度模型的最终波场数据。
优选地，所述速度模型的内部各个网格点的各个时间的波场值的表达式为：
Figure BDA0001697936770000041
其中，
Figure BDA0001697936770000042
u、w分别为弹性波场所对应的水平位移和垂直位移；
Figure BDA0001697936770000043
α,β分别为纵、横波速度，t为时间，x和z是二维介质情况下沿X轴和Z轴的坐标。
优选地，所述波场吸收边界方程通过以下方式获得：对所述弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程；将所述傅里叶变换后的弹性波动方程进行因式分解与化简，获得化简方程；将所述化简方程转换至时间域，获得波场吸收边界方程。
优选地，傅里叶变换后的弹性波动方程为：
Figure BDA0001697936770000044
优选地，所述化简方程为：
Figure BDA0001697936770000045
优选地，通过方程(4a)-(4d)计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，
速度模型底边所对应的吸收边界方程为：
Figure BDA0001697936770000046
速度模型顶边所对应的吸收边界方程为：
Figure BDA0001697936770000047
速度模型左边所对应的吸收边界方程为：
Figure BDA0001697936770000051
速度模型右边所对应的吸收边界方程为：
Figure BDA0001697936770000052
其中，
Figure BDA0001697936770000053
u、w分别为弹性波场所对应的水平位移和垂直位移；
Figure BDA0001697936770000054
Figure BDA0001697936770000055
Figure BDA0001697936770000056
r1，r2、r3，r4为波场吸收系数，1≥r1≥0；1≥r2≥0，1≥r3≥0；1≥r4≥0，α、β分别为纵、横波速度。
优选地，通过方程(5a)-(5d)计算所述速度模型四个角点所对应的各个时刻的波场值，
右下角角点的波场计算方程为：
Figure BDA0001697936770000057
左上角角点的波场计算方程为：
Figure BDA0001697936770000058
右上角角点的波场计算方程为：
Figure BDA0001697936770000059
左下角角点的波场计算方程为：
Figure BDA00016979367700000510
其中：
Figure BDA00016979367700000511
u、w分别为弹性波场所对应的水平位移场和垂直位移场；
Figure BDA0001697936770000061
Figure BDA0001697936770000062
α、β分别为纵、横波速度，x和z是二维介质情况下沿X轴和Z轴的坐标；
其中：
Figure BDA0001697936770000063
Figure BDA0001697936770000064
θ为波的入射角。
根据本发明的另一方面，提出了一种弹性波场数值模拟系统，其上存储有计算机程序，其中，所述程序被处理器执行时实现以下步骤：步骤1：建立速度模型，确定所述速度模型的空间采样间隔和时间采样间隔，对所述速度模型进行有限差分网格离散化；步骤2：根据弹性波动方程计算所述速度模型的内部各个网格点的各个时刻的波场值，获得内部波场数据；步骤3：根据角点波场数值计算公式计算所述速度模型的四个角点所对应的各个时刻的波场值，获得角点波场数据；步骤4：根据波场吸收边界方程和波场吸收系数计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据；步骤5：根据所述边界波场数据、所述内部波场数据、所述角点波场数据，获得所述速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5；步骤6：输出所述速度模型的最终波场数据。
优选地，所述速度模型的内部各个网格点的各个时间的波场值的表达式为：
Figure BDA0001697936770000065
其中，
Figure BDA0001697936770000066
u、w分别为弹性波场所对应的水平位移和垂直位移；
Figure BDA0001697936770000071
α,β分别为纵、横波速度，t为时间，x和z是二维介质情况下沿X轴和Z轴的坐标。
优选地，所述波场吸收边界方程通过以下方式获得：对所述弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程；将所述傅里叶变换后的弹性波动方程进行因式分解与化简，获得化简方程；将所述化简方程转换至时间域，获得波场吸收边界方程。
优选地，傅里叶变换后的弹性波动方程为：
Figure BDA0001697936770000072
优选地，所述化简方程为：
Figure BDA0001697936770000073
优选地，通过方程(4a)-(4d)计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，
速度模型底边所对应的吸收边界方程为：
Figure BDA0001697936770000074
速度模型顶边所对应的吸收边界方程为：
Figure BDA0001697936770000075
速度模型左边所对应的吸收边界方程为：
Figure BDA0001697936770000076
速度模型右边所对应的吸收边界方程为：
Figure BDA0001697936770000077
其中，
Figure BDA0001697936770000081
u、w分别为弹性波场所对应的水平位移和垂直位移；
Figure BDA0001697936770000082
Figure BDA0001697936770000083
Figure BDA0001697936770000084
r1，r2、r3，r4为波场吸收系数，1≥r1≥0；1≥r2≥0，1≥r3≥0；1≥r4≥0，α、β分别为纵、横波速度。
优选地，通过方程(5a)-(5d)计算所述速度模型四个角点所对应的各个时刻的波场值，
右下角角点的波场计算方程为：
Figure BDA0001697936770000085
左上角角点的波场计算方程为：
Figure BDA0001697936770000086
右上角角点的波场计算方程为：
Figure BDA0001697936770000087
左下角角点的波场计算方程为：
Figure BDA0001697936770000088
其中：
Figure BDA0001697936770000089
u、w分别为弹性波场所对应的水平位移场和垂直位移场；
Figure BDA00016979367700000810
Figure BDA00016979367700000811
α、β分别为纵、横波速度，x和z是二维介质情况下沿X轴和Z轴的坐标；
其中：
Figure BDA0001697936770000091
Figure BDA0001697936770000092
θ为波的入射角。
本发明具有其它的特性和优点，这些特性和优点从并入本文中的附图和随后的具体实施方式中将是显而易见的，或者将在并入本文中的附图和随后的具体实施方式中进行详细陈述，这些附图和具体实施方式共同用于解释本发明的特定原理。
附图说明
通过结合附图对本发明示例性实施例进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显，其中，在本发明示例性实施例中，相同的参考标号通常代表相同部件。
图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：输出速度模型的最终波场数据。
在一个示例中，速度模型的内部各个网格点的各个时间的波场值的表达式为：
Figure BDA0001697936770000101
其中，
Figure BDA0001697936770000102
u,w分别为弹性波场所对应的水平位移和垂直位移；
Figure BDA0001697936770000103
Figure BDA0001697936770000111
α,β分别为纵、横波速度，t为时间。
在一个示例中，波场吸收边界方程通过以下方式获得：对弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程；将傅里叶变换后的弹性波动方程进行因式分解与化简，获得化简方程；将化简方程转换至时间域，获得波场吸收边界方程。
在一个示例中，傅里叶变换后的弹性波动方程为：
Figure BDA0001697936770000112
在一个示例中，化简方程为：
Figure BDA0001697936770000113
在一个示例中，通过方程(4a)-(4d)计算速度模型的边界上的网格点所对应的各个时刻的波场值，
速度模型底边所对应的吸收边界方程为：
Figure BDA0001697936770000114
速度模型顶边所对应的吸收边界方程为：
Figure BDA0001697936770000115
速度模型左边所对应的吸收边界方程为：
Figure BDA0001697936770000116
速度模型右边所对应的吸收边界方程为：
Figure BDA0001697936770000117
其中，
Figure BDA0001697936770000118
u、w分别为弹性波场所对应的水平位移和垂直位移；
Figure BDA0001697936770000121
Figure BDA0001697936770000122
Figure BDA0001697936770000123
r1，r2、r3，r4为波场吸收系数，1≥r1≥0；1≥r2≥0，1≥r3≥0；1≥r4≥0，α、β分别为纵、横波速度。
在一个示例中，通过方程(5a)-(5d)计算速度模型四个角点所对应的各个时刻的波场值，
右下角角点的波场计算方程为：
Figure BDA0001697936770000124
左上角角点的波场计算方程为：
Figure BDA0001697936770000125
右上角角点的波场计算方程为：
Figure BDA0001697936770000126
左下角角点的波场计算方程为：
Figure BDA0001697936770000127
其中：
Figure BDA0001697936770000128
u、w分别为弹性波场所对应的水平位移场和垂直位移场；
Figure BDA0001697936770000129
Figure BDA00016979367700001210
α、β分别为纵、横波速度，x和z是二维介质情况下沿X轴和Z轴的坐标；
其中：
Figure BDA00016979367700001211
Figure BDA0001697936770000131
θ为波的入射角。
具体地，根据本发明的弹性波场数值模拟方法可以包括：
步骤1：建立速度模型，确定炮点(即激发源)位置，选择地震子波，确定速度模型的空间采样间隔和时间采样间隔，对速度模型进行有限差分网格离散化。
步骤2：在各向同性介质条件下，假设x和z是二维介质情况下沿X轴和Z轴的坐标，X轴的方向向右，Z轴的方向朝下，我们可以利用两个共轭的二阶偏微分方程来描述介质中P波的运动以及垂向偏振的SV波的运动，在这里不考虑沿水平方向偏振的SH波，u、w分别为水平方向与垂直方向的位移，ρ为介质的密度，t为时间，λ和μ为具体介质中的拉梅系数，那么可以得到关于各向同性非均匀介质的弹性波动方程为公式(5)：
Figure BDA0001697936770000132
假定密度ρ为常数，这样就可以把公式(6)看成是随空间位置变化的P波和SV波速度的函数。其中λ和μ与纵、横波速度α(x,z)和β(x,z)的关系为公式(6)：
Figure BDA0001697936770000133
进而将公式(6)改写为矩阵形式，即为公式(1)，计算速度模型的内部各个网格点的各个时刻的波场值为公式(1)，获得内部波场数据。
步骤3：根据角点波场数值计算公式计算速度模型的四个角点的波场值分别为公式(5a)-(5d)，获得角点波场数据，根据炮点位置坐标以及四个角点所对应的位置坐标，计算出sinθ和cosθ，由此得到J矩阵向量中的各项值。
步骤4：根据弹性波动方程与波场吸收系数获得波场吸收边界方程，根据波场吸收边界方程计算对应波场吸收系数的速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据，包括：
对弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程为公式(2)；为了能对公式(2)的左边项进行分解，考虑将公式(2)的左边项变为公式(8)：
Figure BDA0001697936770000141
其中：
Figure BDA0001697936770000142
对(8)式进行整理，得到公式(9)：
Figure BDA0001697936770000143
由(9)式得到(10)式：
Figure BDA0001697936770000144
令：
Figure BDA0001697936770000145
则有公式(11)：
Figure BDA0001697936770000146
进而获得公式(12)：
Figure BDA0001697936770000151
其中：1≥r1≥0；1≥r2≥0。
于是，(10)式可以写成下面的推导式子：
Figure BDA0001697936770000152
可以将(13)式可改写成下面的形式：
I+C(kz/ω)-A(kx/ω)-B(kz/ω)＝0 (14)
其中：
Figure BDA0001697936770000153
这样，通过因式分解，公式(2)便可以将改写成下面的形式：
Figure BDA0001697936770000154
令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：输出速度模型的最终波场数据。
在一个示例中，速度模型的内部各个网格点的各个时间的波场值的表达式为：
Figure BDA0001697936770000191
其中，
Figure BDA0001697936770000192
u、w分别为弹性波场所对应的水平位移和垂直位移；
Figure BDA0001697936770000193
α,β分别为纵、横波速度，t为时间。
在一个示例中，波场吸收边界方程通过以下方式获得：对弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程；将傅里叶变换后的弹性波动方程进行因式分解与化简，获得化简方程；将化简方程转换至时间域，获得波场吸收边界方程。
在一个示例中，傅里叶变换后的弹性波动方程为：
Figure BDA0001697936770000194
在一个示例中，化简方程为：
Figure BDA0001697936770000195
在一个示例中，通过方程(4a)-(4d)计算速度模型的边界上的网格点所对应的各个时刻的波场值，
速度模型底边所对应的吸收边界方程为：
Figure BDA0001697936770000201
速度模型顶边所对应的吸收边界方程为：
Figure BDA0001697936770000202
速度模型左边所对应的吸收边界方程为：
Figure BDA0001697936770000203
速度模型右边所对应的吸收边界方程为：
Figure BDA0001697936770000204
其中，
Figure BDA0001697936770000205
u、w分别为弹性波场所对应的水平位移和垂直位移；
Figure BDA0001697936770000206
Figure BDA0001697936770000207
Figure BDA0001697936770000208
r1，r2、r3，r4为波场吸收系数，1≥r1≥0；1≥r2≥0，1≥r3≥0；1≥r4≥0，α、β分别为纵、横波速度。
在一个示例中，通过方程(5a)-(5d)计算速度模型四个角点所对应的各个时刻的波场值，
右下角角点的波场计算方程为：
Figure BDA0001697936770000209
左上角角点的波场计算方程为：
Figure BDA00016979367700002010
右上角角点的波场计算方程为：
Figure BDA0001697936770000211
左下角角点的波场计算方程为：
Figure BDA0001697936770000212
其中：
Figure BDA0001697936770000213
u、w分别为弹性波场所对应的水平位移场和垂直位移场；
Figure BDA0001697936770000214
Figure BDA0001697936770000215
α、β分别为纵、横波速度，x和z是二维介质情况下沿X轴和Z轴的坐标；
其中：
Figure BDA0001697936770000216
Figure BDA0001697936770000217
θ为波的入射角。
本系统通过不断调整波场吸收系数，有效地消除来自有限的速度模型边界所产生的反射能量，当吸收系数调整到最佳状态时，则可以完全吸收来自模型边界的反射能量，从而使计算得到的模型波场值能准确地反映出弹性波在二维各向同性介质中传播的真实性，进而达到提高弹性波数值模拟结果精度的目的。
以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。

Claims (6)

1.一种弹性波场数值模拟方法，包括：
步骤1：建立速度模型，确定所述速度模型的空间采样间隔和时间采样间隔，对所述速度模型进行有限差分网格离散化；
步骤2：根据弹性波动方程计算所述速度模型的内部各个网格点的各个时刻的波场值，获得内部波场数据；
步骤3：根据角点波场数值计算公式计算所述速度模型的四个角点所对应的各个时刻的波场值，获得角点波场数据；
步骤4：根据波场吸收边界方程和波场吸收系数计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据；
步骤5：根据所述边界波场数据、所述内部波场数据、所述角点波场数据，获得所述速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5；
步骤6：输出所述速度模型的最终波场数据；
其中，所述波场吸收边界方程通过以下方式获得：
对所述弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程；
将所述傅里叶变换后的弹性波动方程进行因式分解与化简，获得化简方程；
将所述化简方程转换至时间域，获得波场吸收边界方程；
其中，通过方程(4a)-(4d)计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，
速度模型底边所对应的吸收边界方程为：
Figure FDA0002728433570000011
速度模型顶边所对应的吸收边界方程为：
Figure FDA0002728433570000021
速度模型左边所对应的吸收边界方程为：
Figure FDA0002728433570000022
速度模型右边所对应的吸收边界方程为：
Figure FDA0002728433570000023
其中，
Figure FDA0002728433570000024
u、w分别为弹性波场所对应的水平位移和垂直位移；
Figure FDA0002728433570000025
Figure FDA0002728433570000026
Figure FDA0002728433570000027
r1，r2、r3，r4为波场吸收系数，1≥r1≥0；1≥r2≥0，1≥r3≥0；1≥r4≥0，α、β分别为纵、横波速度，t为时间，x和z是二维介质情况下沿X轴和Z轴的坐标。
2.根据权利要求1所述的弹性波场数值模拟方法，其中，所述速度模型的内部各个网格点的各个时间的波场值的表达式为：
Figure FDA0002728433570000028
其中，
Figure FDA0002728433570000029
u、w分别为弹性波场所对应的水平位移和垂直位移；
Figure FDA0002728433570000031
α,β分别为纵、横波速度。
3.根据权利要求1所述的弹性波场数值模拟方法，其中，傅里叶变换后的弹性波动方程为：
Figure FDA0002728433570000032
4.根据权利要求1所述的弹性波场数值模拟方法，其中，所述化简方程为：
Figure FDA0002728433570000033
5.根据权利要求1所述的弹性波场数值模拟方法，其中，通过方程(5a)-(5d)计算所述速度模型四个角点所对应的各个时刻的波场值：
右下角角点的波场计算方程为：
Figure FDA0002728433570000034
左上角角点的波场计算方程为：
Figure FDA0002728433570000035
右上角角点的波场计算方程为：
Figure FDA0002728433570000036
左下角角点的波场计算方程为：
Figure FDA0002728433570000037
其中：
Figure FDA0002728433570000038
u、w分别为弹性波场所对应的水平位移场和垂直位移场；
Figure FDA0002728433570000039
Figure FDA00027284335700000310
α、β分别为纵、横波速度，x和z是二维介质情况下沿X轴和Z轴的坐标；
其中：
Figure FDA0002728433570000041
Figure FDA0002728433570000042
θ为波的入射角。
6.一种弹性波场数值模拟系统，其上存储有计算机程序，其中，所述程序被处理器执行时实现以下步骤：
步骤1：建立速度模型，确定所述速度模型的空间采样间隔和时间采样间隔，对所述速度模型进行有限差分网格离散化；
步骤2：根据弹性波动方程计算所述速度模型的内部各个网格点的各个时刻的波场值，获得内部波场数据；
步骤3：根据角点波场数值计算公式计算所述速度模型的四个角点所对应的各个时刻的波场值，获得角点波场数据；
步骤4：根据波场吸收边界方程和波场吸收系数计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据；
步骤5：根据所述边界波场数据、所述内部波场数据、所述角点波场数据，获得所述速度模型的波场数据，判断波场边界的反射情况是否满足预定标准，若是，则执行步骤6，若否，则重新确定波场吸收系数，重复步骤4-5；
步骤6：输出所述速度模型的最终波场数据；
其中，根据所述弹性波动方程获得波场吸收边界方程，根据所述波场吸收边界方程计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据包括：
对所述弹性波动方程进行傅里叶变换，获得傅里叶变换后的弹性波动方程；
将所述傅里叶变换后的弹性波动方程进行因式分解与化简，获得化简方程；
将所述化简方程转换至时间域，获得波场吸收边界方程；
根据所述波场吸收边界方程计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，获得边界波场数据；
其中，通过方程(4a)-(4d)计算所述速度模型的边界上的网格点所对应的各个时刻的波场值，
速度模型底边所对应的吸收边界方程为：
Figure FDA0002728433570000051
速度模型顶边所对应的吸收边界方程为：
Figure FDA0002728433570000052
速度模型左边所对应的吸收边界方程为：
Figure FDA0002728433570000053
速度模型右边所对应的吸收边界方程为：
Figure FDA0002728433570000054
其中，
Figure FDA0002728433570000055
u、w分别为弹性波场所对应的水平位移和垂直位移；
Figure FDA0002728433570000056
Figure FDA0002728433570000057
Figure FDA0002728433570000058
r1，r2、r3，r4为波场吸收系数，1≥r1≥0；1≥r2≥0，1≥r3≥0；1≥r4≥0，α、β分别为纵、横波速度，t为时间，x和z是二维介质情况下沿X轴和Z轴的坐标。
CN201810620610.XA 2018-06-15 2018-06-15 弹性波场数值模拟方法及系统 Active CN110609325B (zh)

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 CN110609325A (zh) 2019-12-24
CN110609325B true 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)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111257930B (zh) * 2020-02-18 2022-03-29 中国石油大学(华东) 一种黏弹各向异性双相介质区域变网格求解算子

Family Cites Families (5)

* Cited by examiner, † Cited by third party
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
SG11201704620WA (en) * 2015-02-13 2017-09-28 Exxonmobil Upstream Res Co Efficient and stable absorbing boundary condition in finite-difference calculations
CN107315192B (zh) * 2016-04-26 2019-07-05 中国石油化工股份有限公司 基于二维各向同性介质的弹性波波场数值的模拟方法
CN105807317B (zh) * 2016-05-06 2019-04-30 中国地质大学（北京） 基于切比雪夫伪谱法的各向异性衰减面波模拟方法
CN107942375B (zh) * 2017-11-17 2019-04-30 河海大学 一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法

Also Published As

Publication number Publication date
CN110609325A (zh) 2019-12-24

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
Gazdag et al. Migration of seismic data by phase shift plus interpolation
Robertsson A numerical free-surface condition for elastic/viscoelastic finite-difference modeling in the presence of topography
Ladopoulos Petroleum and gas reserves exploration by real-time expert seismology and non-linear seismic wave motion'
US20170299745A1 (en) Prestack egs migration method for seismic wave multi-component data
US20210103065A1 (en) Determining properties of a subterranean formation using an acoustic wave equation with a reflectivity parameterization
CN108108331A (zh) 一种基于拟空间域弹性波方程的有限差分计算方法
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
Wu et al. Seismic wave propagation and scattering in heterogeneous crustal waveguides using screen propagators: I SH waves
CN107390270A (zh) 一种基于弹性波逆时偏移ADCIGs的AVA分析方法
Fu et al. Infinite boundary element absorbing boundary for wave propagation simulations
Golubev et al. Simulation of seismic responses from fractured MARMOUSI2 model
CN110609325B (zh) 弹性波场数值模拟方法及系统
US10379245B2 (en) Method and system for efficient extrapolation of a combined source-and-receiver wavefield
Zhao et al. Frequency-domain finite-difference elastic wave modeling in the presence of surface topography
US9791580B2 (en) Methods and systems to separate wavefields using pressure wavefield data
CN107315192B (zh) 基于二维各向同性介质的弹性波波场数值的模拟方法
Nakagawa et al. Three-dimensional elastic wave scattering by a layer containing vertical periodic fractures
Bader Seismic and Newtonian noise modeling for advanced Virgo and Einstein telescope
Meng et al. A 3D Optimized Frequency–Wavenumber (FK), Time–Space Optimized Symplectic (TSOS) Hybrid Method for Teleseismic Wave Modeling
Gao‐Xiang et al. A quantitative analysis method for the seismic geological complexity of near surface
Lee et al. Elastic full-waveform inversion using both the multiparametric approximate hessian and the discrete cosine transform
Li et al. The numerical modeling of frequency domain scalar wave equation based on the rhombus stencil

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