CN115932949B - 一种基于自适应系数频域有限差分的黏声波模拟方法 - Google Patents

一种基于自适应系数频域有限差分的黏声波模拟方法 Download PDF

Info

Publication number
CN115932949B
CN115932949B CN202211697099.6A CN202211697099A CN115932949B CN 115932949 B CN115932949 B CN 115932949B CN 202211697099 A CN202211697099 A CN 202211697099A CN 115932949 B CN115932949 B CN 115932949B
Authority
CN
China
Prior art keywords
frequency domain
finite difference
wave
equation
adaptive coefficient
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
CN202211697099.6A
Other languages
English (en)
Other versions
CN115932949A (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.)
Hohai University HHU
Original Assignee
Hohai University HHU
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 Hohai University HHU filed Critical Hohai University HHU
Priority to CN202211697099.6A priority Critical patent/CN115932949B/zh
Publication of CN115932949A publication Critical patent/CN115932949A/zh
Application granted granted Critical
Publication of CN115932949B publication Critical patent/CN115932949B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本发明公开了一种基于自适应系数频域有限差分(finite difference frequency domain,FDFD)的黏声波模拟方法：构建自适应系数FDFD通用格式；确定黏声波方程在预设参考频率下的波传播速度模型和品质因子模型；建立关于一个网格点内波长个数的自适应FDFD系数查找表，验证查找表的有效性；在针对黏声波波场模拟的自适应系数FDFD通用格式中加入完美匹配层吸收边界条件消除人工边界反射；通过MPI多进程并行和PARDISO直接求解器求解自适应系数FDFD方法生成的多个大型稀疏线性方程组得到不同频率下的黏声波波场。数值模拟实验表明，对于不同的空间采样间隔之比，相对于常用的二阶五点FDFD方法和优化九点FDFD方法，通用格式自适应系数FDFD方法均可在相似的计算量和内存需求下，有效提高黏声波模拟的精度。

Description

一种基于自适应系数频域有限差分的黏声波模拟方法
技术领域
本发明涉及地震波场模拟领域，具体涉及一种基于自适应系数频域有限差分的黏声波模拟方法。
背景技术
由于地下介质的黏弹性，地震波在地下传播时存在着衰减效应。黏声波方程可以有效表征地震波传播的衰减特征，因此基于黏声波方程的波场数值模拟在地震资料处理、反演研究中具有重要应用，如基于黏声波方程的全波形反演和逆时偏移等。黏声波模拟的主要方法包括时域有限差分方法、频域有限差分(finite difference frequency domain,FDFD)方法、有限元法、伪谱法、傅里叶有限差分法等。在这些方法中，频域有限差分方法具有可灵活处理衰减效应、高效模拟多炮或窄带地震资料、避免时间频散、并行效率高等优点，因而在黏声波模拟中具有较多的应用案例。用FDFD方法模拟黏声波方程的主要挑战在于其需要求解大型的稀疏线性方程组，而大型稀疏线性方程组的求解往往需要花费较多的计算量和内存。克服这一挑战的一个常用方法是发展所需一个波长内的网格点数(numberof points per wavelength,NPPW)更少的FDFD方法，从而减少波场模拟所需的网格点数，进而减少相应大型稀疏线性方程组的阶数。
为求解一般的黏声波方程，传统的二阶五点FDFD方法所需NPPW较多。高阶FDFD方法可以有效减少所需的NPPW，但其对应的大型稀疏线性方程组含有较大的带宽和较多的非零元素，因此所需的计算量和内存较大。Arntsen等将Holberg采用十字形优化有限差分系数应用于黏声波方程的FDFD求解，并针对均匀模型应用解析解验证了其有效性，十字形优化FDFD相对高阶FDFD可以减少所需空间差分的长度，但相对紧凑格式的FDFD方法所需计算量和内存仍然较大，同时十字形优化FDFD在低波数区域的精度相对高阶FDFD会有所减少。通过将声波方程FDFD中的实波数替换为黏声波方程的复波数，Amini和Javaherian将Jo等关于声波方程的优化九点FDFD直接应用于黏声波方程，并给出了在一般黏声波模型中的测试，然而Jo等的FDFD只适用于正方形网格，同时Amini和Javaherian并没有通过黏声波方程的解析解或参考解验证相应FDFD方法在一般黏声波模型中的有效性。Chen在提出基于平均导数方法(average-derivative method,ADM)声波方程优化FDFD方法的同时，给出了在黏声波方程中的相应推广，但并未给出一般黏声波模型的数值算例。在将自适应系数FDFD应用到声波全波形反演的同时，Aghamiry等给出了自适应系数FDFD方法在黏声波方程的推广，但同样未给出一般黏声波模型的数值算例，同时Aghamiry等使用的基于旋转坐标系的自适应系数FDFD方法只适用于空间采样间隔相等的情形。
发明内容
目前FDFD黏声波模拟常用的二阶五点方法和优化九点方法在一个波长内的网格点数小于4时误差较大。为了进一步减少FDFD高精度黏声波模拟所需的一个波长内的网格点数从而减少所需的计算量和内存，本发明针对黏声波波场模拟发展了一种适用于不同空间采样间隔之比的通用格式自适应系数FDFD方法。同时，为了验证自适应系数FDFD方法对一般黏声波模型的有效性，本发明针对两个典型的黏声波模型，分别采用解析解和基于高阶FDFD的参考解验证了所提出方法的有效性。
为实现以上功能，本发明设计一种基于自适应系数频域有限差分的黏声波模拟方法，执行如下步骤S1-步骤S6，完成不同频率下的黏声波波场的模拟：
步骤S1：针对黏声波波场，基于黏声波方程，采用自适应系数频域有限差分方法，通过引入网格点内波长个数、空间采样间隔比，建立适用于不同空间采样间隔比的自适应系数频域有限差分通用格式方程；
步骤S2：确定黏声波方程在预设参考频率下的波传播速度模型和品质因子模型，以波传播速度模型和品质因子模型构建黏声波模型；
步骤S3：基于步骤S1所构建的自适应系数频域有限差分通用格式方程，以及步骤S2所确定的黏声波模型，根据黏声波模型的空间采样间隔比，建立根据网格点内波长个数查找频域有限差分自适应系数的查找表；
步骤S4：根据自适应系数频域有限差分通用格式方程的数值频散关系，验证所建立的查找表的有效性；
步骤S5：针对自适应系数频域有限差分通用格式方程，通过引入完美匹配层边界条件，消除人工边界反射，获得含有完美匹配层边界条件的自适应系数频域有限差分通用格式方程；
步骤S6：将含有完美匹配层边界条件的自适应系数频域有限差分通用格式方程应用于所有网格点，采用MPI多线程并行和PARDDISO直接求解器，求解自适应系数频域有限差分方法生成的大型稀疏线性方程组，获得不同频率下的黏声波波场，完成不同频率下的黏声波波场的模拟。
作为本发明的一种优选技术方案：步骤S1的具体方法如下：
针对黏声波波场，基于Kolsky-Futterman模型所表示的黏声波方程如下：
Figure BDA0004022712240000031
式中，u为频率域波场，ω为角频率，
Figure BDA0004022712240000032
vr表示在参考角频率ωr下的参考波传播速度，Q表示品质因子，sgn表示符号函数，i表示虚数单位；
网格点内波长个数、空间采样间隔的计算如下式：
Figure BDA0004022712240000033
r＝Δx/Δz
式中，
Figure BDA0004022712240000034
为沿x方向一个网格点内波长个数，r为空间采样间隔比，Δx和Δz分别表示沿x和z方向的空间采样间隔，λ表示vr对应的波长；
建立适用于不同空间采样间隔比的自适应系数频域有限差分通用格式方程如下式：
Figure BDA0004022712240000035
式中，um,n＝u(mΔx,nΔz)，(vc)m,n＝vc(mΔx,nΔz)，m和n分别表示波场模拟的网格点对应的下标，
Figure BDA0004022712240000036
表示自适应系数频域有限差分通用格式方程的自适应系数，其中ηi的下标i＝1,2,3。
作为本发明的一种优选技术方案：步骤S2中通过预设参考频率下的波传播速度模型获得品质因子模型Q方法如下：
Figure BDA0004022712240000037
式中，vp为预设参考频率下的波传播速度模型。
作为本发明的一种优选技术方案：步骤S3的具体方法如下：
假设黏声波方程频域有限差分的误差主要来源于数值频散，考虑品质因子模型Q的取值趋于无穷大时，黏声波方程的平面波解u(x,z)如下式：
Figure BDA0004022712240000041
式中，u0为常系数，i表示虚数单位，
Figure BDA0004022712240000042
表示实波数，θ表示波传播方向与z轴的夹角；
将平面波解u(x,z)代入自适应系数频域有限差分通用格式方程，并将vc取值vr，则有下式：
Figure BDA0004022712240000043
Figure BDA0004022712240000044
Figure BDA0004022712240000045
Figure BDA0004022712240000046
Figure BDA0004022712240000047
其中η＝[η1,η2,η3]，极小化平面波解的代入误差，所需的自适应系数
Figure BDA0004022712240000048
通过求解如下优化问题得到：/>
Figure BDA0004022712240000049
在黏声波波场模拟之前，对给定的空间采样间隔比率r，建立根据网格点内波长个数
Figure BDA00040227122400000410
查找自适应系数/>
Figure BDA00040227122400000411
的查找表，其中/>
Figure BDA00040227122400000412
以0.001为间隔在0.5到0.001之间进行离散采样，若自适应系数/>
Figure BDA00040227122400000413
不在查找表中，则通过线性插值方法获得。
作为本发明的一种优选技术方案：步骤S4的具体方法如下：
以vph表示自适应系数频域有限差分通用格式方程对应的数值相速度，则自适应系数频域有限差分方法的平面波解表示为下式：
Figure BDA0004022712240000051
为得到自适应系数频域有限差分通用格式方程对应的数值相速度vph和黏声波方程的参考波传播速度vr的偏差，将上述自适应系数频域有限差分方法的平面波解代入自适应系数频域有限差分通用格式方程，获得数值频散方程式如下式：
Figure BDA0004022712240000052
Figure BDA0004022712240000053
Figure BDA0004022712240000054
Figure BDA0004022712240000055
Figure BDA0004022712240000056
其中，
Figure BDA0004022712240000057
vph/vr表示自适应系数频域有限差分通用格式方程的数值频散关系，/>
Figure BDA0004022712240000058
作为本发明的一种优选技术方案：步骤S5的具体方法如下：
将带有完美匹配层吸收边界和震源项的二维黏声波方程表示为下式：
Figure BDA0004022712240000059
其中f(ω)表示频率域震源函数，(xs,zs)表示震源位置，δ(x)表示狄拉克函数，
Figure BDA00040227122400000510
τ∈{x,z}，R表示人工边界反射系数；vmax表示黏声波模型最大速度，Lτ表示沿τ方向完美匹配层的厚度，/>
Figure BDA00040227122400000511
表示网格点沿τ方向到内部区域的距离；
通过对上式中的
Figure BDA00040227122400000512
采用虚拟半网格点和二阶差分进行离散，结合自适应系数频域有限差分通用格式方程，可得含有完美匹配层边界条件的自适应系数频域有限差分通用格式方程如下式：
Figure BDA0004022712240000061
式中，(xm,zn)＝(mΔx,nΔz)，
Figure BDA0004022712240000062
和/>
Figure BDA0004022712240000063
表示克罗内克函数，(ms,ns)表示震源的网格点下标。
作为本发明的一种优选技术方案：步骤S6的具体方法如下：
在不同频率下，将含有完美匹配层边界条件的自适应系数频域有限差分通用格式方程应用于黏声波模型的所有网格点，则可得到关于频率域黏声波波场的大型稀疏线性方程组如下：
K(ω)u＝f(ω)
其中K(ω)和f(ω)分别表示由含有完美匹配层边界条件的自适应系数频域有限差分通用格式方程左端系数构造的随频率变化的大型稀疏系数矩阵、右端值构造的随频率变化的线性方程组右端项，u表示待求解的随频率变化的黏声波波场；采用MPI多进程并行将不同频率下的大型稀疏线性方程组分配给不同的计算机节点，然后每个计算机节点独立调用PARDDISO直接求解器求解所分配频率下的大型稀疏线性方程组。
有益效果：相对于现有技术，本发明的优点包括：
本发明设计了一种基于自适应系数频域有限差分的黏声波模拟方法，通过令FDFD系数随一个波长内的网格点数自适应从而提高FDFD方法的精度，本发明针对黏声波波场模拟发展了一种适用于不同空间采样间隔之比的通用格式自适应系数FDFD方法。相对目前黏声波模拟常用的二阶五点FDFD方法和优化九点FDFD方法，本发明的方法可以将一个波长内所需的网格点数从4减少到2.5，从而可以在相同的精度情况下，有效减少FDFD黏声波波场模拟所需的计算时间和计算机内存。
附图说明
图1是根据本发明实施例提供的一种基于自适应系数频域有限差分的黏声波模拟方法的流程示意图；
图2是空间采样间隔之比r＝1时不同FDFD方法的数值频散曲线对比图；
图3是空间采样间隔之比r＝3时不同FDFD方法的数值频散曲线对比图；
图4是通过MPI并行和PARDISO直接求解器求解由自适应系数FDFD生成的多个大型稀疏线性方程组得到黏声波模拟结果的流程示意图；
图5是均匀黏声波模型使用的r＝1时的自适应FDFD系数；
图6是均匀黏声波模型不同FDFD方法的接收地震记录对比图；
图7是均匀黏声波模型中二阶五点FDFD、ADM优化九点FDFD、通用格式自适应系数九点FDFD波场模拟结果与解析解单道对比图；
图8是Marmousi黏声波模型示意图；
图9是Marmousi黏声波模型使用的r＝3.125时的自适应FDFD系数；
图10是Marmousi黏声波模型不同FDFD方法的接收地震记录对比图；
图11是Marmousi黏声波模型中二阶五点FDFD、ADM优化九点FDFD、通用格式自适应系数九点FDFD模拟结果与参考解的单道对比图。
具体实施方式
下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。
参照图1，本发明实施例提供的一种基于自适应系数频域有限差分(finitedifference frequency domain,FDFD)的黏声波模拟方法，执行如下步骤S1-步骤S6，完成不同频率下的黏声波波场的模拟：
步骤S1：针对黏声波波场，基于黏声波方程，采用自适应系数频域有限差分方法，通过引入网格点内波长个数、空间采样间隔比，建立适用于不同空间采样间隔比的自适应系数频域有限差分通用格式方程。
步骤S1的具体方法如下：
针对黏声波波场，频率域二维无源黏声波方程如下：
Figure BDA0004022712240000071
式中，u为频率域波场，ω为角频率，vc为复速度；
复速度vc采用如下Kolsky-Futterman模型进行描述：
Figure BDA0004022712240000081
其中vr表示在参考角频率ωr下的参考波传播速度，Q表示品质因子，sgn表示符号函数，i表示虚数单位。本发明将参考角频率对应的频率值设为50Hz。需指出的是，上式对应的一维平面波解的形式是
Figure BDA0004022712240000082
如果对应的一维平面波解的形式是/>
Figure BDA0004022712240000083
则上式虚部的符号将改变。
网格点内波长个数、空间采样间隔的计算如下式：
Figure BDA0004022712240000084
r＝Δx/Δz
式中，
Figure BDA0004022712240000085
为沿x方向一个网格点内波长个数，r为空间采样间隔比，Δx和Δz分别表示沿x和z方向的空间采样间隔，λ表示vr对应的波长；
自适应系数频域有限差分通用格式方程可以在传统二阶频域有限差分格式的基础上加入相关的系数随
Figure BDA0004022712240000086
r变化的校正项得到，其中校正项可以按网格点与中心点的距离进行分类选取。
适用于不同空间采样间隔比的自适应系数频域有限差分通用格式方程，即自适应系数九点FDFD通用格式如下式：
Figure BDA0004022712240000087
式中，um,n＝u(mΔx,nΔz)，(vc)m,n＝vc(mΔx,nΔz)，m和n分别表示波场模拟的网格点对应的下标，
Figure BDA0004022712240000088
表示自适应系数频域有限差分通用格式方程的自适应系数，其中ηi的下标i＝1,2,3。
步骤S2：确定黏声波方程在预设参考频率下的波传播速度模型和品质因子模型，以波传播速度模型和品质因子模型构建黏声波模型。
本发明将参考角频率对应的频率值设为50Hz。假设在进行黏声波波场模拟时，参考频率下的波传播速度模型已被提供(如从地震数据中通过全波形反演等手段反演得到速度模型)。关于品质因子模型Q，若在给出波传播速度模型时未同时提供，本发明通过速度模型换算得到品质因子模型Q，公式如下：
Figure BDA0004022712240000091
式中，vp为预设参考频率下的波传播速度模型。
步骤S3：基于步骤S1所构建的自适应系数频域有限差分通用格式方程，以及步骤S2所确定的黏声波模型，根据黏声波模型的空间采样间隔比，建立根据网格点内波长个数查找频域有限差分自适应系数的查找表。
步骤S3的具体方法如下：
本发明假设黏声波频域有限差分的误差主要来源于数值频散，因此在确定自适应系数的过程中，只考虑品质因子趋于无穷大的情况(即Q→∞)。此时黏声波方程的平面波解与声波方程的平面波解具有相同形式，黏声波方程的平面波解u(x,z)如下式：
Figure BDA0004022712240000092
式中，u0为常系数，i表示虚数单位，
Figure BDA0004022712240000093
表示实波数，θ表示波传播方向与z轴的夹角；
将平面波解u(x,z)代入自适应系数频域有限差分通用格式方程，并将vc取值vr，则有下式：
Figure BDA0004022712240000094
Figure BDA0004022712240000095
Figure BDA0004022712240000096
Figure BDA0004022712240000097
Figure BDA0004022712240000098
其中η＝[η1,η2,η3]，极小化平面波解的代入误差，所需的自适应系数
Figure BDA0004022712240000099
通过求解如下优化问题得到：
Figure BDA00040227122400000910
Figure BDA0004022712240000101
上式中的优化问题可通过Polak-Ribière-Polyak共轭梯度法求解：
Figure BDA0004022712240000102
Figure BDA0004022712240000103
ηk-1＝ηk+αkdk
其中k表示迭代的次数，
Figure BDA0004022712240000104
表示梯度向量，dk表示下降方向，αk表示dk对应的迭代步长，/>
Figure BDA0004022712240000105
ηk是自适应系数频域有限差分通用格式方程的自适应系数对应的迭代变量。
在黏声波波场模拟之前，对给定的空间采样间隔比率r，建立根据网格点内波长个数
Figure BDA0004022712240000106
查找自适应系数/>
Figure BDA0004022712240000107
的查找表，其中/>
Figure BDA0004022712240000108
以0.001为间隔在0.5到0.001之间进行离散采样，若自适应系数/>
Figure BDA0004022712240000109
不在查找表中，则通过线性插值方法获得。
步骤S4：根据自适应系数频域有限差分通用格式方程的数值频散关系，验证所建立的查找表的有效性。
步骤S4的具体方法如下：
使用频域有限差分方法数值求解波动方程时容易导致数值频散，即频域有限差分方法对应的相速度与模型真实速度有偏差。以vph表示自适应系数频域有限差分通用格式方程对应的数值相速度，则自适应系数频域有限差分方法的平面波解(和黏声波方程的平面波解不同)表示为下式：
Figure BDA00040227122400001010
为得到自适应系数频域有限差分通用格式方程对应的数值相速度vph和黏声波方程的参考波传播速度vr的偏差，将上述自适应系数频域有限差分方法的平面波解代入自适应系数频域有限差分通用格式方程，获得数值频散方程式如下式：
Figure BDA00040227122400001011
Figure BDA00040227122400001012
Figure BDA0004022712240000111
Figure BDA0004022712240000112
Figure BDA0004022712240000113
Figure BDA0004022712240000114
其中，
Figure BDA0004022712240000115
vph/vr表示自适应系数频域有限差分通用格式方程的数值频散关系，/>
Figure BDA0004022712240000116
与常系数频域有限差分不同的是，自适应系数频域有限差分的自适应系数是随
Figure BDA0004022712240000117
变化的，所以自适应系数频域有限差分并没有解析的数值频散表达式，即vph/vr无法表示为与vr和ω无关的形式。为了给出自适应系数频域有限差分的数值频散关系，本发明假设对给定的/>
Figure BDA0004022712240000118
和r，在/>
Figure BDA0004022712240000119
的预设邻域范围内存在且只存在一个满足数值频散方程式的/>
Figure BDA00040227122400001110
在这个假设下，相应的/>
Figure BDA00040227122400001111
可通过在/>
Figure BDA00040227122400001112
的邻域内极小化数值频散方程式两端差的平方得到。且自适应系数/>
Figure BDA00040227122400001113
通过查找表得到。通过这种方式，可得自适应系数频域有限差分通用格式方程的自适应系数的数值频散关系
Figure BDA00040227122400001114
具体地，为得到数值频散关系vph/vr，求解等价的/>
Figure BDA00040227122400001115
按照通常的频域有限差分数值频散分析，/>
Figure BDA00040227122400001116
通常已知且为[0,0.5]区间的等间隔采样离散向量。通过已知的黏声波模型得到空间采样间隔之比r后，通过在/>
Figure BDA00040227122400001117
的邻域内极小化数值频散方程式两端差的平方得到相应的/>
Figure BDA00040227122400001118
得到所需的/>
Figure BDA00040227122400001119
图2和图3分别给出空间采样间隔比率r＝1和r＝3时自适应系数九点FDFD和二阶五点FDFD、优化ADM九点FDFD的数值频散曲线对比。图2和图3中(a)图为二阶五点FDFD；(b)图为ADM优化九点FDFD；(c)图为通用格式自适应系数九点FDFD；图2和图3结果表明，相对于二阶五点FDFD和ADM优化九点FDFD，通用格式自适应系数九点FDFD可以更为有效的压制正方形网格和非正方形网格下的数值频散。以数值频散误差不超过1％为标准，相对于二阶五点FDFD和ADM优化九点FDFD，新方法可将黏声波方程一个波长内所需的网格点数从4减少到2.5。
步骤S5：针对自适应系数频域有限差分通用格式方程，通过引入完美匹配层(perfectly matched layer,PML)边界条件，消除人工边界反射，获得含有完美匹配层边界条件的自适应系数频域有限差分通用格式方程。
步骤S5的具体方法如下：
本发明采用完美匹配层(perfectly matched layer，PML)吸收边界条件压制人工边界反射。将带有完美匹配层吸收边界和震源项的二维黏声波方程表示为下式：
Figure BDA0004022712240000121
其中f(ω)表示频率域震源函数，(xs,zs)表示震源位置，δ(x)表示狄拉克函数，
Figure BDA0004022712240000122
τ∈{x,z}，R表示人工边界反射系数(本发明取为0.001)；vmax表示黏声波模型最大速度，Lτ表示沿τ方向完美匹配层的厚度，/>
Figure BDA0004022712240000123
表示网格点沿τ方向到内部区域的距离；
通过对上式中的
Figure BDA0004022712240000124
采用虚拟半网格点和二阶差分进行离散，结合自适应系数频域有限差分通用格式方程，可得含有完美匹配层边界条件的自适应系数频域有限差分通用格式方程如下式：
Figure BDA0004022712240000125
式中，(xm,zn)＝(mΔx,nΔz)，
Figure BDA0004022712240000126
和/>
Figure BDA0004022712240000127
表示克罗内克函数，(ms,ns)表示震源的网格点下标。
步骤S6：将含有完美匹配层边界条件的自适应系数频域有限差分通用格式方程应用于所有网格点，采用MPI多线程并行和PARDDISO直接求解器，求解自适应系数频域有限差分方法生成的大型稀疏线性方程组，获得不同频率下的黏声波波场，完成不同频率下的黏声波波场的模拟。
步骤S6的具体方法如下：
在不同频率下，将含有完美匹配层边界条件的自适应系数频域有限差分通用格式方程应用于黏声波模型的所有网格点，则可得到关于频率域黏声波波场的大型稀疏线性方程组如下：
K(ω)u＝f(ω)
其中K(ω)和f(ω)分别表示由含有完美匹配层边界条件的自适应系数频域有限差分通用格式方程左端系数构造的随频率变化的大型稀疏系数矩阵、右端值构造的随频率变化的线性方程组右端项，u表示待求解的随频率变化的黏声波波场；采用MPI多进程并行将不同频率下的大型稀疏线性方程组分配给不同的计算机节点，然后每个计算机节点独立调用PARDDISO直接求解器求解所分配频率下的大型稀疏线性方程组。图4给出了通过MPI并行和PARDISO直接求解器求解由自适应系数频域有限差分方法生成的多个大型稀疏线性方程组得到黏声波模拟结果的流程示意图。
以下为本发明的两个实施例，说明基于一种基于自适应系数频域有限差分的黏声波模拟方法的实现过程。对于所有实施例，FDFD生成的大型稀疏线性方程组均通过PARDISO直接求解器求解，计算平台均为超算平台的30个节点(AMD EPYC 7452@2.35GHz)，每个节点有64核和256GB的内存，不同频率的波场通过MPI多节点并行计算。
首先考虑网格点数为201×201的均匀模型。模型在50Hz参考频率下的速度和品质因子分别为2000m/s和50。模型的采样间隔为Δx＝Δz＝16.7m，各侧PML吸收边界层数均为20，震源选为主频20Hz的雷克子波，震源位置在模型中心。检波器布置在地表所有网格点处。波场模拟的频率范围是0.5Hz-60 Hz且频率采样间隔为0.5Hz。从频域地震记录转换为时域地震记录的时间范围是0到2s且时间采样间隔为1ms。图5给出本实施例使用的部分自适应FDFD系数(r＝1)，表明所使用的自适应FDFD系数是较光滑的。图6给出了均匀黏声波模型不同FDFD方法的接收地震记录对比，其中图(a)为二阶五点FDFD、图(b)为ADM优化九点FDFD、图(c)为自适应系数九点FDFD。图7给出了在(x,z)＝(0m,0m)处上述三种FDFD接收的地震记录与解析解的对比，其中图(a)为绝对振幅比较，图(b)为误差比较。图6和图7表明对这个均匀黏声波模型，相对于二阶五点FDFD和ADM优化九点FDFD，本发明所设计的自适应系数九点FDFD可以更有效地减少黏声波模拟误差。
第二个实施例考虑一个网格点数为737×750的Marmousi黏声波模型。在50Hz参考频率下，模型速度如图8(a)所示，模型的品质因子如图8(b)所示，是通过速度模型换算得到。模型的采样间隔为Δx＝12.5m、Δz＝4m，沿x和z方向的单侧PML吸收边界层数分别为20和40，震源选为主频20Hz的雷克子波，震源位置为(xs,zs)＝(4600m,40m)。检波器布置在地表所有网格点处。波场模拟频率范围是0.125Hz-60Hz且频率采样间隔为0.125Hz。从频域地震记录转换为时域地震记录的时间范围是0-8s且时间采样间隔为1ms。本实施例采用沿x方向72阶、沿z方向8阶的高阶FDFD作为参考解。图9给出了本实施例使用的部分自适应FDFD系数(r＝3.125)。图10给出了各检波器接收的地震记录的对比，其中图(a)为二阶五点FDFD；图(b)为ADM优化九点FDFD；图(c)为通用格式自适应系数九点FDFD；图(d)为参考解(沿x方向72阶、沿z方向8阶的高阶FDFD)。黑色箭头指向差异较大的区域。图11给出了在(x,z)＝(0m,0m)处上述三种FDFD接受的地震记录与参考解(沿x方向72阶、沿z方向8阶的高阶FDFD)的对比，其中图(a)为绝对振幅；图(b)为误差。图10和图11的结果表明，对该Marmousi黏声波模型，相对于二阶五点FDFD和ADM优化九点FDFD，自适应系数九点FDFD可以更为有效的减少模拟误差。同时上述四种FDFD在超算30个节点的计算时间分别为23s、27s、28s、947s，所需内存分别为2.4GB、3.0GB、3.0GB、47.5GB，通用格式自适应系数九点FDFD的计算时间和所需内存与另两种方法相当，远小于参考解所需时间和内存。
上面结合附图对本发明的实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。

Claims (7)

1.一种基于自适应系数频域有限差分的黏声波模拟方法，其特征在于，执行如下步骤S1-步骤S6，完成不同频率下的黏声波波场的模拟：
步骤S1：针对黏声波波场，基于黏声波方程，采用自适应系数频域有限差分方法，通过引入网格点内波长个数、空间采样间隔比，建立适用于不同空间采样间隔比的自适应系数频域有限差分通用格式方程；
步骤S2：确定黏声波方程在预设参考频率下的波传播速度模型和品质因子模型，以波传播速度模型和品质因子模型构建黏声波模型；
步骤S3：基于步骤S1所构建的自适应系数频域有限差分通用格式方程，以及步骤S2所确定的黏声波模型，根据黏声波模型的空间采样间隔比，建立根据网格点内波长个数查找频域有限差分自适应系数的查找表；
步骤S4：根据自适应系数频域有限差分通用格式方程的数值频散关系，验证所建立的查找表的有效性；
步骤S5：针对自适应系数频域有限差分通用格式方程，通过引入完美匹配层边界条件，消除人工边界反射，获得含有完美匹配层边界条件的自适应系数频域有限差分通用格式方程；
步骤S6：将含有完美匹配层边界条件的自适应系数频域有限差分通用格式方程应用于所有网格点，采用MPI多线程并行和PARDDISO直接求解器，求解自适应系数频域有限差分方法生成的大型稀疏线性方程组，获得不同频率下的黏声波波场，完成不同频率下的黏声波波场的模拟。
2.根据权利要求1所述的一种基于自适应系数频域有限差分的黏声波模拟方法，其特征在于，步骤S1的具体方法如下：
针对黏声波波场，基于Kolsky-Futterman模型所表示的黏声波方程如下：
Figure FDA0004022712230000011
式中，u为频率域波场，ω为角频率，
Figure FDA0004022712230000012
vr表示在参考角频率ωr下的参考波传播速度，Q表示品质因子，sgn表示符号函数，i表示虚数单位；
网格点内波长个数、空间采样间隔的计算如下式：
Figure FDA0004022712230000013
r＝Δx/Δz
式中，
Figure FDA0004022712230000021
为沿x方向一个网格点内波长个数，r为空间采样间隔比，Δx和Δz分别表示沿x和z方向的空间采样间隔，λ表示vr对应的波长；
建立适用于不同空间采样间隔比的自适应系数频域有限差分通用格式方程如下式：
Figure FDA0004022712230000022
式中，um,n＝u(mΔx,nΔz)，(vc)m,n＝vc(mΔx,nΔz)，m和n分别表示波场模拟的网格点对应的下标，
Figure FDA0004022712230000023
表示自适应系数频域有限差分通用格式方程的自适应系数，其中ηi的下标i＝1,2,3。
3.根据权利要求2所述的一种基于自适应系数频域有限差分的黏声波模拟方法，其特征在于，步骤S2中通过预设参考频率下的波传播速度模型获得品质因子模型Q方法如下：
Figure FDA0004022712230000024
式中，vp为预设参考频率下的波传播速度模型。
4.根据权利要求3所述的一种基于自适应系数频域有限差分的黏声波模拟方法，其特征在于，步骤S3的具体方法如下：
假设黏声波方程频域有限差分的误差主要来源于数值频散，考虑品质因子模型Q的取值趋于无穷大时，黏声波方程的平面波解u(x,z)如下式：
Figure FDA0004022712230000025
式中，u0为常系数，i表示虚数单位，
Figure FDA0004022712230000026
表示实波数，θ表示波传播方向与z轴的夹角；
将平面波解u(x,z)代入自适应系数频域有限差分通用格式方程，并将vc取值vr，则有下式：
Figure FDA0004022712230000031
Figure FDA0004022712230000032
Figure FDA0004022712230000033
Figure FDA0004022712230000034
Figure FDA0004022712230000035
其中η＝[η1,η2,η3]，极小化平面波解的代入误差，所需的自适应系数
Figure FDA0004022712230000036
通过求解如下优化问题得到：/>
Figure FDA0004022712230000037
在黏声波波场模拟之前，对给定的空间采样间隔比率r，建立根据网格点内波长个数
Figure FDA0004022712230000038
查找自适应系数/>
Figure FDA0004022712230000039
的查找表，其中/>
Figure FDA00040227122300000310
以0.001为间隔在0.5到0.001之间进行离散采样，若自适应系数/>
Figure FDA00040227122300000311
不在查找表中，则通过线性插值方法获得。
5.根据权利要求4所述的一种基于自适应系数频域有限差分的黏声波模拟方法，其特征在于，步骤S4的具体方法如下：
以vph表示自适应系数频域有限差分通用格式方程对应的数值相速度，则自适应系数频域有限差分方法的平面波解表示为下式：
Figure FDA00040227122300000312
为得到自适应系数频域有限差分通用格式方程对应的数值相速度vph和黏声波方程的参考波传播速度vr的偏差，将上述自适应系数频域有限差分方法的平面波解代入自适应系数频域有限差分通用格式方程，获得数值频散方程式如下式：
Figure FDA00040227122300000313
Figure FDA00040227122300000314
Figure FDA0004022712230000041
Figure FDA0004022712230000042
Figure FDA0004022712230000043
其中，
Figure FDA0004022712230000044
vph/vr表示自适应系数频域有限差分通用格式方程的数值频散关系，/>
Figure FDA0004022712230000045
6.根据权利要求5所述的一种基于自适应系数频域有限差分的黏声波模拟方法，其特征在于，步骤S5的具体方法如下：
将带有完美匹配层吸收边界和震源项的二维黏声波方程表示为下式：
Figure FDA0004022712230000046
其中f(ω)表示频率域震源函数，(xs,zs)表示震源位置，δ(x)表示狄拉克函数，
Figure FDA0004022712230000047
R表示人工边界反射系数；vmax表示黏声波模型最大速度，Lτ表示沿τ方向完美匹配层的厚度，/>
Figure FDA00040227122300000412
表示网格点沿τ方向到内部区域的距离；/>
通过对上式中的
Figure FDA0004022712230000048
采用虚拟半网格点和二阶差分进行离散，结合自适应系数频域有限差分通用格式方程，可得含有完美匹配层边界条件的自适应系数频域有限差分通用格式方程如下式：
Figure FDA0004022712230000049
式中，(xm,zn)＝(mΔx,nΔz)，
Figure FDA00040227122300000410
和/>
Figure FDA00040227122300000411
表示克罗内克函数，(ms,ns)表示震源的网格点下标。
7.根据权利要求6所述的一种基于自适应系数频域有限差分的黏声波模拟方法，其特征在于，步骤S6的具体方法如下：
在不同频率下，将含有完美匹配层边界条件的自适应系数频域有限差分通用格式方程应用于黏声波模型的所有网格点，则可得到关于频率域黏声波波场的大型稀疏线性方程组如下：
K(ω)u＝f(ω)
其中K(ω)和f(ω)分别表示由含有完美匹配层边界条件的自适应系数频域有限差分通用格式方程左端系数构造的随频率变化的大型稀疏系数矩阵、右端值构造的随频率变化的线性方程组右端项，u表示待求解的随频率变化的黏声波波场；采用MPI多进程并行将不同频率下的大型稀疏线性方程组分配给不同的计算机节点，然后每个计算机节点独立调用PARDDISO直接求解器求解所分配频率下的大型稀疏线性方程组。
CN202211697099.6A 2022-12-28 2022-12-28 一种基于自适应系数频域有限差分的黏声波模拟方法 Active CN115932949B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211697099.6A CN115932949B (zh) 2022-12-28 2022-12-28 一种基于自适应系数频域有限差分的黏声波模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211697099.6A CN115932949B (zh) 2022-12-28 2022-12-28 一种基于自适应系数频域有限差分的黏声波模拟方法

Publications (2)

Publication Number Publication Date
CN115932949A CN115932949A (zh) 2023-04-07
CN115932949B true CN115932949B (zh) 2023-06-06

Family

ID=85838374

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211697099.6A Active CN115932949B (zh) 2022-12-28 2022-12-28 一种基于自适应系数频域有限差分的黏声波模拟方法

Country Status (1)

Country Link
CN (1) CN115932949B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117494533B (zh) * 2024-01-02 2024-03-12 电子科技大学 一种涉及磁边界的四分量波导端口特性求解方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107479092A (zh) * 2017-08-17 2017-12-15 电子科技大学 一种基于方向导数的频率域高阶声波方程正演模拟方法
CN112817044A (zh) * 2020-10-16 2021-05-18 中国石油大学(华东) 频率域声波方程的差分系数优化方法、系统

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2548729B (en) * 2014-12-31 2021-02-17 Landmark Graphics Corp Seismic elastic wave simulation for tilted tranversely isotropic media using adaptive Lebedev staggered grid

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107479092A (zh) * 2017-08-17 2017-12-15 电子科技大学 一种基于方向导数的频率域高阶声波方程正演模拟方法
CN112817044A (zh) * 2020-10-16 2021-05-18 中国石油大学(华东) 频率域声波方程的差分系数优化方法、系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
二维黏滞声波方程的优化组合型紧致有限差分数值模拟;汪勇;徐佑德;高刚;桂志先;陈英;王亚楠;;石油地球物理勘探(第06期);全文 *

Also Published As

Publication number Publication date
CN115932949A (zh) 2023-04-07

Similar Documents

Publication Publication Date Title
Manolis A comparative study on three boundary element method approaches to problems in elastodynamics
Min et al. Improved frequency-domain elastic wave modeling using weighted-averaging difference operators
CN107783190A (zh) 一种最小二乘逆时偏移梯度更新方法
CN115932949B (zh) 一种基于自适应系数频域有限差分的黏声波模拟方法
CN110109177B (zh) 基于旋转时空双变网格有限差分法的地震波正演模拟方法
CN110857998B (zh) 一种基于lowrank有限差分的弹性逆时偏移方法及系统
Le Bouteiller et al. A discontinuous Galerkin fast-sweeping eikonal solver for fast and accurate traveltime computation in 3D tilted anisotropic media
Wang et al. Application of a Chebyshev collocation method to solve a parabolic equation model of underwater acoustic propagation
CN110488354B (zh) 一种q补偿的起伏地表棱柱波与一次波联合最小二乘逆时偏移成像方法
CN110596754A (zh) 一种三维TTI介质qP波与qSV波波场模拟方法
CN111665556A (zh) 地层声波传播速度模型构建方法
CN111257930B (zh) 一种黏弹各向异性双相介质区域变网格求解算子
CN108828659A (zh) 地震波场延拓方法及装置
CN105403919B (zh) 一种逆时偏移成像方法及装置
Ren et al. Time-dispersion correction for arbitrary even-order Lax-Wendroff methods and the application on full-waveform inversion
CN110658558A (zh) 吸收衰减介质叠前深度逆时偏移成像方法及系统
Kucherov et al. High‐order absorbing boundary conditions incorporated in a spectral element formulation
AlSalem et al. Embedded boundary methods for modeling 3D finite-difference Laplace-Fourier domain acoustic-wave equation with free-surface topography
Gonza´ lez et al. Weak coupling of multibody dynamics and block diagram simulation tools
CN113311484B (zh) 利用全波形反演获取粘弹性介质弹性参数的方法及装置
Yu et al. Simulations of ground motions under plane wave incidence in 2D complex site based on the spectral element method (SEM) and multi-transmitting formula (MTF): SH problem
CN113221392B (zh) 一种非均匀粘声波在无限域内传播模型的构建方法
Dan et al. Double fictitious background media formulation for the Helmholtz equation in inhomogeneous media
CN111665550A (zh) 地下介质密度信息反演方法
CN113139266B (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