CN112906267A - 频率域粘弹性波正演模拟方法、装置、设备及存储介质 - Google Patents

频率域粘弹性波正演模拟方法、装置、设备及存储介质 Download PDF

Info

Publication number
CN112906267A
CN112906267A CN202110177528.6A CN202110177528A CN112906267A CN 112906267 A CN112906267 A CN 112906267A CN 202110177528 A CN202110177528 A CN 202110177528A CN 112906267 A CN112906267 A CN 112906267A
Authority
CN
China
Prior art keywords
wave
frequency
dimensional
wave number
domain
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.)
Pending
Application number
CN202110177528.6A
Other languages
English (en)
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 University of Petroleum Beijing
Original Assignee
China University of Petroleum Beijing
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 University of Petroleum Beijing filed Critical China University of Petroleum Beijing
Priority to CN202110177528.6A priority Critical patent/CN112906267A/zh
Publication of CN112906267A publication Critical patent/CN112906267A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本说明书实施例提供了一种频率域粘弹性波正演模拟方法、装置、设备及存储介质，该方法包括：获取频率波数域一维粘弹性波动方程，在将其转换为等效积分弱形式后，离散化为有限元控制方程；加载目标工区的三维点力源和吸收边界条件至有限元控制方程；对加载有三维点力源和吸收边界条件的有限元控制方程进行求解，获得目标工区的频率波数域波场；根据自适应波数采样策略确定频率波数域波场对应的波数集合；对波数集合中每个波数作为输入时对应的频率波数域波场，进行非等间隔二维傅里叶反变换，获得对应的频率域三维弹性波波场。本说明书实施例可以提高频率域粘弹性波正演模拟的精度。

Description

频率域粘弹性波正演模拟方法、装置、设备及存储介质
技术领域
本说明书涉及粘弹性波正演模拟技术领域，尤其是涉及一种频率域粘弹性波正演模拟方法、装置、设备及存储介质。
背景技术
针对垂向变化粘弹性介质的波传播模拟在勘探地震领域有着十分广泛的应用。比如，研究波在层状介质中的反射和透射规律、在积分方程法中获取背景场等均涉及到一维模型的全波场计算。在早期，半解析法被广泛应用于研究地震波在一维层状介质中的传播特征。根据Muller(1985)，用于研究水平层状介质波传播模拟和合成人工地震记录的半解析法主要可以分为射线理论、全波形理论、反射率法以及波数累加法。在这些方法中，反射率法应用最为广泛。该方法由Thomseon(1950)最早提出，之后诸多学者对该方法进行了改良和完善(Kennett,1980；Schimide和Tango,1986；Mallick和Frazer)。反射率法的最大优势在于它能给出三维点源在均匀层状模型中激发的全波场解析表达式。然而，尽管反射律法发展十分完善，但由于该方法涉及复杂的汉克尔变换，因此导致不同算法的计算精度和计算效率差距显著。此外，半解析法仅适用于层状均匀介质，对于具有垂向变化特征的粘弹性各向异性介质，通常难以推导其解析表达式。在这种情况下，数值方法显得尤为重要。
目前，大多数的数值模拟方法均针对二维和三维情况。在处理具有复杂分布的实际介质或地形条件时，这些二维或三维算法具有显著优势。然而，采用二维和三维正演模拟方法对垂向变化粘弹性介质进行波传播模拟，则计算代价过高。为此，JafarGandomi和Takenaka(2007)提出了一种模拟一维衰减介质平面波场的有限差分方法，该方法利用广义Zener模型描述介质衰减，并通过在τ-p域求解波动方程获得各个波场分量。之后，JafarGandomi和Takenaka(2013)将该方法进行了拓展，使其能够处理速度变化剧烈的一维模型。然而，这种方法仅能获得平面波解，而无法得到完整的波场解；此外，该方法只局限于各向同性介质。针对垂向变化各向异性粘弹性介质的数值方法对研究地震波在复杂一维介质中的传播规律具有重要价值，而目前仍被普遍忽略。因此，如何提高垂向变化粘弹性介质的波传播模拟的精度，已成为目前亟待解决的技术问题。
发明内容
本说明书实施例的目的在于提供一种频率域粘弹性波正演模拟方法、装置及设备，以提高频率域粘弹性波正演模拟的精度。
为达到上述目的，一方面，本说明书实施例提供了一种频率域粘弹性波正演模拟方法，包括：
获取频率波数域一维粘弹性波动方程；
将所述频率波数域一维粘弹性波动方程转换为等效积分弱形式；
将频率波数域一维粘弹性波动方程的等效积分弱形式离散化为有限元控制方程；
加载目标工区的三维点力源和吸收边界条件至所述有限元控制方程；
对加载有所述三维点力源和所述吸收边界条件的有限元控制方程进行求解，获得所述目标工区的频率波数域波场；
根据自适应波数采样策略确定所述频率波数域波场对应的波数集合；
对所述波数集合中每个波数作为输入时对应的频率波数域波场，进行非等间隔二维傅里叶反变换，获得对应的频率域三维弹性波波场。
本说明书一实施例中，所述频率波数域一维粘弹性波动方程，包括：
Figure BDA0002941127010000021
其中，
Figure BDA0002941127010000022
分别为σ,e,f,L,u关于x和y方向的二维傅里叶变换，L为偏微分算子矩阵，且
Figure BDA0002941127010000023
σ为应力向量；x,y,z分别为空间三维方向；f为震源矢量；ρ为介质密度；ω为角频率，u为频率域位移矢量；e为应变向量；C为粘弹性介质的复刚度矩阵；T为矩阵转置；
Figure BDA0002941127010000024
分别为关于x,y,z 方向的偏导。
本说明书一实施例中，所述频率波数域一维粘弹性波动方程的等效积分弱形式为：
Figure BDA0002941127010000031
其中，Γ为计算区域；
Figure BDA0002941127010000032
为
Figure BDA0002941127010000033
的变分；i为虚数单位；kx为x方向的波数；ky为y方向的波数。
本说明书一实施例中，所述有限元控制方程包括：
Figure BDA0002941127010000034
其中，K为全局刚度矩阵；M为全局质量矩阵；F为全局震源矢量；K,M,F的单元矩阵分别为
Figure BDA0002941127010000035
Figure BDA0002941127010000036
N为形函数矩阵；Γe为与第e个离散单元对应的计算区域。
本说明书一实施例中，所述加载有所述三维点力源和所述吸收边界条件的有限元控制方程，包括：
Figure BDA0002941127010000037
其中，F的单元矩阵为
Figure BDA0002941127010000038
Figure BDA0002941127010000039
为f的方向余弦，
Figure BDA00029411270100000310
为f与z轴夹角，θ为f在xy平面的投影与x轴的夹角；δ为狄拉克函数；A为震源幅值；R(f)为频率域雷克子波函数；f为雷克子波的频率；zs为震源的z轴坐标；D为阻尼矩阵，且D＝DMM；DM为质量阻尼系数，且DM＝DMmaxZ(z′)3；DMmax为质量阻尼系数极大值，且等于入射波的角频率；z′为刚性弱化SRM区域的计算点到内边界的距离；
Figure BDA00029411270100000311
C满足
Figure BDA00029411270100000312
C(z′)为刚性弱化区域的计算点到内边界的距离为z′时对应的粘弹性介质的复刚度矩阵；α(z′)＝αmaxZ(z′)p；
Figure BDA00029411270100000313
p为常数；kinc为入射波的波数；LSRM为SRM吸收层的厚度；rC为SRM吸收层内边界与外边界弹性刚度的比值；C0为SRM区域内边界的刚度矩阵；Re[C0]为C0的实部；e为自然常数。
本说明书一实施例中，所述自适应波数采样策略包括：
将波数区间
Figure BDA00029411270100000314
内的三角单元的平均尺寸设置为：
Figure BDA00029411270100000315
最大波数设置为：kmax＝10kc；
在波数区间
Figure BDA0002941127010000041
内采用非规则网格，网格密度由
Figure BDA0002941127010000042
至
Figure BDA0002941127010000043
逐渐降低；
三角单元的最大尺寸设置为：(Δk)max＝0.9kc；
其中，Δk为三角单元的平均尺寸；Nk为波数k的采样点数量；max为最大值函数；int为取整函数；rmax为最大炮检距；λmax为最大波长；
Figure BDA0002941127010000044
Vmin为最小速度。
本说明书一实施例中，所述频率域三维弹性波波场包括：
Figure BDA0002941127010000045
其中，u(ξ,ζ)为(ξ,ζ)下的频率域三维弹性波波场；n为二维波数域中第n个离散节点的编号，n＝1,2,…,Nd；Nd为离散节点总数；Wn为第n个离散节点处的积分权系数；Ne为二维波数域中划分出的离散单元数量；C＝1/(2π)；m＝1,2,…,6分别为第e个离散单元的6个节点编号；
Figure BDA0002941127010000046
为第e个三角单元中第m个离散节点的权系数；
Figure BDA0002941127010000047
为第e个三角单元中第m个离散节点的位移分量，
Figure BDA0002941127010000048
为第n个离散节点的位移分量。
另一方面，本说明书实施例还提供了一种频率域粘弹性波正演模拟装置，包括：
获取模块，用于获取频率波数域一维粘弹性波动方程；
转换模块，用于将所述频率波数域一维粘弹性波动方程转换为等效积分弱形式；
离散模块，用于将频率波数域一维粘弹性波动方程的等效积分弱形式离散化为有限元控制方程；
加载模块，用于加载目标工区的三维点力源和吸收边界条件至所述有限元控制方程；
求解模块，用于对加载有所述三维点力源和所述吸收边界条件的有限元控制方程进行求解，获得所述目标工区的频率波数域波场；
确定模块，用于根据自适应波数采样策略确定所述频率波数域波场对应的波数集合；
变换模块，用于对所述波数集合中每个波数作为输入时对应的频率波数域波场，进行非等间隔二维傅里叶反变换，获得对应的频率域三维弹性波波场。
另一方面，本说明书实施例还提供了一种计算机设备，包括存储器、处理器、以及存储在所述存储器上的计算机程序，所述计算机程序被所述处理器运行时，执行上述方法的指令。
另一方面，本说明书实施例还提供了一种计算机存储介质，其上存储有计算机程序，所述计算机程序被计算机设备的处理器运行时，执行上述方法的指令。
由以上本说明书实施例提供的技术方案可见，在本说明书实施例中，可以根据自适应波数采样策略确定频率波数域波场对应的波数集合，并对波数集合中每个波数作为输入时对应的频率波数域波场，进行非等间隔二维傅里叶反变换，获得对应的频率域三维弹性波波场；如此，不仅考虑了频率波数域粘弹性波场的振荡效应，并且辅之以自适应波数采样策略，从而实现三维点力源在任意一维各向异性粘弹性介质中的波场模拟；因此，与现有技术相比，本说明书实施例能够大大提升震源附近粘弹性波场的模拟精度。
附图说明
为了更清楚地说明本说明书实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本说明书中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。在附图中：
图1示出了本说明书一些实施例中频率域粘弹性波正演模拟方法的流程图；
图2示出了本说明书一实施例中刚性弱化(Strength Reduction Method，SRM) 吸收边界示意图；
图3a～图3d示出了本说明书一实施例中二维波数域波场的变化特征图；
图4示出了说明书一实施例中二维波数空间线单元剖分示意图；
图5示出了本说明书一实施例中各向异性全空间模型示意图；
图6示出了本说明书一实施例中各向异性全空间模型的解析解与数值解对比图；
图7示出了本说明书一实施例中层状各向异性粘弹性模型示意图；
图8a～图8d分别示出了本说明书一实施例中层状各向异性粘弹性模型的频率域位移场分量示意图；
图9示出了本说明书一些实施例中频率域粘弹性波正演模拟装置的结构框图；
图10示出了本说明书一些实施例计算机设备的结构框图。
【附图标记说明】
91、获取模块；
92、转换模块；
93、离散模块；
94、加载模块；
95、求解模块；
96、确定模块；
97、变换模块；
1102、计算机设备；
1104、处理器；
1106、存储器；
1108、驱动机构；
1110、输入/输出模块；
1112、输入设备；
1114、输出设备；
1116、呈现设备；
1118、图形用户接口；
1120、网络接口；
1122、通信链路；
1124、通信总线。
具体实施方式
为了使本技术领域的人员更好地理解本说明书中的技术方案，下面将结合本说明书实施例中的附图，对本说明书实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本说明书一部分实施例，而不是全部的实施例。基于本说明书中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都应当属于本说明书保护的范围。
参考图1所示，本说明书实施例的频率域粘弹性波正演模拟方法可以包括以下步骤：
S101、获取频率波数域一维粘弹性波动方程。
S102、将所述频率波数域一维粘弹性波动方程转换为等效积分弱形式。
S103、将频率波数域一维粘弹性波动方程的等效积分弱形式离散化为有限元控制方程。
S104、加载目标工区的三维点力源和吸收边界条件至所述有限元控制方程。
S105、对加载有所述三维点力源和所述吸收边界条件的有限元控制方程进行求解，获得所述目标工区的频率波数域波场。
S106、根据自适应波数采样策略确定所述频率波数域波场对应的波数集合。
S107、对所述波数集合中每个波数作为输入时对应的频率波数域波场，进行非等间隔二维傅里叶反变换，获得对应的频率域三维弹性波波场。
在本说明书实施例中，可以根据自适应波数采样策略确定频率波数域波场对应的波数集合，并对波数集合中每个波数作为输入时对应的频率波数域波场，进行非等间隔二维傅里叶反变换，获得对应的频率域三维弹性波波场；如此，不仅考虑了频率波数域粘弹性波场的振荡效应，并且辅之以自适应波数采样策略，从而实现三维点力源在任意一维各向异性粘弹性介质中的波场模拟；因此，与现有技术相比，本说明书实施例能够大大提升震源附近粘弹性波场的模拟精度。
在本说明一些实施例中，在频率域中，根据弹性波运动方程、本构方程以及应变 -位移关系可以表示为：
Figure BDA0002941127010000071
其中：σ＝(σxx，σyy，σzz，σyz，σxz，σxy)T
e＝(exx，eyy，ezz，eyz，exz，exy)T
Figure BDA0002941127010000072
L为偏微分算子矩阵；σ为应力向量；f为震源矢量；u为频率域位移矢量；e为应变向量；C为粘弹性介质的复刚度矩阵；x，y，z分别为空间三维方向；ρ为介质密度；ω为角频率；ρ为介质密度；ω为角频率，u为频率域位移矢量；e为应变向量； T为矩阵转置；在应力向量σ中，每个元素的第一个下标表示曲面的法向量方向，每个元素的第二个下标表示应力方向；例如，以σxy为例，第一个下标x表示曲面的法向量方向为x方向，第二个下标y表示应力方向为y方向。类似的，在应变向量e中，每个元素的第一个下标表示曲面的法向量方向，每个元素的第二个下标表示应变方向；例如，以exy为例，第一个下标x表示曲面的法向量方向为x方向，第二个下标y表示应变方向为y方向。
对于沿深度方向的一维介质，C和ρ仅与z有关。因此，将频率波数域一维粘弹性波动方程关于x方向和y方向进行一维空间傅里叶变换，可以得到如下所示的频率波数域一维粘弹性波动方程：
Figure BDA0002941127010000081
其中，
Figure BDA0002941127010000082
分别为σ,e,f,L,u关于y方向的二维傅里叶变换。
在本说明书实施例中，将所述频率波数域一维粘弹性波动方程转换为等效积分弱形式可以方便后续进行离散化处理。在本说明一些实施例中，所述将所述频率波数域一维粘弹性波动方程转换为等效积分弱形式，可以包括：
1)、确定
Figure BDA0002941127010000083
的余量形式，即余量
Figure BDA0002941127010000084
2)、根据Galerkin加权余量法，取
Figure BDA0002941127010000085
为权函数，R满足积分方程：
Figure BDA0002941127010000086
3)、将
Figure BDA0002941127010000087
代入
Figure BDA0002941127010000088
并结合格林公式，可以得到频率波数域一维粘弹性波动方程的等效积分弱形式为：
Figure BDA0002941127010000089
其中，Γ为计算区域；
Figure BDA00029411270100000810
Figure BDA00029411270100000814
为
Figure BDA00029411270100000811
的变分； i为虚数单位；kx为x方向的波数；ky为y方向的波数。
在本说明一些实施例中，所述将频率波数域一维粘弹性波动方程的等效积分弱形式离散化为有限元控制方程，可以包括：
1)、将深度方向的计算区域Γ剖分为NE个线单元，(这里以线单元分析为例)，则
Figure BDA00029411270100000812
可以写成如下离散形式：
Figure BDA00029411270100000813
在每个线单元内，
Figure BDA0002941127010000091
可以近似为关于z的二次插值函数：
Figure BDA0002941127010000092
其中，
Figure BDA0002941127010000093
为第e个线单元中第n个顶点处的位移矢量；Nn(n＝1,2,3)为关于z的一维形函数；I为3×3的单位阵。
2)、将上述
Figure BDA0002941127010000094
代入上述的离散形式：
Figure BDA0002941127010000095
可以得到有限元控制方程为：
Figure BDA0002941127010000096
其中，K为全局刚度矩阵；M为全局质量矩阵；F为全局震源矢量；u为频率域位移矢量；K,M,F的单元矩阵分别为
Figure BDA0002941127010000097
Figure BDA0002941127010000098
N为形函数矩阵；Γe为与第e个离散单元对应的计算区域。
在本说明书一些实施例中，三维点力源的加载可以选择雷克子波作为震源函数。频率域雷克子波函数可以表示为
Figure BDA0002941127010000099
其中，fm为雷克子波的主频； f为雷克子波的频率。于是，在某一指定位置(0,0,zs)激发的点力源(一个震源矢量) 可以表示为：
f＝AR(f)δ(x)δ(y)δ(z-zs)ns
其中，
Figure BDA00029411270100000910
为f的方向余弦；
Figure BDA00029411270100000911
为f与z轴夹角；θ为f 在xy平面的投影与x轴的夹角；δ为狄拉克函数；A为震源幅值；zs为震源的z轴坐标。
将f＝AR(f)δ(x)δ(y)δ(z-zs)ns关于x方向和y方向做二维空间傅里叶变换，可以得到
Figure BDA00029411270100000912
将
Figure BDA00029411270100000913
代入上述的
Figure BDA00029411270100000914
可以得到：
Figure BDA00029411270100000915
即在加载三维点力源后，全局震源矢量F的单元矩阵变为：
Figure BDA0002941127010000101
在本说明书一些实施例中，为了消除人工边界引起的虚假反射，可以在计算区域边界处施加SRM吸收边界条件(如图2所示)。在此吸收边界条件中，阻尼矩阵D可以由质量矩阵M和刚度矩阵K的线性组合表示为D＝DMM。
式中，DM为质量阻尼系数，且DM＝DMmaxZ(z′)3；DMmax为质量阻尼系数极大值，且等于入射波的角频率；z′为SRM区域的计算点到内边界的距离；
Figure BDA0002941127010000102
且其取值范围为[0,1]，即在SRM内边界上取值为0，外边界取值为1。
此外，为了增强SRM边界的吸收效果，还需要减小吸收层介质的刚度，即粘弹性介质的复刚度矩阵需满足：
Figure BDA0002941127010000103
其中，C(z′)为SRM区域的计算点到内边界的距离为x′时对应的粘弹性介质的复刚度矩阵，α(z′)＝αmaxZ(z′)p；
Figure BDA0002941127010000104
P为常数，一般取值为3；kinc为入射波的波数；LSRM为SRM吸收层的厚度；rC为SRM吸收层内边界与外边界弹性刚度的比值；C0为SRM区域内边界的刚度矩阵；Re[C0]为C0的实部；e为自然常数。
在加入阻尼后，有限元控制方程可以表示为：
Figure BDA0002941127010000105
如此，加载有所述三维点力源和所述吸收边界条件的有限元控制方程可以表示为：
Figure BDA0002941127010000106
其中，全局震源矢量F的单元矩阵变为：
Figure BDA0002941127010000107
在本说明书一些实施例中，可以采用直接解法等求解
Figure BDA0002941127010000108
从而可以获得频率波数域波场(即频率波数域的粘弹性波波场)。
在获得频率波数域波场
Figure BDA0002941127010000109
后，还需将不同波数(kx,ky)下计算的频率波数域波场
Figure BDA00029411270100001010
通过非等间隔二维反傅里叶反变换到空间域中，即
Figure BDA00029411270100001011
根据频率波数域波场
Figure BDA00029411270100001012
的变化特征(如图3a～图3d所示；在图3a和图3c中，
Figure BDA00029411270100001013
为波场
Figure BDA00029411270100001014
的实部；在图3b和图3d中，
Figure BDA00029411270100001015
为波场
Figure BDA00029411270100001016
的虚部)，
Figure BDA00029411270100001017
存在一个截断波数kc，即
Figure BDA00029411270100001018
其中，Vmi为最小速度，在二维波数范围
Figure BDA00029411270100001019
内，
Figure BDA00029411270100001020
具有振荡变化特征，且在距离曲线
Figure BDA00029411270100001021
越近的波数点上，其振荡越显著；而当
Figure BDA00029411270100001022
时，
Figure BDA00029411270100001023
则以一定的速度单调衰减至零，且衰减速度随炮检距的增大而增大。其中，Ω表示有效二维波数区域。
基于此，在本说明书一些实施例中，利用三角单元对二维波数空间进行离散化(如图4所示)可以采用如下波数采样策略确定区间
Figure BDA0002941127010000111
内的三角单元的平均尺寸：
Figure BDA0002941127010000112
式中，Δk为三角单元的平均尺寸；Nk为波数k的采样点数量；max为最大值函数；int为取整函数；rmax为最大炮检距；λmax为最大波长。
另外，为提高震源附近的计算精度，可以设置最大波数为kymax＝10kc，并在波数区间
Figure BDA0002941127010000113
内采用非规则网格，网格密度由
Figure BDA0002941127010000114
至
Figure BDA0002941127010000115
逐渐降低；并设定此区域三角单元的最大尺度设置为(Δk)max＝0.9kc。
在本说明书一些实施例中，基于形函数构建非等间隔二维傅里叶变换算法，可以获得频率域三维弹性波波场。为了方便处理，
Figure BDA0002941127010000116
可以重新写成如下形式：
Figure BDA0002941127010000117
其中，C＝1/(2π)且ξ＝x,ζ＝y,η＝kx,μ＝ky,γ＝ix,λ＝iy；u为位移矢量的任一分量，
Figure BDA0002941127010000118
和u互为傅里叶变换；
Figure BDA0002941127010000119
为虚数单位；π为圆周率。
在常规技术中，通常采用二维快速傅里叶变换计算。尽管这种做法计算效率较高，但由于二维快速傅里叶变换技术仅适用于矩形网格，因此无法准确描述频率-波数域 位移场的变化特征，从而导致数值精度较低。为了改善二维傅里叶变换的精度，这里 采用三角网格对波场进行精确刻画，具体步骤如下：
将二维波数区域
Figure BDA00029411270100001111
离散为Ne个三角单元，得到非等间隔二维傅里叶变换离散公式为：
Figure BDA00029411270100001112
或
Figure BDA00029411270100001113
式中，Ie为第e个三角单元的单元积分，具体表达形式为
Figure BDA0002941127010000121
为了更好地描述频率波数域位移场的振荡特性，这里仅对核函数
Figure BDA0002941127010000122
进行数值近似，从而保留振荡指数项exp(γη+λμ)参与单元积分。一般而言，对于在大部分情况下，二次插值便能满足精度需求，故有：
Figure BDA0002941127010000123
式中，m＝1,2,…,6分别为第e个三角单元中第m个离散节点编号；
Figure BDA0002941127010000124
为第e个三角单元中第m个离散节点的位移分量；
Figure BDA0002941127010000125
为二维形函数，即：
其中，
Figure BDA0002941127010000126
将
Figure BDA0002941127010000127
代入
Figure BDA0002941127010000128
中，可以得到：
Figure BDA0002941127010000129
其中，
Figure BDA00029411270100001210
为第e个三角单元中第m个离散节点的权系数。
考虑到三角单元形状的任意性，直接由
Figure BDA00029411270100001211
来求解积分权系数
Figure BDA00029411270100001212
比较困难。为此，可以利用如下方法将积分变量(η,μ)替换为(L1,L2)，从而将面积分公式
Figure BDA00029411270100001213
转为积分变量相对独立的二重积分。具体的：
首先，根据关系式
Figure BDA00029411270100001214
得到
Figure BDA00029411270100001215
然后，利用Jacobi变换，得到变量替换公式为
Figure BDA0002941127010000131
将变量替换公式代入积分权系数的
Figure BDA0002941127010000132
中，便得到
Figure BDA0002941127010000133
的简单计算公式：
Figure BDA0002941127010000134
其中：
χ1＝b2γ-a2λ
χ2＝a1λ-b1γ
χ3＝η3γ+μ3λ
在此基础上，将
Figure BDA0002941127010000135
中，便能获得频率域位移场的计算公式：
Figure BDA0002941127010000136
其中，u(ξ,ζ)为(ξ,ζ)下的频率域三维弹性波波场；n为二维波数域中第n个离散节点的编号，n＝1,2,…,Nd；Nd为离散节点总数；Wn为第n个离散节点处的积分权系数；Ne为二维波数域中划分出的离散单元数量；C＝1/(2π)；m＝1,2,…,6分别为第e个离散单元的6个节点编号；
Figure BDA0002941127010000137
为第e个三角单元中第m个离散节点的权系数；
Figure BDA0002941127010000138
为第e个三角单元中第m个离散节点的位移分量，
Figure BDA0002941127010000139
为第n个离散节点的位移分量。
获得的频率域三维弹性波波场，可以用于研究地震波在目标地球介质中的传播规律，从而为后续对目标地球介质的油气资源进行勘探开发提供参考依据。
图5示出了本说明书一示例性实施例中的各向异性全空间介质模型。基于此模型，利用上述的频率域粘弹性波正演模拟方法，可以获得频率波数域波场；然后在自适应波数采样策略下，对不同波数下的频率波数域波场进行非等间隔二维空间反傅里叶变换，即可以获得对应的频率域三维弹性波场。图6中示出了各向异性全空间介质模型的数值解(图6中的虚线)和解析解(图6中的实线)。在图6的(a)中，Re(uxx)为波场uxx的实部；在图6的(b)中，Im(uxx)为波场uxx的虚部；在图6的(c)中，Re(uxz) 为波场uxz的实部；在图6的(d)的中，Im(uxz)为波场uxz的虚部。从图6中可以看出，基于本说明书实施例的频率域粘弹性波正演模拟方法所得到的数值解，与对应的解析解吻合度很高，从而验证了本说明书实施例的频率域粘弹性波正演模拟方法具有很高的精度。此外，基于本说明书实施例的频率域粘弹性波正演模拟方法，对于如图7 所示的层状各向异性衰减模型，其对应的频率域粘弹性波正演模拟结果可以如图8a～图8d所示。在图8a中，Re(uxx)为波场uxx的实部；在图8b中，Im(uxx)为波场uxx的虚部；在图8c中，Re(uxz)为波场uxz的实部；在图8d中，Im(uxz)为波场uxz的虚部。
虽然上文描述的过程流程包括以特定顺序出现的多个操作，但是，应当清楚了解，这些过程可以包括更多或更少的操作，这些操作可以顺序执行或并行执行(例如使用并行处理器或多线程环境)。
参考图9所示，与上述的频率域粘弹性波正演模拟方法对应，本说明书一些实施例中还提供了频率域粘弹性波正演模拟装置，其可以包括：
获取模块91，可以用于获取频率波数域一维粘弹性波动方程；
转换模块92，可以用于将所述频率波数域一维粘弹性波动方程转换为等效积分弱形式；
离散模块93，可以用于将频率波数域一维粘弹性波动方程的等效积分弱形式离散化为有限元控制方程；
加载模块94，可以用于加载目标工区的三维点力源和吸收边界条件至所述有限元控制方程；
求解模块95，可以用于对加载有所述三维点力源和所述吸收边界条件的有限元控制方程进行求解，获得所述目标工区的频率波数域波场；
确定模块96，可以用于根据自适应波数采样策略确定所述频率波数域波场对应的波数集合；
变换模块97，可以用于对所述波数集合中每个波数作为输入时对应的频率波数域波场，进行非等间隔二维傅里叶反变换，获得对应的频率域三维弹性波波场。
为了描述的方便，描述以上装置时以功能分为各种单元分别描述。当然，在实施本说明书时可以把各单元的功能在同一个或多个软件和/或硬件中实现。
本说明书的实施例还提供一种计算机设备。如图10所示，在本说明书一些实施例中，所述计算机设备1102可以包括一个或多个处理器1104，诸如一个或多个中央处理单元(CPU)或图形处理器(GPU)，每个处理单元可以实现一个或多个硬件线程。计算机设备1102还可以包括任何存储器1106，其用于存储诸如代码、设置、数据等之类的任何种类的信息，一具体实施方式中，存储器1106上并可在处理器1104上运行的计算机程序，所述计算机程序被所述处理器1104运行时，可以执行根据上述方法的指令。非限制性的，比如，存储器1106可以包括以下任一项或多种组合：任何类型的RAM，任何类型的ROM，闪存设备，硬盘，光盘等。更一般地，任何存储器都可以使用任何技术来存储信息。进一步地，任何存储器可以提供信息的易失性或非易失性保留。进一步地，任何存储器可以为计算机设备1102的固定或可移除部件。在一种情况下，当处理器1104执行被存储在任何存储器或存储器的组合中的相关联的指令时，计算机设备1102可以执行相关联指令的任一操作。计算机设备1102还包括用于与任何存储器交互的一个或多个驱动机构1108，诸如硬盘驱动机构、光盘驱动机构等。
计算机设备1102还可以包括输入/输出模块1110(I/O)，其用于接收各种输入(经由输入设备1112)和用于提供各种输出(经由输出设备1114)。一个具体输出机构可以包括呈现设备1116和相关联的图形用户接口1118(GUI)。在其他实施例中，还可以不包括输入/输出模块1110(I/O)、输入设备1112以及输出设备1114，仅作为网络中的一台计算机设备。计算机设备1102还可以包括一个或多个网络接口1120，其用于经由一个或多个通信链路1122与其他设备交换数据。一个或多个通信总线1124将上文所描述的部件耦合在一起。
通信链路1122可以以任何方式实现，例如，通过局域网、广域网(例如，因特网)、点对点连接等、或其任何组合。通信链路1122可以包括由任何协议或协议组合支配的硬连线链路、无线链路、路由器、网关功能、名称服务器等的任何组合。
本申请是参照本说明书一些实施例的方法、设备(***)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理器的处理器以产生一个机器，使得通过计算机或其他可编程数据处理器的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理器以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理器上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
在一个典型的配置中，计算机设备包括一个或多个处理器(CPU)、输入/输出接口、网络接口和内存。
内存可能包括计算机可读介质中的非永久性存储器，随机存取存储器(RAM)和/或非易失性内存等形式，如只读存储器(ROM)或闪存(flash RAM)。内存是计算机可读介质的示例。
计算机可读介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机的存储介质的例子包括，但不限于相变内存(PRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器(DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体或其他内存技术、只读光盘只读存储器(CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带，磁盘存储或其他磁性存储设备或任何其他非传输介质，可用于存储可以被计算机设备访问的信息。按照本文中的界定，计算机可读介质不包括暂存电脑可读媒体 (transitory media)，如调制的数据信号和载波。
本领域技术人员应明白，本说明书的实施例可提供为方法、***或计算机程序产品。因此，本说明书实施例可采用完全硬件实施例、完全软件实施例或结合软件和硬件方面的实施例的形式。而且，本说明书实施例可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本说明书实施例可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构等等。也可以在分布式计算环境中实践本说明书实施例，在这些分布式计算环境中，由通过通信网络而被连接的远程处理器来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。
本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于***实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本说明书实施例的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。
以上所述仅为本申请的实施例而已，并不用于限制本申请。对于本领域技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本申请的权利要求范围之内。

Claims (10)

1.一种频率域粘弹性波正演模拟方法，其特征在于，包括：
获取频率波数域一维粘弹性波动方程；
将所述频率波数域一维粘弹性波动方程转换为等效积分弱形式；
将频率波数域一维粘弹性波动方程的等效积分弱形式离散化为有限元控制方程；
加载目标工区的三维点力源和吸收边界条件至所述有限元控制方程；
对加载有所述三维点力源和所述吸收边界条件的有限元控制方程进行求解，获得所述目标工区的频率波数域波场；
根据自适应波数采样策略确定所述频率波数域波场对应的波数集合；
对所述波数集合中每个波数作为输入时对应的频率波数域波场，进行非等间隔二维傅里叶反变换，获得对应的频率域三维弹性波波场。
2.如权利要求1所述的频率域粘弹性波正演模拟方法，其特征在于，所述频率波数域一维粘弹性波动方程，包括：
Figure FDA0002941127000000011
其中，
Figure FDA0002941127000000012
分别为σ,e,f,L,u关于x和y方向的二维傅里叶变换，L为偏微分算子矩阵，且
Figure FDA0002941127000000013
σ为应力向量；x,y,z分别为空间三维方向；f为震源矢量；ρ为介质密度；ω为角频率，u为频率域位移矢量；e为应变向量；C为粘弹性介质的复刚度矩阵；T为矩阵转置；
Figure FDA0002941127000000014
分别为关于x,y,z方向的偏导。
3.如权利要求2所述的频率域粘弹性波正演模拟方法，其特征在于，所述频率波数域一维粘弹性波动方程的等效积分弱形式为：
Figure FDA0002941127000000015
其中，Γ为计算区域；
Figure FDA0002941127000000021
Figure FDA0002941127000000022
为
Figure FDA0002941127000000023
的变分；i为虚数单位；kx为x方向的波数；ky为y方向的波数。
4.如权利要求3所述的频率域粘弹性波正演模拟方法，其特征在于，所述有限元控制方程包括：
Figure FDA0002941127000000024
其中，K为全局刚度矩阵；M为全局质量矩阵；F为全局震源矢量；K,M,F的单元矩阵分别为
Figure FDA0002941127000000025
Figure FDA0002941127000000026
N为形函数矩阵；Γe为与第e个离散单元对应的计算区域。
5.如权利要求4所述的频率域粘弹性波正演模拟方法，其特征在于，所述加载有所述三维点力源和所述吸收边界条件的有限元控制方程，包括：
Figure FDA0002941127000000027
其中，F的单元矩阵为
Figure FDA0002941127000000028
Figure FDA0002941127000000029
为f的方向余弦，
Figure FDA00029411270000000210
为f与z轴夹角，θ为f在xy平面的投影与x轴的夹角；δ为狄拉克函数；A为震源幅值；R(f)为频率域雷克子波函数；f为雷克子波的频率；zs为震源的z轴坐标；D为阻尼矩阵，且D＝DMM；DM为质量阻尼系数，且DM＝DMmaxZ(z′)3；DMmax为质量阻尼系数极大值，且等于入射波的角频率；z′为刚性弱化SRM区域的计算点到内边界的距离；
Figure FDA00029411270000000211
C满足
Figure FDA00029411270000000212
C(z′)为刚性弱化区域的计算点到内边界的距离为z′时对应的粘弹性介质的复刚度矩阵；α(z′)＝αmaxZ(z′)p；
Figure FDA00029411270000000213
p为常数；kinc为入射波的波数；LSRM为SRM吸收层的厚度；rC为SRM吸收层内边界与外边界弹性刚度的比值；C0为SRM区域内边界的刚度矩阵；Re[C0]为C0的实部；e为自然常数。
6.如权利要求5所述的频率域粘弹性波正演模拟方法，其特征在于，所述自适应波数采样策略包括：
将波数区间
Figure FDA0002941127000000031
内的三角单元的平均尺寸设置为：
Figure FDA0002941127000000032
最大波数设置为：kmax＝10kc；
在波数区间
Figure FDA0002941127000000033
内采用非规则网格，网格密度由
Figure FDA0002941127000000034
至
Figure FDA0002941127000000035
逐渐降低；
三角单元的最大尺寸设置为：(Δk)max＝0.9kc；
其中，Δk为三角单元的平均尺寸；Nk为波数k的采样点数量；max为最大值函数；int为取整函数；rmax为最大炮检距；λmax为最大波长；
Figure FDA0002941127000000036
Vmin为最小速度。
7.如权利要求1所述的频率域粘弹性波正演模拟方法，其特征在于，所述频率域三维弹性波波场包括：
Figure FDA0002941127000000037
其中，u(ξ,ζ)为(ξ,ζ)下的频率域三维弹性波波场；n为二维波数域中第n个离散节点的编号，n＝1,2,…,Nd；Nd为离散节点总数；Wn为第n个离散节点处的积分权系数；Ne为二维波数域中划分出的离散单元数量；C＝1/(2π)；m＝1,2,…,6分别为第e个离散单元的6个节点编号；
Figure FDA0002941127000000038
为第e个三角单元中第m个离散节点的权系数；
Figure FDA0002941127000000039
为第e个三角单元中第m个离散节点的位移分量，
Figure FDA00029411270000000310
为第n个离散节点的位移分量。
8.一种频率域粘弹性波正演模拟装置，其特征在于，包括：
获取模块，用于获取频率波数域一维粘弹性波动方程；
转换模块，用于将所述频率波数域一维粘弹性波动方程转换为等效积分弱形式；
离散模块，用于将频率波数域一维粘弹性波动方程的等效积分弱形式离散化为有限元控制方程；
加载模块，用于加载目标工区的三维点力源和吸收边界条件至所述有限元控制方程；
求解模块，用于对加载有所述三维点力源和所述吸收边界条件的有限元控制方程进行求解，获得所述目标工区的频率波数域波场；
确定模块，用于根据自适应波数采样策略确定所述频率波数域波场对应的波数集合；
变换模块，用于对所述波数集合中每个波数作为输入时对应的频率波数域波场，进行非等间隔二维傅里叶反变换，获得对应的频率域三维弹性波波场。
9.一种计算机设备，包括存储器、处理器、以及存储在所述存储器上的计算机程序，其特征在于，所述计算机程序被所述处理器运行时，执行根据权利要求1-7任意一项所述方法的指令。
10.一种计算机存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被计算机设备的处理器运行时，执行根据权利要求1-7任意一项所述方法的指令。
CN202110177528.6A 2021-02-07 2021-02-07 频率域粘弹性波正演模拟方法、装置、设备及存储介质 Pending CN112906267A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110177528.6A CN112906267A (zh) 2021-02-07 2021-02-07 频率域粘弹性波正演模拟方法、装置、设备及存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110177528.6A CN112906267A (zh) 2021-02-07 2021-02-07 频率域粘弹性波正演模拟方法、装置、设备及存储介质

Publications (1)

Publication Number Publication Date
CN112906267A true CN112906267A (zh) 2021-06-04

Family

ID=76123062

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110177528.6A Pending CN112906267A (zh) 2021-02-07 2021-02-07 频率域粘弹性波正演模拟方法、装置、设备及存储介质

Country Status (1)

Country Link
CN (1) CN112906267A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113591423A (zh) * 2021-09-30 2021-11-02 北京智芯仿真科技有限公司 有损耗有频散介质下的集成电路全波电磁仿真方法及***

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113591423A (zh) * 2021-09-30 2021-11-02 北京智芯仿真科技有限公司 有损耗有频散介质下的集成电路全波电磁仿真方法及***

Similar Documents

Publication Publication Date Title
Seriani et al. Spectral element method for acoustic wave simulation in heterogeneous media
De La Puente et al. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes-IV. Anisotropy
Lombard et al. Free and smooth boundaries in 2-D finite-difference schemes for transient elastic waves
Packo et al. Numerical simulations of elastic wave propagation using graphical processing units—comparative study of high-performance computing capabilities
Withers et al. Memory-efficient simulation of frequency-dependent Q
CN107798156B (zh) 一种频率域2.5维粘弹性波数值模拟方法及装置
CN109143339B (zh) 基于横波应力不变量的弹性逆时偏移成像方法及装置
Wu Efficient modeling of gravity fields caused by sources with arbitrary geometry and arbitrary density distribution
Kataoka et al. A parallel iterative partitioned coupling analysis system for large-scale acoustic fluid–structure interactions
Franczyk Using the Morris sensitivity analysis method to assess the importance of input variables on time-reversal imaging of seismic sources
CN114114403B (zh) 一种基于分数阶拉氏算子的各向异性衰减介质模拟方法
Bamer et al. A Newmark space-time formulation in structural dynamics
CN112906267A (zh) 频率域粘弹性波正演模拟方法、装置、设备及存储介质
CN108828659B (zh) 基于傅里叶有限差分低秩分解的地震波场延拓方法及装置
CN114139335A (zh) 基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法
Soare et al. Calibration and fast evaluation algorithms for homogeneous orthotropic polynomial yield functions
Petrov et al. Higher-order grid-characteristic schemes for the acoustic system
US10454713B2 (en) Domain decomposition using a multi-dimensional spacepartitioning tree
AlSalem et al. Embedded boundary methods for modeling 3D finite-difference Laplace-Fourier domain acoustic-wave equation with free-surface topography
Wang et al. A stable node-based smoothed finite element method with transparent boundary conditions for the elastic wave scattering by obstacles
CN112966412A (zh) 频率域2.5维地震波正演模拟方法、装置及设备
CN113221392B (zh) 一种非均匀粘声波在无限域内传播模型的构建方法
Kim et al. Two-dimensional grid optimization for sedimentation velocity analysis in the analytical ultracentrifuge
Ratcliffe et al. Measuring rotational degrees of freedom using a laser doppler vibrometer
Gordon et al. A compact three‐dimensional fourth‐order scheme for elasticity using the first‐order formulation

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