CN109782338B - 一种频率域地震正演模拟的高精度差分数值法 - Google Patents

一种频率域地震正演模拟的高精度差分数值法 Download PDF

Info

Publication number
CN109782338B
CN109782338B CN201910003548.4A CN201910003548A CN109782338B CN 109782338 B CN109782338 B CN 109782338B CN 201910003548 A CN201910003548 A CN 201910003548A CN 109782338 B CN109782338 B CN 109782338B
Authority
CN
China
Prior art keywords
order
dispersion
frequency domain
equation
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
Application number
CN201910003548.4A
Other languages
English (en)
Other versions
CN109782338A (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.)
Shenzhen Institute of Information Technology
Original Assignee
Shenzhen Institute of Information Technology
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 Shenzhen Institute of Information Technology filed Critical Shenzhen Institute of Information Technology
Priority to CN201910003548.4A priority Critical patent/CN109782338B/zh
Publication of CN109782338A publication Critical patent/CN109782338A/zh
Application granted granted Critical
Publication of CN109782338B publication Critical patent/CN109782338B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Geophysics And Detection Of Objects (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提出一种频率域地震正演模拟的高精度差分数值法，本方法对于频率域波动方程高阶项的离散，先在25点模板的水平（垂直）方向分别进行单倍和双倍步长的二阶离散，后在垂直（水平）方向利用对称性和泰勒展开式消除二阶项，得到两种四阶精度的逼近。接着，把这两种四阶精度的逼近与传统的四阶中心差分格式进行加权组合，权系数和为1。对于低阶项的离散，先利用25点模板四个角处的16个网格点进行离散，得到四阶精度逼近，接着再与水平和竖直方向的四阶中心离散进行加权组合，权系数的和为1。最后对得到的差分方案进行频散分析，并基于频散极小化原理得到各个权系数的值。本方法达到了四阶收敛，有效减少了数值频散，提升正演模拟的精度与效率。

Description

一种频率域地震正演模拟的高精度差分数值法
技术领域
本发明涉及地震波模拟领域，更具体地，涉及一种频率域地震正演模拟的高精度差分数值法。
背景技术
地震波正演模拟是研究地震波传播规律的重要工具，也是反演地下层结构图像的基础，被广泛地应用于地球物理、海洋技术、油气勘探等科学技术领域。近几十年来，地震波正演模拟技术一直是地震勘探领域的重要研究内容。地震波正演模拟通常在时间域或频率域进行。相对于时间域模拟，频率域模拟更具优势。频率域波动方程正演过程只需要几个频点的信息，大大提高了计算效率。此外，在频率域可以更加灵活地模拟地震波衰减效应。频率域模拟进一步分为有限差分模拟和有限元模拟等两类主要方法。相对于有限元法，有限差分法计算复杂度较低，易于执行，故在地震勘探领域有着广泛的应用。当前，地震波频率域差分模拟方法主要是基于混合交错网格设计差分方案，并通过极小化数值频散机制得到混合交错网格上的加权组合系数。相对于传统有限差分方案，交错网格方案可有效减小数值频散，提升数值模拟精度。
当前，频率域地震正演模拟的混合交错网格差分法主要有9点方案、13点、17点方案、25点方案。通常情况下，使用较多网格点数的差分方案，其数值频散也较小。其中，25点差分方案展现了优良的性能。但是，当前25点差分方案尚存在以下几个问题。首先，多数方案的差分格式是2阶收敛，收敛精度不高。其次，有一部分声称的高阶方案本质上是伪高阶方案，即方程的高阶项采用的是高阶离散，但是低阶项采用的却是二阶离散，因此总体上还是二阶收敛精度。最后，有些基于旋转坐标轴的差分方案虽是高阶，但添加边界条件后，如PML(Perfectly Matched Layer)，相应的差分格式与原方程并非点态相容(pointwiseconsistent)，一定程度上影响了模拟精度。
发明内容
本发明为克服上述现有技术所述的至少一种缺陷(不足)，提供一种频率域地震正演模拟的高精度差分数值法，以提升地震正演数值精度和效率。
为解决上述技术问题，本发明的技术方案如下：
一种频率域地震正演模拟的高精度差分数值法，包括以下步骤：
S1、对二维标量频率域波动方程高阶项离散，首先在25点模板的水平(垂直)方向分别进行单倍和双倍步长的二阶离散，然后在垂直(水平)方向利用对称性和泰勒展开式消除二阶项，得到两种四阶精度的逼近；
S2、把这两种四阶精度的逼近与经典的四阶中心差分格式进行加权组合；其中，为了确保收敛性，权系数的和为1；
S3、对于低阶项的离散，先利用25点模板四个角处的16个网格点进行离散，得到四阶精度逼近，接着再与水平和竖直方向的四阶中心离散进行加权组合，权系数的和为1；
S4、对得到的差分方案进行频散分析，并基于频散极小化原理得到各个权系数的值
优选的，已知二维标量频率域波动方程为：
Figure GDA0002498536500000021
这里，未知量u表示压力场，
Figure GDA0002498536500000022
为波数，其中f和v分别表示波的频率和波的传播速度；
为了在有限的计算区域内进行模拟仿真，对方程(1)进行添加PML边界层，得到
Figure GDA0002498536500000023
其中
Figure GDA0002498536500000024
C＝exey，
Figure GDA0002498536500000025
ω为角频率，σx定义如下：
Figure GDA0002498536500000026
其中f0为震源主频，a0＝1.79为一经验常数，LPML为PML的厚度，lx为到PML和内部区域间界面的距离；σy的定义与σx类似；其中，在内部区域Ω0中，ex＝ey＝1，在PML区域Ω1中ex＝1，在PML区域Ω2中ey＝1；
记um,n为函数u在网格点(xm,yn)处的离散值，定义符号
Figure GDA0002498536500000031
如下：
Figure GDA0002498536500000032
其中j∈{-2,-1,1,2},t＝1,2,h为离散步长；定义符号：
Figure GDA0002498536500000033
Figure GDA0002498536500000034
继续定义符号Φxu为：
Figure GDA0002498536500000035
其中α1+α2+α3＝1Φxu即为高阶项
Figure GDA0002498536500000036
的离散逼近；
类似地，分别定义符号
Figure GDA0002498536500000037
Φx如下：
Figure GDA0002498536500000038
Figure GDA0002498536500000039
Figure GDA00024985365000000310
Figure GDA00024985365000000311
Φyu即为高阶项
Figure GDA00024985365000000312
的离散逼近；则频率域波动方程(2)的高阶项离散逼近为：
-Φxu-Φyu (3)
定义符号
Figure GDA00024985365000000313
Ι1，Ι2如下，其中j＝{-2,-1,0,1,2}；
Figure GDA0002498536500000041
Figure GDA0002498536500000042
I1(Ck2u)＝(Ck2u)m,n,
Figure GDA0002498536500000043
Figure GDA0002498536500000044
则频率域波动方程(2)的低阶项离散逼近为
I(Ck2u):＝β1I1(Ck2u)+β2Ι2(Ck2u)+β3Ι3(Ck2u), (4)
其中β1+β2+β3＝1，综合(3)式与(4)式，得到频率域波动方程(2)的4阶差分逼近；
-Φxu-Φyu-I(Ck2u)＝0 (5)
在PML内部区域，25点四阶差分方案(5)简化为
Figure GDA0002498536500000045
其中，Um+j,n+l为未知量u在(xm+j,yn+l)处的离散值，j,l＝-2,-1,0,1,2，
Figure GDA0002498536500000046
Figure GDA0002498536500000047
令U(x,y):＝e-ik(xcosθ+ysinθ)，则U(x,y)是频率域波动方程的面波解，其中θ为波传播方向与y轴的夹角；记λ为波长，G为单个波长内的网格点数，即G＝λ/h；又由于λ＝v/f，故有kh＝2π/G；把离散面波解
Figure GDA0002498536500000048
带入差分格式(6)，并利用欧拉公式e-ix＝cosx+isinx式进行变换，得到频散方程如下4T1P2Q2+4T2(P1Q2+P2Q1)+2T3(P2+Q2)+4T4P1Q1+2T5(P1+Q1)+T4＝0, (7)
其中
Figure GDA0002498536500000051
Figure GDA0002498536500000052
把Tj(j＝1,2,...,6)的值带入(7)式，并用数值波数kN替换k，得到
Figure GDA0002498536500000053
其中
Figure GDA0002498536500000054
R＝a1M+a2N+a3O,
M＝[(P2+Q2)-16(P1+Q1)+30],
N＝[2P2Q2-4(P1Q2+P2Q1)-(P2+Q2)+4(P1+Q1)],
O＝[4(P1Q2+P2Q1)-4(P2+Q2)-32P1Q1+16(P1+Q1)]
联合kh＝2π/G和(8)式得到
Figure GDA0002498536500000055
记
Figure GDA0002498536500000056
则求参数α1,α2,α3,β1,β2,β3的值转化为求最优化问题
Figure GDA0002498536500000057
其中，
Figure GDA0002498536500000058
由对称性得到，而G≥2由Nyqusit采样定理得到；
令
Figure GDA0002498536500000059
并带入α1＝1-α2-α3,β1＝1-β2-β3,得到方程(11)
Figure GDA00024985365000000510
其中P1,Q1,P2,Q2,O,M,N是关于θ,G的函数；分别在θ和G的取值范围内采样l＝10和r＝100个点，即：
Figure GDA0002498536500000061
把(12)中的l*r＝1000个点带入方程(11)得到超定线性系统
Figure GDA0002498536500000062
其中，
Figure GDA0002498536500000063
Figure GDA0002498536500000064
Figure GDA0002498536500000065
Figure GDA0002498536500000066
Figure GDA0002498536500000067
利用最小二乘法求解线性系统(13)，即
Figure GDA0002498536500000071
得到α2＝0.0266,α3＝0.1923,β1＝0.0970,β2＝0.1902,则α1＝1-α2-α3＝0.7811,
β1＝1-β2-β3＝0.7128
把参数α1,α2,α3,β1,β2,β3的值带入(5)式，得到最终的频率波动方程正演模拟的4阶25点差分方案。
与现有技术相比，本发明技术方案的有益效果是：
本发明提出了一种频率域地震波正演模拟的25点四阶高精度差分数值法，以提升地震正演数值精度和效率。本发明的有益效果如下：
(1)本发明构造了一种新的、不同于当前25点方案的频率域波动方程正演方法。对于方程的高阶项和低阶项，分别设计多种离散逼近格式，然后对它们进行有机加权组合，并基于频散极小化机制，利用最小二乘法求得加权系数。该方法易于实施，且有效地减少了数值频散，提升了正演模拟的计算效率。
(2)本技术方案的差分数值法是真正的高(四)阶精度差分方案，有效地提升了数值模拟精度。对于方程高阶项的离散，先在25点模板的一个方向(水平或垂直方向)分别进行单倍和双倍步长的二阶离散，然后再在另一个方向利用对称性和泰勒展开式消除二阶项，得到两种四阶精度的逼近。接着，把这两种四阶精度的逼近与传统的四阶中心差分格式进行加权组合。对于低阶项的离散，先利用25点模板四个角处的16个网格点进行离散，得到四阶精度逼近，接着再与水平和竖直方向的四阶中心离散进行加权组合。在理论和数值上，本方案都被证明是4阶收敛的。
(3)本技术方案的差分格式在方程添加了PML边界层后，仍与带PML的方程点态相容，这与旋转差分等格式不同。点态相容性是离散方案性质的一个重要指标，其确保了数值模拟的收敛性和精度。
(4)本技术方案的差分数值法实用性好，对光滑性要求不高。目前，一些高阶方案对解的光滑性和方程右端项的要求苛刻，因而限制了其实际应用。本技术方案没有这方面的限制，特别是对于右端项，因此本方案实用性更高。
附图说明
图1为加密算法执行流程图PML示意图；
图2为
Figure GDA0002498536500000081
离散示意图；
图3为
Figure GDA0002498536500000082
的离散示意图；
图4为低阶项Ck2u的离散示意图；
图5为本发明差分格式的标准化相速度曲线示意图；
图6为基于本差分方案的频率域地震波正演模拟(均匀介质)示意图；6(a)为f＝5，6(b)为f＝20；
图7为盐丘体速度场模型示意图；
图8为基于本差分方案的频率域地震波正演模拟(盐丘体模型，非均匀介质)示意图；
图9为盐丘体模型的波场快照示意图。9(a)为t＝50ms，9(b)为t＝100ms。
具体实施方式
对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。下面结合附图和实施例对本发明的技术方案做进一步的说明。
本发明提出一种频率域地震波正演模拟的25点四阶高精度差分数值法，以进一步提升地震正演数值精度和效率。该方案不仅真正具有4阶收敛精度，同时，还与添加PML边界层后的方程保持点态相容，且对光滑性要求不高。在本方法中，对于频率域波动方程高阶项的离散，首先在25点模板的水平(垂直)方向分别进行单倍和双倍步长的二阶离散，然后在垂直(水平)方向利用对称性和泰勒展开式消除二阶项，得到两种四阶精度的逼近。接着，把这两种四阶精度的逼近与传统的四阶中心差分格式进行加权组合。其中，为了确保收敛性，权系数的和为1。对于低阶项的离散，先利用25点模板四个角处的16个网格点进行离散，得到四阶精度逼近，接着再与水平和竖直方向的四阶中心离散进行加权组合。同理，权系数的和也为1。最后，对得到的差分方案进行频散分析，并基于频散极小化原理得到各个权系数的值。仿真实验显示，本数值方法达到了四阶收敛，且有效地减少了数值频散，提升了正演模拟的精度与效率。
本发明是为了进一步提升频率域地震正演模拟的精度和效率，为地球物理、油气勘探等科学技术领域提供重要的应用支持。
其具体的实现过程为：已知二维标量频率域波动方程为，
Figure GDA0002498536500000091
这里，未知量u表示压力场，
Figure GDA0002498536500000092
为波数，其中f和v分别表示波的频率和波的传播速度。
为了在有限的计算区域内进行模拟仿真，对方程(1)进行添加PML边界层，得到
Figure GDA0002498536500000093
其中
Figure GDA0002498536500000094
C＝exey，
Figure GDA0002498536500000095
ω为角频率，σx定义如下：
Figure GDA0002498536500000096
其中f0为震源主频，a0＝1.79为一经验常数，LPML为PML的厚度，lx为到PML和内部区域间界面的距离。σy的定义与σx类似。PML的示意图见图1。其中，在内部区域Ω0中，ex＝ey＝1，在PML区域Ω1中ex＝1，在PML区域Ω2中ey＝1。
记um,n为函数u在网格点(xm,yn)处的离散值，定义符号
Figure GDA0002498536500000097
如下
Figure GDA0002498536500000098
其中j∈{-2,-1,1,2},t＝1,2,h为离散步长。定义符号
Figure GDA0002498536500000099
Figure GDA0002498536500000101
继续定义符号Φxu为：
Figure GDA0002498536500000102
其中α1+α2+α3＝1.Φxu即为高阶项
Figure GDA0002498536500000103
的离散逼近，离散示意图见图2。
类似地，分别定义符号
Figure GDA0002498536500000104
Φyu如下：
Figure GDA0002498536500000105
Figure GDA0002498536500000106
Figure GDA0002498536500000107
Figure GDA0002498536500000108
Φyu即为高阶项
Figure GDA0002498536500000109
的离散逼近，离散示意图见图3。则频率域波动方程(2)的高阶项离散逼近为：
-Φxu-Φyu (3)
定义符号
Figure GDA00024985365000001010
Ι1，Ι2，Ι3如下(j＝{-2,-1,0,1,2})
Figure GDA00024985365000001011
Figure GDA00024985365000001012
Figure GDA00024985365000001013
Figure GDA00024985365000001014
则频率域波动方程(2)的低阶项离散逼近为
I(Ck2u):＝β1I1(Ck2u)+β2Ι2(Ck2u)+β3Ι3(Ck2u), (4)
其中β1+β2+β3＝1.低阶项Ck2u的离散示意图见图4，综合(3)式与(4)式，得到频率域波动方程(2)的4阶差分逼近。
-Φxu-Φyu-I(Ck2u)＝0 (5)
在PML内部区域，25点四阶差分方案(5)可简化为
Figure GDA0002498536500000111
其中，Um+j,n+l(j,l＝-2,-1,0,1,2)为未知量u在(xm+j,yn+l)处的离散值，
Figure GDA0002498536500000112
Figure GDA0002498536500000113
令U(x,y):＝e-ik(xcosθ+ysinθ)，则U(x,y)是频率域波动方程的面波解，其中θ为波传播方向与y轴的夹角。记λ为波长，G为单个波长内的网格点数，即G＝λ/h.又由于λ＝v/f，故有kh＝2π/G。把离散面波解
Figure GDA0002498536500000114
带入差分格式(6)，并利用欧拉公式e-ix＝cosx+isinx式进行变换，得到频散方程如下4T1P2Q2+4T2(P1Q2+P2Q1)+2T3(P2+Q2)+4T4P1Q1+2T5(P1+Q1)+T4＝0, (7)
其中
Figure GDA0002498536500000115
Figure GDA0002498536500000116
把Tj(j＝1,2,...,6)的值带入(7)式，并用数值波数kN替换k，得到
Figure GDA0002498536500000117
其中
Figure GDA0002498536500000121
R＝a1M+a2N+a3O,
M＝[(P2+Q2)-16(P1+Q1)+30],
N＝[2P2Q2-4(P1Q2+P2Q1)-(P2+Q2)+4(P1+Q1)],
O＝[4(P1Q2+P2Q1)-4(P2+Q2)-32P1Q1+16(P1+Q1)].
联合kh＝2π/G和(8)式得到
Figure GDA0002498536500000122
记
Figure GDA0002498536500000123
则求参数α1,α2,α3,β1,β2,β3的值转化为求最优化问题
Figure GDA0002498536500000124
其中，
Figure GDA0002498536500000125
由对称性得到，而G≥2由Nyqusit采样定理得到。
令
Figure GDA0002498536500000126
并带入α1＝1-α2-α3,β1＝1-β2-β3,得到方程(11)
Figure GDA0002498536500000127
其中P1,Q1,P2,Q2,O,M,N是关于θ,G的函数。分别在θ和G的取值范围内采样l＝10和r＝100个点，即
Figure GDA0002498536500000128
把(12)中的l*r＝1000个点带入方程(11)得到超定线性系统
Figure GDA0002498536500000131
其中，
Figure GDA0002498536500000132
Figure GDA0002498536500000133
Figure GDA0002498536500000134
Figure GDA0002498536500000135
Figure GDA0002498536500000136
利用最小二乘法求解线性系统(13)，即
Figure GDA0002498536500000141
得到α2＝0.0266,α3＝0.1923,β1＝0.0970,β2＝0.1902,则α1＝1-α2-α3＝0.7811,
β1＝1-β2-β3＝0.7128.
把参数α1,α2,α3,β1,β2,β3的值带入(5)式，得到最终的频率波动方程正演模拟的4阶25点差分方案。
仿真实验
(1)数值频散分析。首先给出本发明25点4阶差分格式的标准相速度曲线，即以单位波长内网格点数的倒数1/G为横坐标，以数值波数和精确波数的比值kN/k为纵坐标的曲线。图5显示了4条标准化相速度曲线，分别对应θ＝0°,15°,30°,45°；可以看出，标准化相速度曲线接近直线y＝1，且即使每个波长内仅取2个网格点(达到Nyqusit采样极限)时，kN/k的值仍接近于1，表明相应的数值频散较小，有效地提升了数值模拟精度。
(2)收敛阶测试。考虑求解右端项非零的频率域波动方程(14)
Figure GDA0002498536500000142
其中，右端项g(x,y)＝(2n2-1)sin(nkx)sin(nky)，则方程有精确解
Figure GDA0002498536500000143
利用本方案提出的差分数值法求解(14)，表1列出了模拟误差(在|·|∞的意义下)和数值收敛阶。可以看出，计算的数值收敛阶接近于4，表明本差分方案是4阶精度收敛的。
表1本差分方案的数值收敛阶
Figure GDA0002498536500000144
Figure GDA0002498536500000151
(3)频率域正演模拟测试。设置方程(2)的右端项为点震源，下面针对均匀传播介质和非均匀传播介质两种情况进行正演模拟。对于均匀介质情形，计算区域为2160m×2160m，波速v＝2000，频率分别取f＝5,20，对仿真结果的实部进行可视化，如图6所示。
对于均匀介质情形，采用盐丘体速度场模型进行测试。盐丘体模型模拟了海底地质结构，考虑矩形区域2160m×4720m，相应的速度场如图7所示。地震波在盐丘体中最小传播速度为1524m/s，最大传播速度为4480m/s。把点震源放置在(1080m,2360m)处，取f＝20，基于本差分方案对频率域波动方程进行正演模拟，模拟结果的实部可视化如图8所示。可见，波的传播从震源处开始，穿过整个盐丘体区域。接着，生成盐丘体模型对应的波长快照。波长快照给出了具体时刻波前(wavefront)的位置，其形状取决于速度的分布。基于本差分方案，模拟60个单频率波(f＝1,2,...,59,60)的传播情况，在进行波场快照时，取200个采样点，且时间间隔为0.005s。图9显示了t＝50ms(毫秒)和t＝100ms时的快照。频率域正演模拟测试显示了本差分方案的优良性能。
显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (4)

1.一种频率域地震正演模拟的高精度差分数值法，其特征在于，包括以下步骤：
S1、对二维标量频率域波动方程高阶项离散，首先在25点模板的水平或垂直方向分别进行单倍和双倍步长的二阶离散，然后在垂直或水平方向利用对称性和泰勒展开式消除二阶项，得到两种四阶精度的逼近；
S2、把这两种四阶精度的逼近与经典的四阶中心差分格式进行加权组合；其中，权系数的和为1；
S3、对于低阶项的离散，先利用25点模板四个角处的16个网格点进行离散，得到四阶精度逼近，接着再与水平和竖直方向的四阶中心离散进行加权组合，权系数的和为1；
S4、对得到的差分方案进行频散分析，并基于频散极小化原理得到各个权系数的值。
2.根据权利要求1所述的一种频率域地震正演模拟的高精度差分数值法，其特征在于，已知二维标量频率域波动方程为：
Figure FDA0002498536490000011
这里，未知量u表示压力场，
Figure FDA0002498536490000012
为波数，其中f和v分别表示波的频率和波的传播速度；
为了在有限的计算区域内进行模拟仿真，对方程(1)进行添加PML边界层，得到
Figure FDA0002498536490000013
其中
Figure FDA0002498536490000014
C＝exey，
Figure FDA0002498536490000015
ω为角频率，σx定义如下：
Figure FDA0002498536490000016
其中f0为震源主频，a0为一经验常数，LPML为PML的厚度，lx为到PML和内部区域间界面的距离；σy的定义与σx类似；
记um,n为函数u在网格点(xm,yn)处的离散值，定义符号
Figure FDA0002498536490000021
如下：
Figure FDA0002498536490000022
其中j∈{-2,-1,1,2},t＝1,2,h为离散步长；定义符号：
Figure FDA0002498536490000023
Figure FDA0002498536490000024
继续定义符号Φxu为：
Figure FDA0002498536490000025
其中α1+α2+α3＝1，Φxu即为高阶项
Figure FDA0002498536490000026
的离散逼近；
类似地，分别定义符号
Figure FDA0002498536490000027
Figure FDA0002498536490000028
Φyu如下：
Figure FDA0002498536490000029
Figure FDA00024985364900000210
Figure FDA00024985364900000211
Figure FDA00024985364900000212
Φyu即为高阶项
Figure FDA00024985364900000213
的离散逼近；则频率域波动方程(2)的高阶项离散逼近为：
-Φxu-Φyu (3)
定义符号
Figure FDA00024985364900000214
Ι1，Ι2，Ι3如下，其中j＝{-2,-1,0,1,2}；
Figure FDA0002498536490000031
Figure FDA0002498536490000032
I1(Ck2u)＝(Ck2u)m,n,
Figure FDA0002498536490000033
Figure FDA0002498536490000034
则频率域波动方程(2)的低阶项离散逼近为
I(Ck2u):＝β1I1(Ck2u)+β2Ι2(Ck2u)+β3Ι3(Ck2u), (4)
其中β1+β2+β3＝1，综合(3)式与(4)式，得到频率域波动方程(2)的4阶差分逼近；
-Φxu-Φyu-I(Ck2u)＝0 (5)
在PML内部区域，25点四阶差分方案(5)简化为
Figure FDA0002498536490000035
其中，Um+j,n+l为未知量u在(xm+j,yn+l)处的离散值，j,l＝-2,-1,0,1,2，系数T1至T6为：
Figure FDA0002498536490000036
Figure FDA0002498536490000037
令U(x,y):＝e-ik(xcosθ+ysinθ)，则U(x,y)是频率域波动方程的面波解，其中θ为波传播方向与y轴的夹角；记λ为波长，G为单个波长内的网格点数，即G＝λ/h；又由于λ＝v/f，故有kh＝2π/G；
把离散面波解
Figure FDA0002498536490000038
带入差分格式(6)，并利用欧拉公式e-ix＝cosx+isinx式进行变换，得到频散方程如下：
4T1P2Q2+4T2(P1Q2+P2Q1)+2T3(P2+Q2)+4T4P1Q1+2T5(P1+Q1)+T4＝0, (7)
其中
Figure FDA0002498536490000041
Figure FDA0002498536490000042
把Tj的值带入(7)式，j＝1,2,...,6；并用数值波数kN替换k，得到
Figure FDA0002498536490000043
其中
Figure FDA0002498536490000044
R＝a1M+a2N+a3O,
M＝[(P2+Q2)-16(P1+Q1)+30],
N＝[2P2Q2-4(P1Q2+P2Q1)-(P2+Q2)+4(P1+Q1)],
O＝[4(P1Q2+P2Q1)-4(P2+Q2)-32P1Q1+16(P1+Q1)]
联合kh＝2π/G和(8)式得到
Figure FDA0002498536490000045
记
Figure FDA0002498536490000046
则求参数α1,α2,α3,β1,β2,β3的值转化为求最优化问题
Figure FDA0002498536490000047
其中，
Figure FDA00024985364900000410
由对称性得到，而G≥2由Nyqusit采样定理得到；
令
Figure FDA0002498536490000048
并带入α1＝1-α2-α3,β1＝1-β2-β3,得到方程(11)
Figure FDA0002498536490000049
其中P1,Q1,P2,Q2,O,M,N是关于θ,G的函数；分别在θ和G的取值范围内采样l＝10和r＝100个点，即
Figure FDA0002498536490000051
把(12)中的l*r＝1000个点带入方程(11)得到超定线性系统
Figure FDA0002498536490000052
其中，
Figure FDA0002498536490000053
Figure FDA0002498536490000054
Figure FDA0002498536490000055
Figure FDA0002498536490000056
Figure FDA0002498536490000057
利用最小二乘法求解线性系统(13)，即
Figure FDA0002498536490000061
得到α2＝0.0266,α3＝0.1923,β1＝0.0970,β2＝0.1902,则α1＝1-α2-α3＝0.7811,
β1＝1-β2-β3＝0.7128
把参数α1,α2,α3,β1,β2,β3的值带入(5)式，得到最终的频率波动方程正演模拟的4阶25点差分方案。
3.根据权利要求2所述的一种频率域地震正演模拟的高精度差分数值法，其特征在于，经验常数a0＝1.79。
4.根据权利要求2或3所述的一种频率域地震正演模拟的高精度差分数值法，其特征在于，在PML内部区域Ω0中，ex＝ey＝1，在PML区域Ω1中ex＝1，在PML区域Ω2中ey＝1。
CN201910003548.4A 2019-01-03 2019-01-03 一种频率域地震正演模拟的高精度差分数值法 Active CN109782338B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910003548.4A CN109782338B (zh) 2019-01-03 2019-01-03 一种频率域地震正演模拟的高精度差分数值法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910003548.4A CN109782338B (zh) 2019-01-03 2019-01-03 一种频率域地震正演模拟的高精度差分数值法

Publications (2)

Publication Number Publication Date
CN109782338A CN109782338A (zh) 2019-05-21
CN109782338B true CN109782338B (zh) 2020-06-30

Family

ID=66499782

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910003548.4A Active CN109782338B (zh) 2019-01-03 2019-01-03 一种频率域地震正演模拟的高精度差分数值法

Country Status (1)

Country Link
CN (1) CN109782338B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102183790A (zh) * 2011-02-12 2011-09-14 中国石油大学(华东) 基于时空双变网格的弹性波正演模拟技术
CN103901472A (zh) * 2014-03-31 2014-07-02 中国科学院地质与地球物理研究所 一种频率域正演方法及装置
CN107479092A (zh) * 2017-08-17 2017-12-15 电子科技大学 一种基于方向导数的频率域高阶声波方程正演模拟方法
CN107976710A (zh) * 2017-11-17 2018-05-01 河海大学 一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法
CN108108331A (zh) * 2017-12-13 2018-06-01 国家深海基地管理中心 一种基于拟空间域弹性波方程的有限差分计算方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CA2972028C (en) * 2015-02-13 2019-08-13 Exxonmobil Upstream Research Company Efficient and stable absorbing boundary condition in finite-difference calculations

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102183790A (zh) * 2011-02-12 2011-09-14 中国石油大学(华东) 基于时空双变网格的弹性波正演模拟技术
CN103901472A (zh) * 2014-03-31 2014-07-02 中国科学院地质与地球物理研究所 一种频率域正演方法及装置
CN107479092A (zh) * 2017-08-17 2017-12-15 电子科技大学 一种基于方向导数的频率域高阶声波方程正演模拟方法
CN107976710A (zh) * 2017-11-17 2018-05-01 河海大学 一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法
CN108108331A (zh) * 2017-12-13 2018-06-01 国家深海基地管理中心 一种基于拟空间域弹性波方程的有限差分计算方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
2D声波频率域数值模拟中几种有限差分方法的对比分析;马晓娜 等;《地球物理学进展》;20151231;第30卷(第2期);第878-888页 *
频率域波动方程的加权五点有限差分数值模拟;程东升 等;《中国科技信息》;20181031(第19期);第96-97,99页 *

Also Published As

Publication number Publication date
CN109782338A (zh) 2019-05-21

Similar Documents

Publication Publication Date Title
Helbig Foundations of Anisotropy for Exploration Seismics: Section I. Seismic Exploration
Pavlis et al. The generalized earthquake-location (GENLOC) package: An earthquake-location library
Bennett Array design by inverse methods
Ivannikova et al. Regional 3‐D modeling of ground electromagnetic field due to realistic geomagnetic disturbances
Gao et al. A new approach for extracting the amplitude spectrum of the seismic wavelet from the seismic traces
Jekeli Spectral methods in geodesy and geophysics
DiNapoli Fast field program for multilayered media
CN110058303A (zh) 声波各向异性逆时偏移混合方法
Jing et al. Dispersion-relation preserving stereo-modeling method beyond Nyquist frequency for acoustic wave equation
Salehipour et al. A higher order discontinuous Galerkin, global shallow water model: Global ocean tides and aquaplanet benchmarks
CN109782338B (zh) 一种频率域地震正演模拟的高精度差分数值法
CN113609700A (zh) 一种山区桥梁完全非平稳风场的模拟方法
Frederick et al. Numerical methods for multiscale inverse problems
Li et al. Decoupling approximation of P-and S-wave phase velocities in orthorhombic media
Yuen et al. Geophysical applications of multidimensional filtering with wavelets
CN113885101B (zh) 一种基于椭球模型构建重力梯度基准图方法
CN115508898A (zh) G-s变换的接地长导线源瞬变电磁快速正反演方法及系统
Mu et al. Efficient pure qP-wave simulation and reverse time migration imaging for vertical transverse isotropic (VTI) media
Tu et al. A wavenumber integration model of underwater acoustic propagation in arbitrary horizontally stratified media based on a spectral method
Barucq et al. Space-time Trefftz-discontinuous Galerkin approximation for elasto-acoustics
CN111506871A (zh) 基于频域特性构建海域垂线偏差模型的径向基函数格网法
Nie et al. Prestack Seismic Inversion via Nonconvex L1-2 Regularization
De Hoop et al. Semiclassical inverse spectral problem for seismic surface waves in isotropic media: part I. Love waves
Feng et al. A finite volume MHD code in spherical coordinates for background solar wind
CN113267830A (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