CN116882214B - 基于dfl粘弹性方程的瑞雷波数值模拟方法及系统 - Google Patents

基于dfl粘弹性方程的瑞雷波数值模拟方法及系统 Download PDF

Info

Publication number
CN116882214B
CN116882214B CN202311146141.XA CN202311146141A CN116882214B CN 116882214 B CN116882214 B CN 116882214B CN 202311146141 A CN202311146141 A CN 202311146141A CN 116882214 B CN116882214 B CN 116882214B
Authority
CN
China
Prior art keywords
representing
wave
dfl
cpml
equation
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
CN202311146141.XA
Other languages
English (en)
Other versions
CN116882214A (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.)
Sanya Offshore Oil And Gas Research Institute Of Northeast Petroleum University
Original Assignee
Sanya Offshore Oil And Gas Research Institute Of Northeast Petroleum University
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 Sanya Offshore Oil And Gas Research Institute Of Northeast Petroleum University filed Critical Sanya Offshore Oil And Gas Research Institute Of Northeast Petroleum University
Priority to CN202311146141.XA priority Critical patent/CN116882214B/zh
Publication of CN116882214A publication Critical patent/CN116882214A/zh
Application granted granted Critical
Publication of CN116882214B publication Critical patent/CN116882214B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • 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)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Geometry (AREA)
  • Operations Research (AREA)
  • Evolutionary Computation (AREA)
  • Algebra (AREA)
  • Computer Hardware Design (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本申请公开了基于DFL粘弹性方程的瑞雷波数值模拟方法及系统，属于地震波模拟技术领域，方法包括以下步骤：根据地震波在衰减介质中的传播特性，建立DFL粘弹性方程；根据面波产生条件和消除人工边界反射方法，设置真空自由表面条件和CPML边界条件；以所述DFL粘弹性方程、真空形式和CPML作为模型，采用SGPS法数值实现瑞雷波正演模拟。本申请以DFL粘弹性方程描述体波的传播，真空形式作为自由表面边界条件，CPML为吸收边界条件，采用SGPS法对浅层各向同性衰减介质进行了包括瑞雷波和体波在内的全波场模拟，该方法既继承了真空形式的简单易实现特点，又避免了使用真空形式时困扰基于GSLS方程的数值频散问题。

Description

基于DFL粘弹性方程的瑞雷波数值模拟方法及系统
技术领域
本申请属于地震波模拟技术领域，具体涉及基于DFL粘弹性方程的瑞雷波数值模拟方法及系统。
背景技术
在层状介质中，瑞雷波具有相速度随频率变化而变化的频散特性，这使得瑞雷波在浅地表结构探测中具有广泛的应用，如反演地下介质参数和识别地层厚度。此外，若在全波形反演中，不考虑瑞雷波信息会导致数据匹配性差，最终降低速度反演结果的可靠性。另一方面，瑞雷波频散曲线对横波速度和地层厚度比较敏感，通过多道面波分析方法可以反演横波速度剖面和识别地层厚度。因此，准确模拟瑞雷波对理解地震波传播机理和建立可靠的地下模型至关重要。
由于瑞雷波产生条件严格，而真实地下介质的粘滞性对瑞雷波的影响实际上比体波更为复杂，因此需要选取合适的粘滞波动方程。在勘探地球物理领域中，通常应用带有自由边界条件的广义标准线性体(GSLS)模型描述体波和瑞雷波在地下介质中的传播。
但在现有的技术中，GSLS应用有限差分法求解时，具有便于时域并行处理的优点，但该模型也存在一定的问题。首先，GSLS内包含了记忆变量，从而导致大量额外的计算时间和内存，特别是在3D弹性问题中；其次，品质因子是一组松弛时间的函数，两者间的对应关系需要通过优化方法确定，这无形之中又增加了计算成本；GSLS中控制振幅衰减和速度频散项是耦合在一起的，不利于振幅补偿和相位校正。真空法具有实现过程简单，计算效率较高的特点，但该方法只适用于二阶精度的差分格式，在四阶差分格式中，计算易发生不稳定现象。因此，带有真空自由边界条件的GSLS并不适用于高阶差分格式。
发明内容
本申请旨在解决现有技术的不足，提出一种适用于DFL粘弹性方程的瑞雷波数值模拟方法及系统。采用DFL粘弹性方程与真空自由表面条件相结合的方法对瑞雷波进行高精度数值模拟。
为实现上述目的，本申请提供了如下方案：
基于DFL粘弹性方程的瑞雷波数值模拟方法，包括以下步骤：
根据地震波在衰减介质中的传播特性，建立DFL粘弹性方程；
根据面波产生条件和消除人工边界反射方法，设置真空自由表面条件和CPML边界条件；
以所述DFL粘弹性方程、真空形式和CPML作为模型，采用SGPS法数值实现瑞雷波正演模拟。
优选的，所述DFL粘弹性方程表示为：
；
其中，
；
和/>表示x和z方向的速度分量；/>和/>表示应力分量；/>表示介质密度；表示x和z方向的震源项；/>或/>，表示P波或S波；/>和/>表示传播速度和衰减强度；/>和/>表示在参考频率/>处的P波和S波速度；/>和/>表示P波和S波的品质因子；t表示时刻；/>表示拉普拉斯算子。
优选的，设置所述真空自由表面条件的方法包括：
将自由表面设置在通过正应力和速度水平分量的采样位置；
在所述自由表面上部添加真空层，并将所述真空层中的密度设置为零。
优选的，设置所述CPML边界条件的方法包括：
定义微分算子和复频移拉伸函数分别为：
式中，表示空间方向，取/>控制/>方向的衰减，取表示CPML内的点到CPML内边界的距离；/>表示CPML区域的厚度；/>表示正整数；/>表示理论反射系数，取/>控制表面波的吸收，取的正实数；/>控制低频分量的吸收，取/>为子波的主频；i表示虚数；/>表示圆角频率；vmax表示最大速度；
对作傅里叶逆变换得到：
其中，是/>的傅里叶逆变换，其表达式为：
其中，为单位脉冲函数，/>为Heaviside阶跃函数；
令：
则表示为：
其中第二项的卷积，在第个时间步写为积分形式：
由于是定义在半个时间步上，故写为如下离散形式：
其中，，表示积分变量；
根据递归方法，将离散形式写为如下递推形式：
即，CPML应用于所述DFL粘弹性方程表示为：
其中，；，控制j方向的衰减；/>表示CPML内的点到CPML内边界的距离；/>表示CPML区域的厚度；/>表示正整数；/>表示理论反射系数；/>，控制表面波的吸收；/>表示正实数；/>，控制低频分量的吸收；/>表示子波的主频；/>表示空间方向。
优选的，所述实现瑞雷波正演模拟的过程包括：
计算空间导数和分数阶拉普拉斯算子；
基于所述空间导数和所述分数阶拉普拉斯算子更新质点的速度和应力，以及吸收层质点的速度和应力，得到所述模拟结果。
优选的，所述空间导数的计算方法为：
采用交错伪谱法计算所述空间导数：
；
其中，和/>分别表示一维傅里叶正变换和逆变换，上标/>表示空间相位左移/右移；
，上标/>表示当前时刻和前一时刻。
优选的，所述分数阶拉普拉斯算子的计算方法为：
利用二维傅里叶变换计算分数阶拉普拉斯算子：
；
时间导数采用差分法计算：
；
其中，和/>分别表示二维傅里叶正变换和逆变换；/>和/>表示空间方向；/>表示时间步长。
优选的，所述质点的速度和应力更新公式为：
；
所述吸收层质点的速度和应力更新公式为：
；
其中，分别表示x和z方向的震源项。
本申请还提供了基于DFL粘弹性方程的瑞雷波数值模拟系统，包括：方程构建模块、边界条件设置模块和数值实现模块；
所述方程构建模块用于根据地震波在衰减介质中的传播特性，建立DFL粘弹性方程；
所述边界条件设置模块用于根据面波产生条件和消除人工边界反射方法，设置真空自由表面条件和CPML边界条件；
所述数值实现模块用于以DFL粘弹性方程、真空形式和CPML作为模型，采用SGPS法数值实现瑞雷波正演模拟。
与现有技术相比，本申请的有益效果为：
本申请以DFL粘弹性方程描述体波在地下介质中的传播规律，真空形式作为自由表面边界条件，卷积完美匹配层吸收边界反射，交错伪谱法数值模拟衰减介质中瑞雷波的传播。该方法既继承了真空形式的简单易实现特点，又避免了使用真空形式时困扰基于GSLS方程(采用有限差分法求解)的数值色散问题。
附图说明
为了更清楚地说明本申请的技术方案，下面对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本申请实施例的方法流程示意图；
图2为本申请实施例的波场和介质参数在交错网格上的分布图；
图3为本申请实施例的均匀半空间模型偏移距600m的解析解与数值解对比，(a)为水平分量对比，(b)为垂直分量对比；
图4为本申请实施例的均匀半空间模型在t=1s的波场快照，(a)和(b)空间步长2.5m和时间步长0.5ms的SGFD法模拟的和/>分量；(c)和(d)参数同SGFD法的SGPS模拟的/>和/>分量；(e)和(f)减小参数的SGFD法模拟的/>和/>分量；
图5为本申请实施例的自由表面的局部放大图，(a)和(b)空间步长2.5m和时间步长0.5ms的SGFD法模拟的和/>分量；(c)和(d)参数同SGFD法的SGPS模拟的/>和/>分量；(e)和(f)减小参数的SGFD法模拟的/>和/>分量；
图6为本申请实施例的均匀半空间模型偏移距600m的波形曲线对比，(a)为分量对比，(b)为/>分量对比；
图7为本申请实施例的均匀半空间模型的波场快照，(a)为弹性/>分量，(b)为弹性/>分量，(c)为粘弹性/>分量，(d)为粘弹性/>分量；
图8为本申请实施例的粘弹性介质均匀半空间模型炮集记录(a)和频散能量图(b)；
图9为本申请实施例的粘弹性两层模型的速度场快照，(a)为分量，(b)为/>分量；
图10为本申请实施例的粘弹性两层模型的炮集记录(a)和频散能量图(b)；
图11为本申请实施例的两层模型(a)和t=0.35s的波场快照(b)；
图12为本申请实施例的不同Q值的波形曲线比较，(a)为和情况的比较图，(b)为/>和/>情况下的比较图；
图13为本申请实施例的模拟时长为10s的能量衰减曲线；
图14为本申请实施例的Marmousi模型示意图，(a)为P波速度模型，(b)为S波速度模型，(c)为模型，(d)为/>模型；
图15为本申请实施例的分量的波场快照和炮集记录，(a)为本实施例方案的波场快照，(b)为本实施例方案的炮集记录，(c)为SIM模拟的波场快照，(d)为SIM模拟的炮集记录，(e)为(a)和(c)之间的差值，(f)为(b)和(d)之间的差值；
图16为本申请实施例的无面波和有面波的波场快照和炮集记录，(a)为本实施例方案的波场快照，(b)为本实施例方案的炮集记录，(c)为SIM模拟的波场快照，(d)为SIM模拟的炮集记录，(e)为(a)和(c)之间的差值，(f)为(b)和(d)之间的差值；
图17为本申请实施例的系统结构示意图。
具体实施方式
下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。
为使本申请的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本申请作进一步详细的说明。
实施例一
在本实施例中，如图1所示，基于DFL粘弹性方程的瑞雷波数值模拟方法，包括以下步骤：
S1.根据地震波在衰减介质中的传播特性，建立DFL粘弹性方程。
Xing和Zhu在Kjartansson的常Q频散关系的基础上提出粘性介质声波的新频散关系：
其中是参考频率/>处对应的参考相速度；/>，新频散关系中波数/>的指数均为常数，不再是空间位置的函数。
由于包含P波和S波分量的弹性波也满足上述近似的频散关系，将上述频散关系进行推广，得到粘性介质弹性波的频散关系：
其中表示纵波或横波，取/>。
由此得到波数域的复模量：
将复模量代入弹性波的应力-应变关系：
其中和/>为波数域中的应力张量，/>和/>为波数域中的应变张量，i、j、k代表不同的空间方向。
将上式代入耦合的应力张量关系式：
对上式作傅里叶逆变换，并对时间求导，得到：
再根据速度-应力关系、速度-应变关系：
得到解耦分数阶拉普拉斯算子(DFL)粘弹性波动方程为：
其中，
和/>表示x和z方向的速度分量；/>和/>表示应力分量；/>表示介质密度；表示x和z方向的震源项；/>或/>，表示P波或S波；/>和/>表示传播速度和衰减强度；/>和/>表示在参考频率/>处的P波和S波速度；/>和/>表示P波和S波的品质因子；t表示时刻；/>表示拉普拉斯算子；/>表示符号函数；/>表示应变张量。
DFL粘弹性方程中包含了空间独立的拉普拉斯算子，因此它可以准确地描述地震波在衰减介质中的传播。此外，它还能解耦振幅衰减和相位频散项，有利于地震成像中振幅的补偿。此外，分数阶拉普拉斯算子通常采用伪谱法求解，而伪谱法被视为空间域达到无穷阶的有限差分法，可有效抑制数值频散；且在相同模拟精度下，其所需空间采样点数比有限差分法少，所需内存也相对较小。
S2.设置真空自由表面条件和CPML边界条件。
S201.设置真空自由表面界条件的方法包括：将自由表面设置在通过正应力和速度水平分量的采样位置；在自由表面上部添加真空层，并将真空层中的密度设置为零，速度设置为接近于零。
在本实施例中，自由表面是介质和空气接触的不连续分界面，介质内部的应力应变在这个界面消失，其突变性形成了瑞雷波，并使得瑞雷波在浅地表传播。真空自由表面条件的基本思想是在自由表面以上添加真空层，真空层中的拉梅系数和密度设置为零，真空层中的网格点采用与内部网格点相同的格式计算，不需要单独设置自由表面上正应力为零等条件。
在二维交错网格中，将自由表面设置在通过正应力和速度水平分量/>的采样位置，即自由表面位于整数网格线j上，如图2中的加粗黑线。自由表面以上的密度设置为零，而速度接近于零(以避免除以零)以近似真空层，即
式中，表示横向任意的网格剖分点；/>表示自由表面所在的位置；/>表示真空层中网格点距离自由表面的网格点数。
S202.由于计算成本和硬件能力的限制，数值模拟通常在有限域内进行。因此，适当的吸收边界条件是获得准确数值模拟结果的必要条件。本实施例采用卷积完全匹配层(CPML)来消除人工边界的非物理反射。CPML的主要思想是在PML中引入卷积，并以递归的方式计算这个卷积。
定义微分算子和复频移拉伸函数分别为：
式中表示空间方向，取/>控制/>方向的衰减，取表示CPML内的点到CPML内边界的距离；/>表示CPML区域的厚度；/>一般取为2或3；/>表示理论反射系数，取/>控制表面波的吸收，取的正实数；/>控制低频分量的吸收，取一般取为子波的主频；i表示虚数；/>表示圆角频率；vmax表示最大速度；
对作傅里叶逆变换得到：
其中，是/>的傅里叶逆变换，其表达式为：
其中，为单位脉冲函数，/>为Heaviside阶跃函数。
若令
则可表示为：
上式第二项的卷积，在第个时间步可写为如下积分形式：
由于是定义在半个时间步上，故上式可写为如下离散形式：
/>
其中，表示积分变量。
根据递归方法，上式可写为如下递推形式：
因此，CPML应用于所述DFL粘弹性波动方程可表示为：
其中，；，控制j方向的衰减；/>表示某点到卷积完全匹配层内边界的距离；/>表示吸收层厚度；/>表示正整数；/>表示理论反射系数；/>，控制表面波的吸收；/>表示正实数；/>，控制低频分量的吸收；/>表示子波的主频；/>表示空间方向。
S3.以DFL粘弹性方程、真空形式和CPML作为模型，采用SGPS法数值实现瑞雷波正演模拟。
伪谱法能有效地将指数算子转化为乘法算子，因此被广泛应用于求解分数阶拉普拉斯算。在实施例中，分别使用有限差分和交错伪谱法(SGPS)来计算时间导数和分数阶拉普拉斯算子。由图2可知，在交错网格中同一个物理参数的不同分量定义在不同的网格点，因此我们使用算术平均和谐平均方案来计算模型参数。具体计算如下：
其中，表示整网格点处的密度；/>和/>表示x和z方向半网格点的密度；/>和/>为整网格点和半网格点处的算子。
数值实现的过程包括：计算空间导数和分数阶拉普拉斯算子；更新质点的速度和应力，以及吸收层质点的速度和应力，得到瑞雷波数值模拟结果。
空间导数的计算方法为：采用交错伪谱法计算所述空间导数：
其中，和/>分别表示一维傅里叶正变换和逆变换，上标/>表示空间相位左移/右移；/>，上标/>表示当前时刻和前一时刻。
分数阶拉普拉斯算子的计算为：利用二维傅里叶变换计算分数阶拉普拉斯算子：
时间导数采用差分法计算：
其中，和/>分别表示二维傅里叶正变换和逆变换；/>表示时间步长；/>和/>表示空间方向。
更新质点的速度和应力：
更新吸收层质点的速度和应力：
其中，分别表示x和z方向的震源项。循环迭代计算，直到满足迭代次数，得到数值模拟结果。
实施例二
1、对于弹性介质，在本实施例中，首先应用弹性均匀半空间模型验证方案的准确性。弹性均匀半空间模型大小为2000×800个网格点，空间网格间距为1.25m，吸收层为20个网格点，参考频率200Hz处的相速度为，密度为/>，最大模拟时间为1s，时间间隔为0.25ms。震源为垂直方向点力源，震源函数为主频20Hz的Ricker子波，震源和接收器分别位于(0m，0m)和(600m，0m)。对于弹性均匀介质，通过Cagniard-De Hoop技术求解带自由表面条件的格林函数得到解析解。图3为偏移距600m处解析解与数值解的对比结果，虚线表示数值解，实线表示解析解。在图3中，虚线与实线匹配良好，说明所提方案具有较好的准确性。
其次，验证方案抑制瑞雷波的数值频散效应。弹性模型大小为1000×400个网格，空间网格间距为2.5m，时间步长为0.5ms，其他参数与上述模型相同。图4给出t=1s时不同数值方法得到的波场快照，从上到下分别对应网格间距为2.5m和时间步长为0.5ms的交错有限差分(SGFD)法、相同参数设置的交错伪谱(SGPS)法和缩小参数的SGFD法模拟结果；第一和第二列分别表示和/>分量。为了更清晰地观测自由表面的瑞雷波，在图5中展示了自由表面的放大图。在图5a和图5b中，箭头处可以观察到明显的数值伪影，这说明SGFD法在空间步长2.5 m和时间步长0.5 ms时，存在严重的数值频散。当减小SGFD方法的参数（空间和时间步长分别为1.25 m和0.25 ms）时，数值伪影明显减小（如图5e和5f所示）。图5c和5d的波前与图5e和5f的波前效果是基本相同，这表明SGPS方法能够以两倍于SGFD方法的空间和时间步长的设置实现精确的波场外推。图6是偏移距600m处记录的0.2s-0.8s的波形曲线，与波场快照的观察到的情况类似，在空间步长2.5m和时间步长0.5ms的SFFD法得到的波形曲线（虚线）上可以观察到明显的振铃现象；当空间和时间步长减小时，瑞雷波波形畸变明显减小(实线)，与空间步长2.5m和时间步长0.5ms的SGPS法得到(点线)的波形几乎相同。图4和图6表明，该方法有效抑制了真空法与GSLS相结合(SGFD求解时)时产生的数值频散。
2、将本申请应用于粘弹性介质。模型大小及参数与弹性介质第一个模型相同，品质因子取为，震源放置在(1250m，0m)。由于SGFD方法数值求解分数阶拉普拉斯算子比较困难，这里采用SGPS方法进行所数值模拟。图7是弹性和粘弹性均匀半空间模型t=0.5s的波场快照对比，其中RW、P、S和S*分别表示瑞雷波、P波、S波和与自由表面相关的S*波。与弹性波场(图7a和图7b)相比，粘弹性波场(图7c和图7d)的能量衰减严重。
下面验证本申请在均匀粘弹性介质中的准确性。将震源点设置在(0m，0m)，图8是粘弹性均匀半空间模型的炮集记录(图8a)和频散能量图(图8b)，其中黑点为瑞雷波理论相速度。由图8a可知，瑞雷波的能量最强，且随着偏移距增大而逐渐减小。对偏移量为250m-1500m和记录时间为2s的炮集记录(图8a)作相移变换法，得到瑞雷波的频散能量图(图8b)。在图8b中，瑞雷波相速度随频率增加而增加，并在20Hz处相速度为1055m/s；瑞雷波的频散能量与理论相位速度匹配良好，除频率低于10Hz时略有差异，说明本申请能够准确模拟粘弹性介质的瑞雷波。
实施例三
在本实施例中，首先设计一个两层速度递增模型，验证方案处理非均匀粘弹介质的能力。两层模型大小为2000×800个网格点，空间网格间距为1.25m，吸收层为20个网格点，在深度250m处存在一水平界面。震源为垂直方向点力源，位于自由表面的中心，震源函数为主频20Hz的Ricker子波。模拟时长为2s，时间间隔为0.25ms。各层参数如表1所示。图9是t= 0.45s速度场的快照。由图9可知，在自由表面附近，瑞雷波能量强于体波；地震波在两层介质中传播时，由于分界面的存在，在分界面处还形成了反射纵波(RPP)、反射横波(RPS)、透射纵波(TPP)和透射横波(TPS)。数值模拟结果生成了符合勘探地球物理规律的波场，说明该方案适用于非均匀粘弹性介质。
表1
/>
下面通过频散能量曲线验证该方案的准确性。图10为两层模型的炮集记录和频散能量图，其中炮集记录的偏移量为125m- 1000m。在图9a中，可以清晰地观察到瑞雷波的频散和振幅衰减现象，且随着偏移距离的增大而加剧。在图9b可以看到，由于层状介质固有的频散特性，瑞雷波具有多阶模式；在高频范围(f≥20Hz)，基阶模式的相速度随频率的增加而略有增加；频散能量与Knopoff方法计算的解析结果(黑点)相吻合，说明所提方案对于非均匀粘弹介质的模拟结果是准确的。
其次，设计了一组模型来分析衰减对瑞雷波的影响，其中速度保持不变，只考虑Q的变化(参数见图11左列)。假设在两层模型中，/>和/>分别表示第一、二层的品质因子。Model 1作为参考，Model 2用来考察第一层/>对瑞雷波的影响，而Model 3用来考察第二层/>对瑞雷波的影响。图11b是t=0.35s的波场快照。与Model 1的波场快照相比，Model 2的瑞雷波和波场能量随/>的变化发生改变，而未受到/>变化的影响。为具体分析，抽取(875m，0m)处的波形曲线(图12)。如图12所示，瑞雷波的幅值随第一层/>的增大而增大(图12a)，不随第二层/>的变化而变化(图12b)；衰减导致瑞雷波波形发生延迟，这种延迟随第一层/>的增大而减小(图12a)，而基本不随第二层/>的变化而变化(图12b)。产生上述现象是由于瑞雷波的能量沿深度呈指数衰减，并且主要集中在一个波长范围内。
为了考察本申请在大时刻数值模拟中的稳定性，将模拟时间增加到10s。从能量衰减曲线(图13)可以看出，所提方法在10s内仍然保持稳定。
实施例四
为研究瑞雷波在高度异构介质中的传播特性，选取具有真实结构复杂度的Marmousi模型进行模拟测试。模型大小为663×234个网格点，空间网格间距为2m，时间步长为0.1ms，其它参数如图14所示。在 (664 m, 0 m)点处加载一个垂向点力源，震源为主频20hz的Ricker子波，模拟时间为2s。以镜像法(SIM)的模拟结果作为参考。图15为分量的波场快照和炮集记录，其中从上到下分别对应本申请的模拟结果、镜像法（SIM）的模拟结果、本申请与SIM模拟结果的差值。图15表明：本申请的模拟结果与SIM模拟结果差异很小(图15e和15f)。因此，本申请能够有效模拟复杂的非均匀介质。图16为无面波和含有面波的对比。与无面波的波场快照相比(图16a)，含有面波的波场在自由表面附近振幅明显增大(见图16c)；与无面波的炮集记录相比(图16b)，瑞雷波的在炮集记录中呈现扫帚状（图16d箭头），炮集记录残差(图16f)反映了面波对波场的影响。
实施例五
在本实施例中，如图17所示，基于DFL粘弹性方程的瑞雷波数值模拟系统，包括：方程构建模块、边界条件设置模块和数值实现模块；
方程构建模块用于根据地震波在衰减介质中的传播特性，建立DFL粘弹性波动方程。
Xing和Zhu在Kjartansson的常Q频散关系的基础上提出粘性介质声波的新频散关系：
/>
其中是参考频率/>处对应的参考相速度；/>，新频散关系中波数/>的指数均为常数，不再是空间位置的函数。
由于包含P波和S波分量的弹性波也满足上述近似的频散关系，将上述频散关系进行推广，得到粘性介质弹性波的频散关系：
其中表示纵波或横波，取/>。
由此得到波数域的复模量：
将复模量代入弹性波的应力-应变关系：
其中和/>为波数域中的应力张量，/>和/>为波数域中的应变张量，i、j、k代表不同的空间方向。
将上式代入耦合的应力张量关系式：
对上式作傅里叶逆变换，并对时间求导，得到：
再根据速度-应力关系、速度-应变关系：
得到解耦分数阶拉普拉斯算子(DFL)粘弹性波动方程为：
其中，
和/>表示x和z方向的速度分量；/>和/>示应力分量；/>表示介质密度；表示x和z方向的震源项；/>或/>，表示P波或S波；/>和/>表示传播速度和衰减强度；/>和/>表示在参考频率/>处的P波和S波速度；/>和/>表示P波和S波的品质因子；t表示时刻；/>表示拉普拉斯算子；/>表示符号函数；/>表示应变张量。
DFL粘弹性方程中包含了空间独立的拉普拉斯算子，因此它可以准确地描述地震波在衰减介质中的传播。此外，它还能解耦振幅衰减和相位频散项，有利于地震成像中振幅的补偿。此外，分数阶拉普拉斯算子通常采用伪谱法求解，而伪谱法被视为空间域达到无穷阶的有限差分法，可有效抑制数值频散；且在相同模拟精度下，其所需空间采样点数比有限差分法少，所需内存也相对较小。
边界条件设置模块用于设置真空自由表面条件和CPML边界条件。
设置真空自由表面界条件的方法包括：将自由表面设置在通过正应力和速度水平分量的采样位置；在自由表面上部添加真空层，并将真空层中的密度设置为零，速度设置为接近于零。
在本实施例中，自由表面是介质和空气接触的不连续分界面，介质内部的应力应变在这个界面消失，其突变性形成了瑞雷波，并使得瑞雷波在浅地表传播。真空自由表面条件的基本思想是在自由表面以上添加真空层，真空层中的拉梅系数和密度设置为零，真空层中的网格点采用与内部网格点相同的格式计算，不需要单独设置自由表面上正应力为零等条件。
在二维交错网格中，将自由表面设置在通过正应力和速度水平分量/>的采样位置，即自由表面位于整数网格线j上，如图2中的加粗黑线。自由表面以上的密度设置为零，而速度接近于零(以避免除以零)以近似真空层，即
式中，表示横向任意的网格剖分点；/>表示自由表面所在的位置；/>表示真空层中网格点距离自由表面的网格点数。
由于计算成本和硬件能力的限制，数值模拟通常在有限域内进行。因此，适当的吸收边界条件是获得准确数值模拟结果的必要条件。本实施例采用CPML来消除人工边界的非物理反射。CPML的主要思想是在PML中引入卷积，并以递归的方式计算这个卷积。
定义微分算子和复频移拉伸函数（Kuzuoglu和 Mittra）分别为：
式中表示空间方向，取/>控制/>方向的衰减，取表示CPML内的点到CPML内边界的距离；/>表示CPML区域的厚度；/>一般取为2或3；/>表示理论反射系数，取/>控制表面波的吸收，取的正实数；/>控制低频分量的吸收，取一般取为子波的主频；i表示虚数；/>表示圆角频率；vmax表示最大速度；
对作傅里叶逆变换得到：
其中，是/>的傅里叶逆变换，其表达式为：
其中，为单位脉冲函数，/>为Heaviside阶跃函数。
若令
则可表示为：
上式第二项的卷积，在第个时间步可写为如下积分形式：
由于是定义在半个时间步上，故上式可写为如下离散形式：
其中，表示积分变量；
根据递归方法，上式可写为如下递推形式：
因此，CPML应用于所述DFL粘弹性波动方程可表示为：
/>
其中，；，控制j方向的衰减；/>表示某点到卷积完全匹配层内边界的距离；/>表示吸收层厚度；/>表示正整数；/>表示理论反射系数；/>，控制表面波的吸收；/>表示正实数；/>，控制低频分量的吸收；/>表示子波的主频；/>表示空间方向。
数值实现模块以DFL粘弹性方程、真空自由表面边界条件和CPML边界条件作为模型，采用交错伪谱法求解得到瑞雷波数值模拟结果。
伪谱法能有效地将指数算子转化为乘法算子，因此被广泛应用于求解分数阶拉普拉斯算。在实施例中，分别使用有限差分和SGPS法来计算时间导数和分数阶拉普拉斯算子。由图2可知，在交错网格中同一个物理参数的不同分量定义在不同的网格点，因此我们使用算术平均和谐平均方案来计算模型参数。具体计算如下：
其中，表示整网格点处的密度；/>和/>表示x和z方向半网格点的密度；
和/>为整网格点和半网格点处的算子。
数值实现的过程包括：计算空间导数和分数阶拉普拉斯算子；更新质点的速度和应力，以及吸收层质点的速度和应力，得到瑞雷波数值模拟结果。
空间导数的计算方法为：采用交错伪谱法计算所述空间导数：
其中，和/>分别表示一维傅里叶正变换和逆变换，上标/>表示空间相位左移/右移；/>，上标/>表示当前时刻和前一时刻。
分数阶拉普拉斯算子的计算为：利用二维傅里叶变换计算分数阶拉普拉斯算子：
时间导数采用差分法计算：
其中，和/>分别表示二维傅里叶正变换和逆变换；/>表示时间步长；/>和/>表示空间方向。
更新质点的速度和应力：
更新吸收层质点的速度和应力：
其中，分别表示x和z方向的震源项。循环迭代计算，直到满足迭代次数，得到数值模拟结果。
以上所述的实施例仅是对本申请优选方式进行的描述，并非对本申请的范围进行限定，在不脱离本申请设计精神的前提下，本领域普通技术人员对本申请的技术方案做出的各种变形和改进，均应落入本申请权利要求书确定的保护范围内。

Claims (6)

1.基于DFL粘弹性方程的瑞雷波数值模拟方法，其特征在于，包括以下步骤：
根据地震波在衰减介质中的传播特性，建立解耦分数阶拉普拉斯算子DFL粘弹性方程；
根据面波产生条件和消除人工边界反射方法，设置真空自由表面条件和卷积完全匹配层CPML边界条件；
以所述DFL粘弹性方程、所述真空自由表面条件和CPML边界条件作为模型，采用交错伪谱SGPS法进行瑞雷波正演模拟；
设置所述CPML边界条件的方法包括：
定义微分算子和复频移拉伸函数分别为：
式中，/>表示空间方向，取/>；/>控制/>方向的衰减，取；/>表示CPML内的点到CPML内边界的距离，/>；/>表示CPML区域的厚度；/>表示正整数；/>表示理论反射系数，取/>；/>控制表面波的吸收，取，/>表示正实数，/>；/>控制低频分量的吸收，取，/>为子波的主频；i表示虚数；/>表示圆角频率；vmax表示最大速度；
对作傅里叶逆变换得到：
其中，/>是/>的傅里叶逆变换，其表达式为：
其中，/>为单位脉冲函数，/>为Heaviside阶跃函数；
令：
则/>表示为：
其中第二项的卷积，在第/>个时间步写为积分形式：
由于/>是定义在半个时间步上，故写为如下离散形式：
其中，
，/>，；/>表示积分变量；
根据递归方法，将离散形式写为如下递推形式：
即，CPML应用于所述DFL粘弹性方程表示为：
其中，，/>，/>；/>，，/>，控制j方向的衰减；/>表示CPML内的点到CPML内边界的距离；/>表示CPML区域的厚度；/>表示正整数；/>表示理论反射系数；，控制表面波的吸收；/>表示正实数，/>；/>，控制低频分量的吸收；/>表示子波的主频；/>表示空间方向；
所述实现瑞雷波正演模拟的过程包括：
计算空间导数和分数阶拉普拉斯算子；
基于所述空间导数和所述分数阶拉普拉斯算子更新质点的速度和应力，以及吸收层质点的速度和应力，得到所述模拟结果；
所述质点的速度和应力更新公式为：
；
所述吸收层质点的速度和应力更新公式为：
；
其中，分别表示x和z方向的震源项，/>表示时间步长。
2.根据权利要求1所述基于DFL粘弹性方程的瑞雷波数值模拟方法，其特征在于，所述DFL粘弹性方程表示为：
；
其中，
；
和/>表示x和z方向的速度分量；/>和/>表示应力分量；/>表示介质密度；/>表示x和z方向的震源项；/>或/>，表示P波或S波；/>和/>表示传播速度和衰减强度；/>和/>表示在参考频率/>处的P波和S波速度；/>和/>表示P波和S波的品质因子；t表示时刻；/>表示拉普拉斯算子。
3.根据权利要求1所述基于DFL粘弹性方程的瑞雷波数值模拟方法，其特征在于，设置所述真空自由表面条件的方法包括：
将自由表面设置在通过正应力和速度水平分量的采样位置；
在所述自由表面上部添加真空层，并将所述真空层中的密度设置为零。
4.根据权利要求1所述基于DFL粘弹性方程的瑞雷波数值模拟方法，其特征在于，所述空间导数的计算方法为：
采用交错伪谱法计算所述空间导数：
；
其中，和/>分别表示一维傅里叶正变换和逆变换，上标/>表示空间相位左移/右移；
，上标/>表示当前时刻和前一时刻。
5.根据权利要求1所述基于DFL粘弹性方程的瑞雷波数值模拟方法，其特征在于，所述分数阶拉普拉斯算子的计算方法为：
利用二维傅里叶变换计算分数阶拉普拉斯算子：
；
时间导数采用差分法计算：
；
其中，和/>分别表示二维傅里叶正变换和逆变换；/>和/>表示空间方向；/>表示时间步长。
6.基于DFL粘弹性方程的瑞雷波数值模拟系统，所述模拟系统应用权利要求1-5任一项所述的模拟方法，其特征在于，包括：方程构建模块、边界条件设置模块和数值实现模块；
所述方程构建模块用于根据地震波在衰减介质中的传播特性，建立DFL粘弹性方程；
所述边界条件设置模块用于根据面波产生条件和消除人工边界反射方法，设置真空自由表面条件和CPML边界条件；
所述数值实现模块用于以所述DFL粘弹性方程、所述真空自由表面条件和CPML边界条件作为模型，采用交错伪谱SGPS法进行瑞雷波正演模拟。
CN202311146141.XA 2023-09-07 2023-09-07 基于dfl粘弹性方程的瑞雷波数值模拟方法及系统 Active CN116882214B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202311146141.XA CN116882214B (zh) 2023-09-07 2023-09-07 基于dfl粘弹性方程的瑞雷波数值模拟方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202311146141.XA CN116882214B (zh) 2023-09-07 2023-09-07 基于dfl粘弹性方程的瑞雷波数值模拟方法及系统

Publications (2)

Publication Number Publication Date
CN116882214A CN116882214A (zh) 2023-10-13
CN116882214B true CN116882214B (zh) 2023-12-26

Family

ID=88266672

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202311146141.XA Active CN116882214B (zh) 2023-09-07 2023-09-07 基于dfl粘弹性方程的瑞雷波数值模拟方法及系统

Country Status (1)

Country Link
CN (1) CN116882214B (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008059991A (ja) * 2006-09-01 2008-03-13 Canon Inc プラズマ処理装置及びプラズマ処理方法
CN102421938A (zh) * 2009-05-15 2012-04-18 株式会社岛津制作所 表面波等离子体cvd设备以及成膜方法
CN102749643A (zh) * 2011-04-22 2012-10-24 中国石油天然气股份有限公司 一种波动方程正演的瑞利面波频散响应计算方法及其装置
CN106499917A (zh) * 2016-11-01 2017-03-15 中国人民解放军国防科学技术大学 一种非限位的多功能粘弹性支撑结构群

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8428922B2 (en) * 2010-02-05 2013-04-23 Seiko Epson Corporation Finite difference level set projection method on multi-staged quadrilateral grids
JP5670832B2 (ja) * 2011-05-24 2015-02-18 富士通株式会社 シミュレーション方法、シミュレーション装置及びシミュレーションプログラム
US20180210101A1 (en) * 2016-12-30 2018-07-26 China Petroleum & Chemical Corporation System and method for seismic inversion

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008059991A (ja) * 2006-09-01 2008-03-13 Canon Inc プラズマ処理装置及びプラズマ処理方法
CN102421938A (zh) * 2009-05-15 2012-04-18 株式会社岛津制作所 表面波等离子体cvd设备以及成膜方法
CN102749643A (zh) * 2011-04-22 2012-10-24 中国石油天然气股份有限公司 一种波动方程正演的瑞利面波频散响应计算方法及其装置
CN106499917A (zh) * 2016-11-01 2017-03-15 中国人民解放军国防科学技术大学 一种非限位的多功能粘弹性支撑结构群

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于分数阶拉普拉斯算子的黏弹性各向异性介质正演方法研究;温磊;中国优秀硕士学位论文全文数据库 基础科学辑(第02期);第A011-635页 *
基于粘滞波动方程的衰减补偿逆时偏移方法研究;赵学彬;中国优秀硕士学位论文全文数据库 基础科学辑(第02期);第A011-663页 *

Also Published As

Publication number Publication date
CN116882214A (zh) 2023-10-13

Similar Documents

Publication Publication Date Title
Chen et al. Elastic least-squares reverse time migration via linearized elastic full-waveform inversion with pseudo-Hessian preconditioning
CA2726453C (en) Removal of surface-wave noise in seismic data
AU2007302695B2 (en) Iterative inversion of data from simultaneous geophysical sources
Yang et al. Viscoacoustic reverse time migration using a time-domain complex-valued wave equation
Borisov et al. Application of 2D full-waveform inversion on exploration land data
KR20130060231A (ko) 지구 물리학적 데이터의 반복 반전의 아티팩트 감소
Liu et al. An efficient step-length formula for correlative least-squares reverse time migration
CA2972033C (en) Multistage full wavefield inversion process that generates a multiple free data set
Zhao et al. A stable approach for Q-compensated viscoelastic reverse time migration using excitation amplitude imaging condition
CN105652321A (zh) 一种粘声各向异性最小二乘逆时偏移成像方法
Yao et al. An effective absorbing layer for the boundary condition in acoustic seismic wave simulation
Sun et al. Preconditioning least-squares RTM in viscoacoustic media by Q-compensated RTM
Fathalian et al. Q-compensated reverse time migration in tilted transversely isotropic media
Métivier et al. Combining asymptotic linearized inversion and full waveform inversion
Djebbi et al. Frequency domain multiparameter acoustic inversion for transversely isotropic media with a vertical axis of symmetry
Mu et al. Attenuation compensation and anisotropy correction in reverse time migration for attenuating tilted transversely isotropic media
CN109738950A (zh) 基于三维稀疏聚焦域反演的噪声型数据一次波反演方法
Bai et al. Gaussian beam reconstruction of seismic data
CN116882214B (zh) 基于dfl粘弹性方程的瑞雷波数值模拟方法及系统
Greer et al. Improving migration resolution by approximating the least-squares Hessian using nonstationary amplitude and frequency matching
Dantas et al. A multiscale approach to full-waveform inversion using a sequence of time-domain misfit functions
Slob et al. Unified elimination of 1D acoustic multiple reflection
Du et al. Seismic attenuation compensation with spectral-shaping regularization
Lecomte Hybrid modeling with ray tracing and finite difference
Margrave et al. Full waveform inversion with wave equation migration and well control

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