CN107976710B - 一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法 - Google Patents
一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法 Download PDFInfo
- Publication number
- CN107976710B CN107976710B CN201711142912.2A CN201711142912A CN107976710B CN 107976710 B CN107976710 B CN 107976710B CN 201711142912 A CN201711142912 A CN 201711142912A CN 107976710 B CN107976710 B CN 107976710B
- Authority
- CN
- China
- Prior art keywords
- time
- difference
- wave
- implicit
- order
- 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.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 95
- 238000005457 optimization Methods 0.000 title claims abstract description 26
- 238000004088 simulation Methods 0.000 title claims abstract description 22
- 239000006185 dispersion Substances 0.000 claims abstract description 52
- 238000010521 absorption reaction Methods 0.000 claims abstract description 11
- 229910003460 diamond Inorganic materials 0.000 claims abstract description 10
- 239000010432 diamond Substances 0.000 claims abstract description 10
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 5
- 238000004364 calculation method Methods 0.000 claims description 17
- 230000007704 transition Effects 0.000 claims description 8
- 238000004458 analytical method Methods 0.000 claims description 4
- 230000008569 process Effects 0.000 claims description 4
- 238000005070 sampling Methods 0.000 claims description 4
- 230000010354 integration Effects 0.000 claims description 3
- 238000012886 linear function Methods 0.000 claims description 3
- 230000006870 function Effects 0.000 description 13
- 239000012088 reference solution Substances 0.000 description 11
- 150000003839 salts Chemical class 0.000 description 8
- 238000007796 conventional method Methods 0.000 description 5
- 230000005012 migration Effects 0.000 description 4
- 238000013508 migration Methods 0.000 description 4
- 239000000243 solution Substances 0.000 description 3
- 230000007547 defect Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000010835 comparative analysis Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000002123 temporal effect 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
-
- 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/51—Migration
-
- 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
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/70—Other details related to processing
- G01V2210/74—Visualisation of seismic data
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)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
- Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
- Secondary Cells (AREA)
Abstract
本发明公开了一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法，包括以下步骤：(1)读取参数；(2)基于菱形差分算子和二阶时间中心差分，得到时间导数的时间高阶离散格式；(3)采用空间隐式离散格式求解二阶空间导数，基于空间频散关系，优化方法求解隐式差分系数；(4)基于时间高阶离散格式和空间隐式离散格式，得到时空域频散关系；(5)采用线性优化算法对时间高阶离散格式中的差分系数进行求解；(6)利用求取的差分系数，采用混合吸收边界条件对边界反射进行吸收，按照波动方程进行递推，得到任意时刻的波场及整个地震记录；(7)记录波场快照，输出地震记录并结束。本发明的有限差分数值模拟方法精度高，稳定性好，频散误差小。
Description
技术领域
本发明属于地震波场数值模拟技术领域，涉及一种基于求解声波方程的高精度、高效率线性优化隐式时空域有限差分数值模拟方法。
背景技术
目前逆时偏移方法和全波形反演技术在处理如盐丘等复杂构造的成像和储层参数反演方面已经显示出了巨大的潜力，越来越受到地球物理学家的重视。但逆时偏移和全波形反演方法都涉及波动方程的数值求解，在模型偏大尤其三维时，计算量大，耗时长，某种程度上制约着这两种技术的发展。因此高精度、高效率的数值模拟方法就显得尤其重要。
有限差分法相比其他数值模拟方法，具有操作简单、计算效率高、易于并行和所需内存小等诸多优势，在逆时偏移和全波形反演中被广泛应用。但目前有限差分法存在频散严重、精度低和稳定性差等缺点。频散按照离散对象可以分为时间差分频散和空间差分频散。
发明内容
本发明的目的是为克服传统隐式空间有限差分方法时间精度低、稳定性差和易于频散等问题而提供一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法，具有更高精度和更好稳定性，本发明在保持空间隐式差分精度的前提下，提出了显著压制时间频散的线性优化时空域策略(即权利请求书中的第5步骤)，可以为逆时偏移和全波形反演方法提供频散更小、精度更高的波场，同时可以帮助提高声波方程数值模拟的精度和效率。
本发明的一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法，包括以下几个步骤：
(1)读取参数，所述参数包括正演模拟所需要的速度模型参数文件、有限差分算子长度、子波函数及子波主频、正演所采用的时间与空间步长及地震记录时长；
(2)基于菱形差分算子和二阶时间中心差分，得到时间导数的时间高阶离散格式；
(3)借鉴Claerbout(1985)提出的求解二阶空间导数隐式差分格式，在常规显式高阶差分的基础上，分母中引入二阶中心差分格式，构造空间隐式差分格式求解二阶空间导数，用于减小有限差分算子长度；进一步，基于最小二乘法，以L2范数目标函数极小化空间域频散关系，求解优化的隐式差分系数；
(4)基于所述时间高阶离散格式和空间隐式离散格式，带入声波方程，构建同时具有时间高阶和空间隐式的差分递推格式，得到时空域频散关系；
(5)固定空间差分系数，采用优化策略极小化整个递推格式频散关系，建立关于时间差分格式中差分系数的L2范数目标函数，利用线性优化算法对时间差分系数进行求解(详细内容见下文公式(11)～(16))；
(6)采用Liu等(2010)提出的混合吸收边界条件对边界反射能量进行吸收，按照上述介绍的时间高阶和空间隐式差分格式离散波动方程并进行递推，得到任意时刻的波场及整个地震记录；
(7)记录波场快照，输出地震记录并结束。
步骤(2)中，求解二阶时间导数时采用中心差分，形式如下：
其中，为压力场，τ为时间步长，t为时间；为了进一步提高时间离散精度，将菱形差分算子引入到公式(1)中，得出如下的时间高阶离散公式：
其中，h为空间采样步长，v为速度，am,n为有限差分系数；N为菱形算子中的差分算子长度，其值越大，时间差分近似精度越高，但增加的计算量也增大；n和m为菱形算子中的坐标位置序列。
步骤(3)中，采用如下空间隐式离散格式对二阶空间导数进行近似，形式如下：
其中，a′0、a′m和b为隐式差分系数，当b＝0时，公式(3)退化为显式差分格式；h为空间步长，x为方向坐标，M为差分算子长度；公式(3)对应的空间频散关系如下：
其中，β＝kxh，kx为水平方向波数；极小化公式(4)对应的空间频散关系，目标函数为关于差分系数的线性凸函数，最优差分系数通过求解如下公式得到：
其中，n＝1,2,…,M+1，d＝[d1,d2,…,dM,dM+1]T＝[a′1,a′2,…,a′M-1,a′M,b]T，q为介于0和1之间的自然数，控制着积分范围，
步骤(4)中，在求解时间导数时，采用公式(2)计算；求解空间导数时，采用公式(3)计算，将公式(2)和公式(3)带入到声波方程并进行平面波理论分析，过程如下：
均匀模型中，声波方程压力P满足如下平面波方程，
其中，ω是角频率，是虚数单位，kz为垂直方向波数；m、n和l为位置和时间序列；
将公式(7)带入公式(2)并进行化简得到如下公式：
同样，对空间导数公式进行分析，得到如下频散关系：
将上述(8)、(9)和(10)带入声波方程并简化和重组得到如下时空域频散关系，
其中，ω为角频率，kx和kz分别为沿x和z方向的波数，r为库朗数，r＝vτ/h。
步骤(5)中，首先，差分系数a′m和b由公式(5)计算得到，在计算差分系数am,0和am,n时保持固定，此时优化差分系数am,0和am,n是线性优化问题，直接线性求解，定义如下目标函数：
其中，d＝é[d1,…,dN+(N-1)N/2]T＝[a1,0,…,aN,0,a1,1,…,a1,N-1,…,aN-1,1]T，表示差分系数；θ为传播角度，q为介于0和1间的实数，波数kx和kz与波数k之间满足如下关系：
kx＝kcos(θ),kz＝ksin(θ),(13)ψ(β,θ,r)及定义如下
目标函数是关于差分系数am,0和am,n的线性函数，差分系数am,0和am,n通过求解如下公式得到，
步骤(6)中，整个计算区域被划分为内部区域、边界区域和过渡区域，具体的方法如下：
(6-1)利用公式(2)和公式(3)离散时间和空间导数，采用公式(5)和公式(16)计算差分系数，进行双程波方程求解，求解声波方程，计算所有区域处的双程波波场；
(6-2)求解单程波方程，计算边界区域和过渡区域处的单程波波场，计算公式及方法如下：
针对二阶声波方程，以左边界和左上角为例，左边界单程波方程为：
其中，x和z代表坐标；
左上角单程波方程为：
其他边界处及角点处单程波可以类似得到，并进行对应求解；采用上述吸收边界条件策略计算各位置处任意时刻波场值；
(6-3)计算最终波场：
边界处波场为单程波场值，内部区域处波场为双程波场值，中间过渡区域处波场通过如下线性组合得到：
P＝BiPtwo+(1-Bi)Pone, (19)
其中，Ptwo为双程波场，Pone为单程波场，Bi为权重系数，Bi＝(i-1)/Nb，i＝2,…,Nb，Nb为过渡区域吸收边界层数，过渡区域从内到外，Bi由(Nb-1)/Nb变化到1/Nb；如果Nb＝1，则此时混合吸收边界条件退化为单程波吸收边界条件。
本发明与现有技术相比其显著优点在于：
(a)时间精度高，频散小。在使用相同时间步长且步长偏大时，本发明方法可以得到精度更高的波场。
(b)稳定性好，模拟中可以采用更大的时间步长。
(c)时间精度提高的同时，空间隐式精度很好地保持。
附图说明
图1为本发明的一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法工作流程示意图；
图2(a)为常规隐式空域方法频散曲线随时间步长变化图；
图2(b)为本发明方法当N＝2时频散曲线随时间步长变化图；
图2(c)为本发明方法当N＝3时频散曲线随时间步长变化图；
图3为本发明方法稳定性条件曲线与常规方法对比图；
图4(a)为均匀模型常规隐式差分方法模拟震动曲线同参考解对比图；
图4(b)为均匀模型本发明方法在N＝2时模拟震动曲线同参考解对比图；
图4(c)为本发明方法在N＝3时模拟震动曲线同参考解对比图；
图4(d)为本发明方法与常规隐式差分法模拟结果相对于参考解均方根误差对比图；
图5为SEG/EAGE盐丘模型；
图6(a)为SEG/EAGE盐丘模型2.4s时刻参考解波场快照；
图6(b)为SEG/EAGE盐丘模型应用常规隐式差分方法后2.4s时刻波场快照；
图6(c)为SEG/EAGE盐丘模型应用本发明方法后2.4s时刻波场快照；
图7(a)为SEG/EAGE盐丘模型应用常规隐式差分方法后(500m,4500m)处地震记录震动图；
图7(b)为SEG/EAGE盐丘模型应用本发明方法后(500m,4500m)处地震记录震动图。
具体实施方式
下面结合附图和实施例对本发明的具体实施方式作进一步的详细说明。
本发明提供了一种针对声波方程的基于线性优化方法的隐式时空域有限差分数值模拟方法。为解决常规隐式空间差分方法时间精度低、易于时间频散的不足，该方法在常规时间二阶离散格式的基础上，加入了菱形差分算子。结合高阶隐式空间差分方法，本发明进一步得出了同时具有高阶时间精度和隐式空间精度的递推格式，并推导了相应递推格式的时空域频散关系。为进一步压制时间频散，本发明从时空域频散关系出发，采用线性最小二乘优化算法，提出了一种计算时空域差分系数的方法和策略。相比常规隐式差分方法，本发明精度更高，方法更加稳定且计算效率高。
如图1所示，本发明提出的一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法包括如下具体步骤：
(1)读取参数，包括正演模拟所需要的速度模型参数文件、有限差分算子长度(时间导数和空间导数)、子波函数及子波主频、正演所采用的时间和空间步长及地震记录时长等参数；
(2)针对时间导数的时间高阶离散格式构造；
传统方法在求解2阶时间导数时采用中心差分，形式如下：
其中，为压力场，τ为时间步长，t为时间。为了进一步提高时间离散精度，本发明将菱形差分算子引入到公式(1)中，提出了如下的时间高阶离散公式：
其中，h为空间采样步长，v为速度，am,n为有限差分系数。N为菱形算子中的差分算子长度，其值越大，时间差分近似精度越高，但增加的计算量也增大。n和m为菱形算子中的坐标位置序列。(3)针对空间导数，采用隐式差分格式求解以减小有限差分算子长度。为进一步减小频散，基于空间频散关系，通过优化方法求解隐式差分系数；
针对二阶空间导数，为有效减小算子长度，本发明采用如下空间隐式离散格式对空间二阶导数进行近似，形式如下：
其中，a′0、a′m和b为隐式差分系数，h为空间步长，x为方向坐标，M为差分算子长度。常规方法在求取差分系数a′m和b时采用泰勒展开法，精度低，易于频散。本文采用优化方法求解，公式(3)对应的空间频散关系如下：
其中，β＝kxh，kx为水平方向波数。极小化公式(4)左右两边的绝对误差，目标函数为关于差分系数的线性凸函数，可以通过求解如下公式得到：
其中，d＝[d1,d2,…,dM,dM+1]T＝[a′1,a′2,…,a′M-1,a′M,b]T，q为介于0和1之间的自然数，控制着积分范围，
(4)结合时间高阶离散格式和空间隐式离散格式，构建同时具有时间高阶和空间隐式的差分递推格式，得到时空域频散关系；
在求解时间导数时，采用公式(2)计算，求解空间导数时，采用公式(3)计算，将公式(2)和公式(3)带入到声波方程并进行平面波分析，可以得到如下时空域频散关系，
其中，ω为角频率，kx和kz分别为沿x和z方向的波数，r为库朗数，r＝vτ/h。
(5)从时空域频散关系出发，采用线性优化算法对时空域差分系数进行求解；
不同于显式差分方法，隐式差分方法频散关系形式更为复杂，直接同时优化am,0、am,n、a′m和b是强非线性问题。
首先，差分系数a′m和b由公式(5)计算得到，在计算差分系数am,0和am,n时保持固定，此时优化差分系数am,0和am,n是线性优化问题，可以直接线性求解，定义如下目标函数：
其中，d＝[d1,…,dN+(N-1)N/2]T＝[a1,0,…,aN,0,a1,1,…,a1,N-1,…,aN-1,1]T，表示差分系数，θ为传播角度，q为介于0和1间的实数，
目标函数(8)是关于差分系数am,0和am,n的线性函数，差分系数am,0和am,n可以通过求解如下公式得到，
(6)利用求取的差分系数，采用适当的吸收边界条件，进行波场延拓；
求得差分系数后，按照波动方程进行递推，得到任意时刻的波场及整个地震记录。波场延拓中必须要考虑的一个问题是边界反射。本发明采用Liu和Sen(2010)提出的混合吸收边界条件对边界反射进行吸收。
(7)记录波场快照，输出地震记录并终止程序。
以下为本发明的一个实施例，说明基于一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法的实现过程。为了方便区分本发明方法和常规方法，以下内容中，将常规隐式差分方法记做ISD，将本发明方法记做ITSD-1。
图2a至2c为本发明提出的方法与传统方法频散误差对比分析，图中纵坐标为数值模拟速度与真实速度比值，反映了方法的精度。当纵坐标值为1.0时，此时方法无频散。纵坐标偏离1.0的位置越远，则频散越大。图中，速度参数为1500m/s，空间采样步长为15m，算子长度M＝6。图中显示的是传播方向为45度时采用不同时间步长时频散曲线。对比发现，相比传统方法，本发明中提出的方法不同程度的改进了传统方法的精度，压制了数值频散。随着参数N的增大，本发明中提出的方法精度提高，频散变小。
图3为本发明提出的方法与传统方法稳定性对比分析图，图中纵坐标为稳定性因子，当库朗数大于该稳定性因子时，差分方法不稳定。该值越大，方法越稳定。对比发现，传统方法稳定性最差，稳定性因子最小。本发明提出的方法稳定性更好，且随着N的增大，稳定性越来越好。
图4(a)至4(c)为一均匀模型(1200m,2250m)处模拟波场震动图与参考解对比图。模型大小为4500m×4500m，空间步长为15m，时间步长为2ms，采用的子波为主频是18Hz的雷克子波，在模型中心激发。参考解为由Cagniard-de Hoop(De Hoop,1960)方法计算得到。本发明中方法采用M＝6。对比发现，传统方法频散误差明显，与参考解吻合差。本发明提出的方法则显著提高了精度，N＝3时可以得到比N＝2更高的精度。图4(d)为各方法与参考解之间的统计均方根误差，明显本发明提出的方法误差要小。
图5为地震勘探数值模拟中常用到的SEG/EAGE盐丘模型。模型大小为5000m×15000m，空间步长为25m，采用的子波为主频是15Hz的雷克子波，在地表中心激发，时间步长为2.0ms。本发明中方法采用M＝8和N＝2。图6(a)至6(c)为各方法在2.4s时刻时传播波场快照，参考解为常规显式时空域方法采用很小的时间步长和很大的算子长度计算得到。同参考解对比发现，本发明中方法明显具有更高的精度，白色箭头所指区域频散相比传统方法明显小。图7(a)和(b)为各方法在(500m,4500m)处地震记录震动图相比参考解对比图。图中显示的误差值为统计的均方根误差，明显本发明提出的方法均方根误差更小。
为降低频散误差，传统隐式差分方法采用1ms时间步长时计算的均方根误差为1.1824e-4，消耗的CPU时间为403.58s。本发明提出的方法(ITSD-1)采用2ms时间步长时计算的均方根误差为1.2390e-4，消耗的CPU时间为295.27s。可以发现，本发明提出的方法在基本保持了精度的前提下，显著降低了计算时间，效率得到大大提高。
本方法能够显著提高常规模拟方法的精度和效率，能够为计算量巨大的逆时偏移方法和全波形反演方法高效率地提供更为准确的波场信息，具有重要的实际用途。本发明的具体实施方式中未涉及的说明属于本领域公知的技术，可参考公知技术加以实施。
以上具体实施方式及实施例是对本发明提出的一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法技术思想的具体支持，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在本技术方案基础上所做的任何等同变化或等效的改动，均仍属于本发明技术方案保护的范围。
Claims (2)
1.一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法，其特征在于，包括以下几个步骤：
(1)读取参数，所述参数包括正演模拟所需要的速度模型参数文件、有限差分算子长度、子波函数及子波主频、正演所采用的时间与空间步长及地震记录时长；
(2)基于菱形差分算子和二阶时间中心差分，得到时间导数的时间高阶离散格式；
(3)求解二阶空间导数隐式差分格式，在显式高阶差分的基础上，分母中引入二阶中心差分格式，构造空间隐式差分格式求解二阶空间导数，用于减小有限差分算子长度；进一步，基于最小二乘法，以L2范数目标函数极小化空间域频散关系，求解优化的隐式差分系数；
(4)基于所述时间高阶离散格式和空间隐式离散格式，带入声波方程，构建同时具有时间高阶和空间隐式的差分递推格式，得到时空域频散关系；
(5)固定空间差分系数，采用优化策略极小化整个递推格式频散关系，建立关于时间差分格式中差分系数的L2范数目标函数，利用线性优化算法对时间差分系数进行求解；
(6)采用混合吸收边界条件对边界反射能量进行吸收，按照时间高阶和空间隐式差分格式离散波动方程并进行递推，得到任意时刻的波场及整个地震记录；
(7)记录波场快照，输出地震记录并结束；
步骤(2)中，求解二阶时间导数时采用中心差分，形式如下：
其中，为压力场，τ为时间步长，t为时间；为了进一步提高时间离散精度，将菱形差分算子引入到公式(1)中，得出如下的时间高阶离散公式：
其中，h为空间采样步长，v为速度，am,n为有限差分系数；N为菱形算子中的差分算子长度；n和m为菱形算子中的坐标位置序列；
步骤(3)中，采用如下空间隐式离散格式对二阶空间导数进行近似，形式如下：
其中，a′0、a′m和b为隐式差分系数，当b＝0时，公式(3)退化为显式差分格式；h为空间步长，x为方向坐标，M为差分算子长度；公式(3)对应的空间频散关系如下：
其中，β＝kxh，kx为水平方向波数；极小化公式(4)对应的空间频散关系，目标函数为关于差分系数的线性凸函数，最优差分系数通过求解如下公式得到：
其中，n＝1,2,L,M+1，d＝[d1,d2,L,dM,dM+1]T＝[a′1,a′2,L,a′M-1,a′M,b]T，q为介于0和1之间的自然数，控制着积分范围，
步骤(4)中，在求解时间导数时，采用公式(2)计算；求解空间导数时，采用公式(3)计算，将公式(2)和公式(3)带入到声波方程并进行平面波理论分析，过程如下：
均匀模型中，声波方程压力P满足如下平面波方程，
其中，ω是角频率，是虚数单位，kz为垂直方向波数；m、n和l为位置和时间序列；
将公式(7)带入公式(2)并进行化简得到如下公式：
同样，对空间导数公式进行分析，得到如下频散关系：
将上述(8)、(9)和(10)带入声波方程并简化和重组得到如下时空域频散关系，
其中，ω为角频率，kx和kz分别为沿x和z方向的波数，r为库朗数，r＝vτ/h；
步骤(5)中，首先，差分系数a′m和b由公式(5)计算得到，在计算差分系数am,0和am,n时保持固定，此时优化差分系数am,0和am,n是线性优化问题，直接线性求解，定义如下目标函数：
其中，d＝[d1,L,dN+(N-1)N/2]T＝[a1,0,L,aN,0,a1,1,L,a1,N-1,L,aN-1,1]T，表示差分系数；θ为传播角度，q为介于0和1间的实数，波数kx和kz与波数k之间满足如下关系：
kx＝kcos(θ),kz＝ksin(θ), (13)ψ(β,θ,r)及定义如下：
目标函数是关于差分系数am,0和am,n的线性函数，差分系数am,0和am,n通过求解如下公式得到，
2.根据权利要求1所述的基于声波方程的线性优化隐式时空域有限差分数值模拟方法，其特征在于，步骤(6)中，整个计算区域被划分为内部区域、边界区域和过渡区域，具体的方法如下：
(6-1)利用公式(2)和公式(3)离散时间和空间导数，采用公式(5)和公式(16)计算差分系数，进行双程波方程求解，求解声波方程，计算所有区域处的双程波波场；
(6-2)求解单程波方程，计算边界区域和过渡区域处的单程波波场，计算公式及方法如下：
针对二阶声波方程，以左边界和左上角为例，左边界单程波方程为：
其中，x和z代表坐标；
左上角单程波方程为：
其他边界处及角点处单程波可以类似得到，并进行对应求解；采用上述吸收边界条件策略计算各位置处任意时刻波场值；
(6-3)计算最终波场：
边界处波场为单程波场值，内部区域处波场为双程波场值，中间过渡区域处波场通过如下线性组合得到：
P＝BiPtwo+(1-Bi)Pone, (19)
其中，Ptwo为双程波场，Pone为单程波场，Bi为权重系数，Bi＝(i-1)/Nb，i＝2,L,Nb，Nb为过渡区域吸收边界层数，过渡区域从内到外，Bi由(Nb-1)/Nb变化到1/Nb；如果Nb＝1，则此时混合吸收边界条件退化为单程波吸收边界条件。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711142912.2A CN107976710B (zh) | 2017-11-17 | 2017-11-17 | 一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711142912.2A CN107976710B (zh) | 2017-11-17 | 2017-11-17 | 一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107976710A CN107976710A (zh) | 2018-05-01 |
| CN107976710B true CN107976710B (zh) | 2019-05-28 |
Family
ID=62010247
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201711142912.2A Expired - Fee Related CN107976710B (zh) | 2017-11-17 | 2017-11-17 | 一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107976710B (zh) |
Families Citing this family (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110579796A (zh) * | 2018-06-08 | 2019-12-17 | 中国科学院地质与地球物理研究所 | 拓展有限差分稳定性条件的波场模拟方法、装置及设备 |
| CN110858002B (zh) * | 2018-08-23 | 2022-07-15 | 中国科学院地质与地球物理研究所 | 拓展显式差分稳定性条件的波场模拟方法、装置及设备 |
| CN109709602B (zh) * | 2018-11-22 | 2021-01-29 | 中国石油天然气股份有限公司 | 一种远探测声波偏移成像方法、装置及系统 |
| CN109725351A (zh) * | 2018-12-18 | 2019-05-07 | 中国石油天然气集团有限公司 | 一种3d弹性波混合吸收边界条件的确定方法、装置及系统 |
| CN109782338B (zh) * | 2019-01-03 | 2020-06-30 | 深圳信息职业技术学院 | 一种频率域地震正演模拟的高精度差分数值法 |
| CN112285772B (zh) * | 2020-10-07 | 2023-05-19 | 长安大学 | 有限差分数值模拟方法、系统、介质、计算机设备及应用 |
| CN114611349B (zh) * | 2022-03-02 | 2024-09-27 | 中国科学院声学研究所 | 一种基于隐式离散模型的频率域声波仿真方法及系统 |
| CN115034103B (zh) * | 2022-05-11 | 2024-08-16 | 中国石油大学(华东) | 一种基于优化差分系数的正演模拟方法及其设备 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8456953B2 (en) * | 2010-01-18 | 2013-06-04 | Westerngeco L.L.C. | Wave equation illumination |
| CN103630933A (zh) * | 2013-12-09 | 2014-03-12 | 中国石油天然气集团公司 | 基于非线性优化的时空域交错网格有限差分方法和装置 |
| CN104597488A (zh) * | 2015-01-21 | 2015-05-06 | 中国石油天然气集团公司 | 非等边长网格波动方程有限差分模板优化设计方法 |
| CN107179549A (zh) * | 2017-07-11 | 2017-09-19 | 中海石油(中国)有限公司 | 一种时间域声波方程显式有限差分地震响应模拟方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107942375B (zh) * | 2017-11-17 | 2019-04-30 | 河海大学 | 一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法 |
-
2017
- 2017-11-17 CN CN201711142912.2A patent/CN107976710B/zh not_active Expired - Fee Related
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8456953B2 (en) * | 2010-01-18 | 2013-06-04 | Westerngeco L.L.C. | Wave equation illumination |
| CN103630933A (zh) * | 2013-12-09 | 2014-03-12 | 中国石油天然气集团公司 | 基于非线性优化的时空域交错网格有限差分方法和装置 |
| CN104597488A (zh) * | 2015-01-21 | 2015-05-06 | 中国石油天然气集团公司 | 非等边长网格波动方程有限差分模板优化设计方法 |
| CN107179549A (zh) * | 2017-07-11 | 2017-09-19 | 中海石油(中国)有限公司 | 一种时间域声波方程显式有限差分地震响应模拟方法 |
Non-Patent Citations (2)
| Title |
|---|
| 声波方程隐式时空域有限差分正演模拟方法;王恩江，等;《2017中国地球科学联合学术年会》;20171031;732-733 |
| 适于声波方程数值模拟的时间-空间域有限差分算子改进线性算法;梁文全，等;《地震工程学报》;20161031;815-821 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107976710A (zh) | 2018-05-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107976710B (zh) | 一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法 | |
| CN107942375B (zh) | 一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法 | |
| CN105044771B (zh) | 基于有限差分法的三维tti双相介质地震波场数值模拟方法 | |
| CN110542924B (zh) | 一种高精度的纵横波阻抗反演方法 | |
| CN110058303B (zh) | 声波各向异性逆时偏移混合方法 | |
| CN108828668B (zh) | 一种叠前时间偏移数据处理方法及装置 | |
| CN105652320B (zh) | 逆时偏移成像方法和装置 | |
| CN104570106A (zh) | 一种近地表层析速度分析方法 | |
| CN108181653A (zh) | 针对vti介质逆时偏移方法、设备及介质 | |
| CN110531410A (zh) | 一种基于直达波场的最小二乘逆时偏移梯度预条件方法 | |
| Mu et al. | Attenuation compensation and anisotropy correction in reverse time migration for attenuating tilted transversely isotropic media | |
| CN102866426A (zh) | 一种利用avo大角度道集分析岩体油气信息的方法 | |
| CN107966729B (zh) | 一种三维tti介质射线追踪方法及系统 | |
| CN109521470B (zh) | 分析地质构造对地震反演裂缝密度影响的方法 | |
| CN115327626B (zh) | 一种用于三维ati介质的弹性波场矢量分解方法及系统 | |
| Xu et al. | Time-space-domain temporal high-order staggered-grid finite-difference schemes by combining orthogonality and pyramid stencils for 3D elastic-wave propagation | |
| CN109738944B (zh) | 基于广角反射的地震采集参数确定方法及装置 | |
| CN114325832B (zh) | 一种裂缝参数和弹性参数同步反演方法及系统 | |
| CN112147700A (zh) | 速度异常区的低频模型构建方法及系统 | |
| CN109725346A (zh) | 一种时间-空间高阶vti矩形网格有限差分方法及装置 | |
| Xu et al. | SUB-TRIANGLE SHOOTING RAY-TRACING IN COMPLEX | |
| CN110261896B (zh) | 稳定的各向异性ti介质正演模拟方法 | |
| CN102508291B (zh) | 横向变速小尺度体反射系数公式应用方法 | |
| CN110568497A (zh) | 一种复杂介质条件下地震初至波旅行时的精确求解方法 | |
| CN103852786B (zh) | 一种应用于陆上地震数据的逆时偏移成像的方法及系统 |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190528 |