CN115081267B - 基于优化的时空变步长有限差分地震波数值模拟方法 - Google Patents

基于优化的时空变步长有限差分地震波数值模拟方法 Download PDF

Info

Publication number
CN115081267B
CN115081267B CN202210551413.3A CN202210551413A CN115081267B CN 115081267 B CN115081267 B CN 115081267B CN 202210551413 A CN202210551413 A CN 202210551413A CN 115081267 B CN115081267 B CN 115081267B
Authority
CN
China
Prior art keywords
stress
seismic
time
wave
representing
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
CN202210551413.3A
Other languages
English (en)
Other versions
CN115081267A (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.)
Inner Mongolia Agricultural University
Original Assignee
Inner Mongolia Agricultural 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 Inner Mongolia Agricultural University filed Critical Inner Mongolia Agricultural University
Priority to CN202210551413.3A priority Critical patent/CN115081267B/zh
Publication of CN115081267A publication Critical patent/CN115081267A/zh
Application granted granted Critical
Publication of CN115081267B publication Critical patent/CN115081267B/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
    • 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
    • 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
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05—Geographic models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
    • 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/10—Noise analysis or noise optimisation
    • 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)
  • Software Systems (AREA)
  • Geometry (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • General Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Computer Graphics (AREA)
  • Evolutionary Computation (AREA)
  • Operations Research (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Computer Hardware Design (AREA)
  • Remote Sensing (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开基于优化的时空变步长有限差分地震波数值模拟方法，包括以下步骤：S1根据研究对象建立模型，并依据介质分布情况，将整个区域划分成若干个子区域，所述子区域包括低速区域和非低速区域；其中，低速区域采用细网格和短时间步长，非低速区域采用粗网格和长时间步长，低速区域与非低速区域间采用的时间步长变化比和空间步长变化比为相同值；S2根据空间步长变化比和空间差分格式的精度，在粗网格和细网格交界处设置两个过渡区，靠近细网格的过渡区为第一过渡区，靠近粗网格的过渡区为第二过渡区；通过优化的时‑空变步长有限差分法消除数值频散和高频非物理震荡波，从而提高地震波数值模拟的精度和稳定性。

Description

基于优化的时空变步长有限差分地震波数值模拟方法
技术领域
本发明属于地震波数值模拟领域，尤其涉及一种基于优化的时空变步长有限差分地震波数值模拟方法。
背景技术
传统的有限差分法采用统一均匀的空间网格步长来离散整个模拟区域，用于时间迭代的时间步长也是统一的。这时空间网格的大小将由低速区决定，用此步长来离散整个计算区域将会在模型高速区导致空间过采样，从而浪费计算内存和效率。根据稳定性条件的关系，时间步长与空间网格步长成正比，因此较细的空间网格步长也需要较短的时间步长。如果将较短的时间步长用于整个计算区域，同样造成计算时间上的浪费。
为了降低地震波数值模拟所需的计算内存和提高计算效率，可以采用时空双变网格技术，即在模型低速区采用较细的网格和较短的时间步长，在模型高速区采用较粗的网格和较长的时间步长。时空双变网格技术的关键在于如何处理时间和空间采样率突变边界处的波场交换问题。
目前，用于实现空间变网格和时间变步长边界处波场交换的方法有插值法和变系数法。插值法容易产生不稳定和数值噪音。变系数法不进行波场插值，而通过在刚进入变网格区域的时候采用重新计算差分系数，得到精细网格所需要的波场值。这种算法比插值法有更高的稳定性。
现有的空间变网格有限差分法和时间变步长有限差分法的差分系数都是基于Taylor级数展开或Langrange插值方法求取的。这些方法不仅容易产生数值频散，而且在长时间计算时容易产生不稳定现象。
发明内容
针对现有方法的不足，本发明提出一种基于优化的时空变步长有限差分地震波数值模拟方法，旨在解决上述存在的问题。
本发明是这样实现的，基于优化的时-空变步长有限差分地震波数值模拟方法，包括以下步骤：
S1根据模型参数的分布情况，将计算区划分成若干个子区域，在低速区域采用细网格和短时间步长，其它区域采用粗网格和长时间步长，其中时间步长变化比和空间步长变化比取相同值；
S2根据空间步长变化比和空间差分格式的精度，在粗、细网格交界处附近设置两个过渡区；
S3采用空间定步长有限差分格式计算第一过渡区外粗、细网格点上的地震波速度和应力的空间导数；
S4当波场传播到第一过渡区内时，采用优化的空间变步长有限差分格式计算过渡区1内细网格点上的地震波速度和应力的空间导数，这时粗网格点上缺失点处的波场值可以利用值插值法求取；
S5利用选择性滤波法消除计算过程中的非物理高频震荡波，用来保证计算的稳定性；
S6采用优化的时间变步长有限差分格式迭代计算第二过渡区内粗网格点上的地震波速度和应力在下一个短时间步长处的值，并存贮在内存中，用于计算第一过渡区内细网格点的地震波速度和应力在下一个短时间步长处的值；
S7采用时间定步长有限差分格式迭代计算粗、细网格点上的地震波速度和应力在下一时刻的值；
S8获得基于优化的时-空变步长有限差分地震波数值模拟结果。
在步骤S1中，空间可变网格和时间变步长的实现可以采用如下方法：1.利用精细网格剖分介质；2.给定某个阈值对剖分的介质进行扫描，对达到要求的区域保留地质参数，而其它区域进行粗网格的重采样，采样率可取任意整数；3.根据网剖分布情况，在粗网格区域采用长时间步长，细网格采用短时间步长，其中时间步长变化比和空间步长变化比取相同值。
在步骤S3中，采用如下2N+1点中心差分格式计算第一过渡区外粗、细网格点上地震波速度和应力的空间导数：
式中，U＝[vx,vz,σxx,σzz,σxz]T，vx表示地震波速度的水平分量，vz表示地震波速度的垂直平分量，σxx和σzz表示正应力，σxz表示剪切应力，Δx为空间步长，aj为中心差分系数，并且有aj＝-a-j。
在步骤S4中，采用如下的2N+1点差分格式计算过渡区1内细网格点上地震波速和应力的空间导数：
式中，U＝[vx,vz,σxx,σzz,σxz]T，vx表示地震波速度的水平分量，vz表示地震波速度的垂直平分量，σxx和σxx表示正应力，σxz表示剪切应力，Δx为细网格的步长，bj为差分系数，并且有aj＝-a-j，l取决于网格点的位置和粗、细网格步长比；
为了方便讨论bj的计算步骤，考虑下面的一维一阶波动方程：
式中，U表示波场变量，t表示时间，x表示坐标，A表示系数；
对方程(5)中的空间导数，有如下2N+1点优化空间的变步长格式为：
式中，U表示波场变量，Δx为细网格的步长，bj为差分系数，l取决于网格点的位置和粗、细网格步长比；
根据Fourier分析，将(6)式变换到波数域并整理后可得：
式中k*和k分别为有效空间波数和准确空间波数；
为了让k*Δx和kΔx的误差最小，要求差分系数的选择保证下式误差积分最小
对(8)式求最小极值，即
求解(9)式确定的线性方程组，便可得到差分系数bj。
在步骤S5中，采用下式对地震波速度和应力进行选择性滤波处理：
Ud(x,z)＝U(x,z)-σDu(x,z) (10)
式中，U＝[vx,vz,σxx,σzz,σxz]T，vx表示地震波速度的水平分量，vz表示地震波速度的垂直平分量，σxx和σxx表示正应力，σxz表示剪切应力，Ud是滤波后的地震波速度和应力，σ是滤波衰减因子，cj是滤波算法的系数满足cj＝c-j。
在步骤S6中，采用下式计算第二过渡区内粗网格点上地震波速度和应力在第一个短时间步长处的值：
式中，U＝[vx,vz,σxx,σzz,σxz]T，vx表示地震波速度的水平分量，vz表示地震波速度的垂直平分量，σxx和σxx表示正应力，σxz表示剪切应力，q表示长、短时间步长变化比，Δt表示短时间步长，Un+1/q表示过度带内粗网格点上(n+1/q)Δt时刻的地震波速度和应力，αj为差分系数，p取决于时间差分精度，由下面的波动方程给出：
式中，U＝[vx,vz,σxx,σzz,σxz]T
其中，ρ是密度，λ和μ为介质的拉梅系数；
为了方便讨论αj的计算步骤，考虑一维一阶波动方程：
式中，U表示波场变量，t表示时间，x表示坐标，A表示系数；
对于方程(14)，采用下式计算第二过渡区内粗网格点上波场变量在第一个短时间步长处的值：
式中，U表示波场变量，q表示长，短时间步长变化比，Δt表示短时间步长，αj为差分系数，p取决于时间差分精度。
根据Fourier分析，将(15)式变换到波数域并整理后可得：
为了确定αj，使下列积分取最小值：
式中σ是权系数。对积分求最小极值得
(dE/dαj)＝0 (18)
求解(18)式确定的线性方程组，便可得到时间差分系数αj；
同理可以计算第二过渡区内粗网格点上的地震波速度和应力在其它短时间步长处的值。
在步骤S7中，采用如下差分格式计算粗、细网格点上的地震波速度和应力在下一时刻的值：
式中，Un+1表示粗、细网格点上(n+1)Δt时刻的地震波速度和应力，Un表示粗、细网格点上nΔt时刻的地震波速度和应力，表示粗网格点上(n-j)Δt时刻地震波速度和应力对时间的导数，Δt表示短时间步长，βj表示差分系数，p取决于时间差分精度，/>由方程(13)给出。
完成上述步骤，获得基于优化的时-空变步长有限差分地震波数值模拟结果。
与现有技术相比，本发明的有益效果是：通过优化的时-空变步长有限差分法消除数值频散和高频非物理震荡波，从而提高地震波数值模拟的精度和稳定性。
附图说明
图1为实施例中空间变步长区域处理示意图；
图2为实施例中时间变步长区域处理示意图；
图3为试验例1中建立的二维速度模型；
图4为试验例1中空间变步长网格剖分示意图；
图5为试验例1中模拟得到的水平分量波场快照；
图6为试验例1中模拟得到的垂直分量波场快照；
图7为试验例1中模拟得到的水平分量地震记录；
图8为试验例1中模拟得到的垂直分量地震记录；
图9为试验例2中建立的二维速度模型；
图10为试验例2中空间变步长网格剖分示意图；
图11为试验例2中模拟得到的水平分量波场快照；
图12为试验例2中模拟得到的垂直分量波场快照。
具体实施方式
为了进一步了解本发明的内容，下面以粗、细网格比r＝2为例，对本发明进行进一步详细说明。应当理解，此处所描述的实施例仅仅用以解释本发明，并不用于限定本发明。
实施例
本发明提供一种优化的时-空变步长有限差分地震波数值模拟方法，包括以下步骤：
S1根据模型参数的分布情况，将计算区划分成若干个子区域，低速区域采用细网格和短时间步长，其它区域采用粗网格和长时间步长，其中时间步长变化比和空间步长变化比取相同值；
S2根据空间步长变化比和空间差分格式的精度，在粗、细网格交界处设置如图1所示的两个过渡区；
S3采用空间定步长有限差分格式计算第一过渡区外粗、细网格点上的地震波速度和应力的空间导数；
S4当地震波传播至第一过渡区内时，采用优化的空间变步长有限差分格式计算第一过渡区内细网格点上的地震波速度和应力的空间导数，这时粗网格中缺失点处的波场值采用插值法求取；
S5利用选择性滤波法消除计算过程中的非物理高频震荡波，用来保证计算的稳定性；
S6采用优化的时间变步长有限差分格式迭代计算图1中第二过渡区内粗网格点上的地震波速度和应力在下一个短时间步长处的值，并存贮在内存中，用于计算第一过渡区内细网格点的地震波速度和应力在下一个短时间步长处的值。
S7采用时间定步长有限差分格式迭代计算粗、细网格点上的地震波速度和应力在下一时刻的值。
在步骤S1中，空间可变网格和时间变步长的实现可以采用如下方法：1.利用精细网格剖分介质；2.给定某个阈值对剖分的介质进行扫描，对达到要求的区域保留地质参数，而其它区域进行粗网格的重采样，采样率为2；3.根据网剖分布情况，在粗网格区域采用长的时间步长，细网格采用短时间步长，时间步长变化比取2。
在步骤S3中，在粗、细网格边界附近设置两个过渡区(图1中的灰色区域)，然后利用如下的7点空间定步长有限差分格式计算第一过渡区外粗、细网格点上地震速度和应力的空间导数；
式中，U＝[vx,vz,σxx,σzz,σxz]T，vx表示地震波速度的水平分量，vz表示地震波速度的垂直平分量，σxx和σxx表示正应力，σxz表示剪切应力，x表示水平坐标，Δx为空间步长，aj表示差分系数，具体的值见表1。
表1
在步骤S4中，采用如下差分格式计算图1中第一过渡区内A点上的地震波速度和应力的空间导数：
其中，U＝[vx,vz,σxx,σzz,σxz]T，vx表示地震波速度的水平分量，vz表示地震波速度的垂直平分量，σxx和σxx表示正应力，σxz表示剪切应力，x表示水平坐标，Δx为空间步长，为差分系数，并且有bj＝-bj，U(x+jΔx)由图1中○上的波场变量给出；
为了方便讨论的计算步骤，考虑下面的一维一阶波动方程：
对方程(5)中的空间导数，有如下2N+1点优化空间的变步长格式为：
式中，U表示波场变量，t表示时间，x表示坐标，A表示系数，U(x+jΔx)由图1中○上的波场变量给出，为差分系数，并且有bj＝-bj；
根据Fourier分析，将(6)式变换到波数域并整理后可得
为了让k*Δx和kΔx的误差最小，要求差分系数的选择保证下式误差积分最小
对积分方程(8)取极小值，即
这时求解方程(9)对应的线性方程组可得优化的差分系数具体的值见表1；同理可以得到B点上的空间差分系数/>具体的值见表1；
过渡带内细网格点B′和C′上的空间导数同样可以采用上面的方法求取，这时粗网格区域中缺失点处的波场值(如图1中C、D、E点上的值)采用Tam and Kurbatskii提出的优化的对称点插值法求取。
在步骤S5中，有如下7点选择性中心滤波算法：
式中，U＝[vx,vz,σxx,σzz,σxz]T，vx表示地震波速度的水平分量，vz表示地震波速度的垂直平分量，σxx和σxx表示正应力，σxz表示剪切应力，Ud是滤波后的地震波速度和应力，σ∈[0,1]是滤波衰减因子，cj为滤波算法的权系数，其中σ＝0.2π时过渡带外网格点上波场变量的滤波衰减系数cj和过渡带内波场变量的滤波衰减系数的值见表2；
表2
在步骤S6中，如图2所示假设已知n时刻的地震波速度和应力，利用第二过渡区内(n-3)、(n-2)、(n-1)、n时刻的地震波速度和应力，采用如下优化的变步长有限差分格式计算第二过渡区内粗网格点上半个时间步长处的地震波速度和应力(图2中×上的值)：
式中，U＝[vx,vz,σxx,σzz,σxz]T，vx表示地震波速度的水平分量，vz表示地震波速度的垂直平分量，σxx和σxx表示正应力，σxz表示剪切应力，Δt表示短时间步长，Un+1/2表示粗网格点上(n+1/2)Δt时刻的地震波速度和应力，αj为差分系数，由下面的波动方程给出：
式中，U＝[vx,vz,σxx,σzz,σxz]T
其中，ρ是密度，λ和μ为介质的拉梅系数；
为了方便讨论αj的计算步骤，考虑一维一阶波动方程：
式中，U表示波场变量，t表示时间，x表示坐标，A表示系数；
对于方程(14)，采用下式计算第二过渡区内粗网格点上波场变量在第一个短时间步长处的值：
将上式变换到波数域并整理后可得：
为了确定αj，使下列积分取最小值：
式中σ是权系数。对积分求最小极值得
(dE/dαj)＝0 (18)
求解上式确定的线性方程组，得到变时间步长差分系数αj，σ＝0.42时αj的值见表3。
表3
在步骤S7中，利用(n-3)、(n-2)、(n-1)、n时刻的地震波速度和应力，采用如下的时间定步长有限差分格式计算(n+1)时刻的波场值：
式中，Δt表示短时间步长，表示(n+1)Δt时刻的地震波速度和应力，βj为差分系数，具体值见表3。
试验例1：地表存在局部低速带时的地震波数值模拟
以理论模型来进行验证本发明的有效性，建立一个如图3所示，地表存在局部低速带的速度模型，模型大小为5000m×3000m，第一层速度为Vp＝2200m/s、Vs＝1300m/s，第二层速度为Vp＝2500m/s、Vs＝1400m/s，第三层速度为Vp＝3000m/s、Vs＝1700m/s，低速带速度为Vp＝1000m/s、Vs＝800m/s，震源位置为(1000m,10m)，接收点位于地下5m处。
首先对模型进行如图4所示的变步长网格剖分，细网格间距取5m，粗网格间距取10m。利用优化的时-空变步长有限差分法模拟地震波场。
图5为1.6s时刻的水平分量波场快照对比图，其中，图5(a)为优化的时-空变步长有限差分法模拟得到的结果、图5(b)为定步长有限差分法模拟得到的结果。
图6为1.6s时刻的垂直分量波场快照对比图，其中，图6(a)为优化的时-空变步长有限差分法模拟得到的结果、图6(b)为定步长有限差分法模拟得到的结果。
图7为水平分量地震记录对比图，其中，图7(a)为优化的时-空变步长有限差分法模拟得到的结果、图7(b)为定步长有限差分法模拟得到的结果。
图8为垂直分量地震记录对比图，其中，图8(a)为优化的时-空变步长有限差分法模拟得到的结果、图8(b)为定步长有限差分法模拟得到的结果。
对比优化的时-空变步长有限差分地震波数值模拟结果和定步长有限差分地震波数值模拟结果可知优化的时-空变步长有限差分法能够有效压制数值频散和高频非物理震荡波，从而提高计算精度。
试验例2：地下存在一个含水溶洞时的地震波数值模拟
以理论模型来进行验证本发明的有效性，建立一个如图9所示地下存在一个含水溶洞的速度模型，模型大小为6000m×4000m，第一层速度为Vp＝2200m/s、Vs＝1400m/s，第二层速度为Vp＝2500m/s、Vs＝1600m/s，第三层速度为Vp＝3000m/s、Vs＝1800m/s，含水溶洞内地震波速度为Vp＝1500m/s，震源位置为(3000m,10m)，接收点位于地下5m处。
首先对模型进行如图10所示的变步长网格剖分，细网格间距取5m，粗网格间距取10m。利用优化的时-空变步长有限差分法模拟地震波场。
图11为1.8s时刻的水平分量波场快照对比图，其中，图11(a)为优化的时-空变步长有限差分法模拟得到的结果、图11(b)为定步长有限差分法模拟得到的结果。
图12为1.8s时刻的垂直分量波场快照对比图，其中，图12(a)为优化的时-空变步长有限差分法模拟得到的结果、图12(b)为定步长有限差分法模拟得到的结果。
对比优化的时-空变步长有限差分地震波数值模拟结果和定步长有限差分地震波数值模拟结果可知优化的时-空变步长有限差分法能够有效压制数值频散和高频非物理震荡波，并且能保证长时间计算的稳定性。

Claims (7)

1.基于优化的时空变步长有限差分地震波数值模拟方法，其特征在于，包括以下步骤：
S1根据研究对象建立模型，并依据介质分布情况，将整个区域划分成若干个子区域，所述子区域包括低速区域和非低速区域；其中，低速区域采用细网格和短时间步长，非低速区域采用粗网格和长时间步长，低速区域与非低速区域间采用的时间步长变化比和空间步长变化比为相同值；
S2根据空间步长变化比和空间差分格式的精度，在粗网格和细网格交界处设置两个过渡区，靠近细网格的过渡区为第一过渡区，靠近粗网格的过渡区为第二过渡区；
S3采用空间定步长有限差分格式，计算第一过渡区外粗网格点和细网格点上的地震波速度和应力的空间导数；
S4当波场传播到第一过渡区内时，采用优化的空间变步长有限差分格式计算第一过渡区内细网格点上的地震波速度和应力的空间导数，利用插值法计算粗网格点上缺失点处的波场值；
S5利用选择性滤波法消除计算过程中的非物理高频震荡波；
S6采用优化的时间变步长有限差分格式迭代计算第二过渡区内粗网格点上的地震波速度和应力在下一个短时间步长处的值，并存贮在内存中，用于计算第一过渡区内细网格点的地震波速度和应力在下一个短时间步长处的值；
S7采用时间定步长有限差分格式迭代计算粗网格点和细网格点上的地震波速度和应力在下一时刻的值；
S8获得基于优化的时空变步长有限差分地震波数值模拟结果。
2.根据权利要求1所述的基于优化的时空变步长有限差分地震波数值模拟方法，其特征在于，在步骤S1中，具体包括以下步骤：
S11利用精细网格剖分介质；
S12给定阈值对剖分的介质进行扫描，对达到要求的区域保留地质参数，剩余区域进行粗网格的重采样，采样率为任意整数；
S13根据网剖分布情况，在粗网格区域采用长的时间步长，细网格区域采用短时间步长，时间步长变化比和空间步长变化比为相同值。
3.根据权利要求1所述的基于优化的时空变步长有限差分地震波数值模拟方法，其特征在于，在步骤S3中，计算第一过渡区外粗网格点和细网格点上的地震波速度和应力的空间导数，采用如下2N+1点中心差分格式：
式中，U＝[vx,vz,σxx,σzz,σxz]T，vx表示地震波速度的水平分量，vz表示地震波速度的垂直平分量，σxx和σzz表示正应力，σxz表示剪切应力，Δx为空间步长，aj为中心差分系数，并且有aj＝-a-j。
4.根据权利要求1所述的基于优化的时空变步长有限差分地震波数值模拟方法，其特征在于，在步骤S4中，计算第一过渡区内细网格点上的地震波速度和应力的空间导数，采用如下2N+1点差分格式：
式中，U＝[vx,vz,σxx,σzz,σxz]T，vx表示地震波速度的水平分量，vz表示地震波速度的垂直平分量，σxx和σxx表示正应力，σxz表示剪切应力，Δx为细网格的步长，bj为差分系数，并且有aj＝-a-j，l取决于网格点的位置和粗、细网格步长比；
其中，bj的计算步骤应用一维一阶波动方程：
式中，U表示波场变量，t表示时间，x表示坐标，A表示系数；
对方程(5)中的空间导数，有如下2N+1点优化空间的变步长格式为：
式中，U表示波场变量，Δx为细网格的步长，bj为差分系数，l取决于网格点的位置和粗、细网格步长比；
根据Fourier分析，将(6)式变换到波数域并整理后可得：
式中，k*和k分别为有效空间波数和准确空间波数；
为了让k*Δx和kΔx的误差最小，要求差分系数的选择保证下式误差积分最小：
对(8)式求最小极值，即
求解(9)式确定的线性方程组，得到差分系数bj。
5.根据权利要求1所述的基于优化的时空变步长有限差分地震波数值模拟方法，其特征在于，在步骤S5中，采用下式对地震波速度和应力进行选择性滤波处理：
Ud(x,z)＝U(x,z)-σDu(x,z) (10)
式中，U＝[vx,vz,σxx,σzz,σxz]T，vx表示地震波速度的水平分量，vz表示地震波速度的垂直平分量，σxx和σxx表示正应力，σxz表示剪切应力，Ud是滤波后的地震波速度和应力，σ是滤波衰减因子，cj是滤波算法的系数满足cj＝c-j。
6.根据权利要求1所述的基于优化的时空变步长有限差分地震波数值模拟方法，其特征在于，在步骤S6中，采用下式计算第二过渡区内粗网格点上地震波速度和应力在第一个短时间步长处的值：
式中，U＝[vx,vz,σxx,σzz,σxz]T，vx表示地震波速度的水平分量，vz表示地震波速度的垂直平分量，σxx和σxx表示正应力，σxz表示剪切应力，q表示长、短时间步长变化比，Δt表示短时间步长，Un+1/q表示过渡带内粗网格点上(n+1/q)Δt时刻的地震波速度和应力，αj为差分系数，p取决于时间差分精度，由下面的波动方程给出：
式中，U＝[vx,vz,σxx,σzz,σxz]T，
其中，ρ是密度，λ和μ为介质的拉梅系数；
αj的计算步骤应用一维一阶波动方程：
式中，U表示波场变量，t表示时间，x表示坐标，A表示系数；
对于方程(14)，采用下式计算第二过渡区内粗网格点上波场变量在第一个短时间步长处的值：
式中，U表示波场变量，q表示长短时间步长变化比，Δt表示短时间步长，αj为差分系数，p取决于时间差分精度；
根据Fourier分析，将(15)式变换到波数域并整理后可得：
为了确定αj，使下列积分取最小值：
式中σ是权系数；
对积分求最小极值得
(dE/dαj)＝0 (18)
求解(18)式确定的线性方程组，便可得到时间差分系数αj；
同理计算第二过渡区内粗网格点上的地震波速度和应力在其它短时间步长处的值。
7.根据权利要求1所述的基于优化的时空变步长有限差分地震波数值模拟方法，其特征在于，在步骤S7中，采用如下差分格式计算粗网格点和细网格点上的地震波速度和应力在下一时刻的值：
式中，Un+1表示粗网格点和细网格点上(n+1)Δt时刻的地震波速度和应力，Un表示粗网格点和细网格点上nΔt时刻的地震波速度和应力，表示粗网格点上(n-j)Δt时刻地震波速度和应力对时间的导数，Δt表示短时间步长，βj表示差分系数，p取决于时间差分精度，/>由方程(13)给出。
CN202210551413.3A 2022-05-18 2022-05-18 基于优化的时空变步长有限差分地震波数值模拟方法 Active CN115081267B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210551413.3A CN115081267B (zh) 2022-05-18 2022-05-18 基于优化的时空变步长有限差分地震波数值模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210551413.3A CN115081267B (zh) 2022-05-18 2022-05-18 基于优化的时空变步长有限差分地震波数值模拟方法

Publications (2)

Publication Number Publication Date
CN115081267A CN115081267A (zh) 2022-09-20
CN115081267B true CN115081267B (zh) 2023-08-29

Family

ID=83249776

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210551413.3A Active CN115081267B (zh) 2022-05-18 2022-05-18 基于优化的时空变步长有限差分地震波数值模拟方法

Country Status (1)

Country Link
CN (1) CN115081267B (zh)

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102183790A (zh) * 2011-02-12 2011-09-14 中国石油大学(华东) 基于时空双变网格的弹性波正演模拟技术
CN104977607A (zh) * 2014-04-09 2015-10-14 中国石油集团东方地球物理勘探有限责任公司 利用变步长网格声波波场模拟的时间域全波形反演方法
CN105277980A (zh) * 2014-06-26 2016-01-27 中石化石油工程地球物理有限公司胜利分公司 高精度空间和时间任意倍数可变网格有限差分正演方法
CN107807392A (zh) * 2017-09-28 2018-03-16 中国海洋石油总公司 一种自适应抗频散的分块时空双变逆时偏移方法
CN108108331A (zh) * 2017-12-13 2018-06-01 国家深海基地管理中心 一种基于拟空间域弹性波方程的有限差分计算方法
CN110109177A (zh) * 2019-06-05 2019-08-09 吉林大学 基于旋转时空双变网格有限差分法的地震波正演模拟方法
WO2019233475A1 (zh) * 2018-06-08 2019-12-12 中国科学院地质与地球物理研究所 拓展有限差分稳定性条件的波场模拟方法、设备及介质
CN112329311A (zh) * 2020-11-10 2021-02-05 中国石油大学(华东) 地震波传播有限差分模拟方法、装置及计算机存储介质
CN114089419A (zh) * 2020-08-24 2022-02-25 中国石油化工集团有限公司 优化的变网格地震正演模拟方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11281825B2 (en) * 2020-06-30 2022-03-22 China Petroleum & Chemical Corporation Computer-implemented method for high speed multi-source loading and retrieval of wavefields employing finite difference models

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102183790A (zh) * 2011-02-12 2011-09-14 中国石油大学(华东) 基于时空双变网格的弹性波正演模拟技术
CN104977607A (zh) * 2014-04-09 2015-10-14 中国石油集团东方地球物理勘探有限责任公司 利用变步长网格声波波场模拟的时间域全波形反演方法
CN105277980A (zh) * 2014-06-26 2016-01-27 中石化石油工程地球物理有限公司胜利分公司 高精度空间和时间任意倍数可变网格有限差分正演方法
CN107807392A (zh) * 2017-09-28 2018-03-16 中国海洋石油总公司 一种自适应抗频散的分块时空双变逆时偏移方法
CN108108331A (zh) * 2017-12-13 2018-06-01 国家深海基地管理中心 一种基于拟空间域弹性波方程的有限差分计算方法
WO2019233475A1 (zh) * 2018-06-08 2019-12-12 中国科学院地质与地球物理研究所 拓展有限差分稳定性条件的波场模拟方法、设备及介质
CN110109177A (zh) * 2019-06-05 2019-08-09 吉林大学 基于旋转时空双变网格有限差分法的地震波正演模拟方法
CN114089419A (zh) * 2020-08-24 2022-02-25 中国石油化工集团有限公司 优化的变网格地震正演模拟方法
CN112329311A (zh) * 2020-11-10 2021-02-05 中国石油大学(华东) 地震波传播有限差分模拟方法、装置及计算机存储介质

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于自适应网格的仿真型有限差分地震波数值模拟;刘志强等;地球物理学报;第59卷(第12期);4654-4665 *

Also Published As

Publication number Publication date
CN115081267A (zh) 2022-09-20

Similar Documents

Publication Publication Date Title
US7675818B2 (en) Method for tomographic inversion by matrix transformation
Amorocho Nonlinear hydrologic analysis
Nyman et al. The display-equalized filter for frequency-time analysis
CN110109177B (zh) 基于旋转时空双变网格有限差分法的地震波正演模拟方法
Beux et al. Exact-gradient shape optimization of a 2-D Euler flow
CN109459789A (zh) 基于振幅衰减与线性插值的时间域全波形反演方法
CN115081267B (zh) 基于优化的时空变步长有限差分地震波数值模拟方法
CN104459791A (zh) 一种基于波动方程的小尺度大模型正演模拟方法
CN110244353B (zh) 一种基于稀疏范数优化算法的地震数据规则化方法
CN108445539B (zh) 一种消除地震子波旁瓣干扰的方法、设备及系统
CN103675905A (zh) 一种优化系数获取方法、装置及相关波场模拟方法、装置
WO2020057286A1 (zh) 波场正演模拟方法及装置
CN109164488B (zh) 一种梯形网格有限差分地震波场模拟方法
RAI et al. Conservative high-order-accurate finite-difference methods for curvilinear grids
CN113341463B (zh) 一种叠前地震资料非平稳盲反褶积方法及相关组件
CN113221392B (zh) 一种非均匀粘声波在无限域内传播模型的构建方法
Rauf et al. Calculation of Lyapunov exponents through nonlinear adaptive filters
Wu et al. Trapezoid-grid finite difference frequency domain method for seismic wave simulation
Darby et al. Application of dynamic programming to the problem of plane wave propagation in a layered medium
CN113671574B (zh) 基于频域的最小二乘共轭梯度迭代的地震正演模拟方法
CN104977608B (zh) 利用固定网格声波波场模拟的时间域全波形反演方法
CN113917526B (zh) 基于非分裂完全匹配层吸收边界的正演模拟方法
CN115390135A (zh) 一种弹性波变网格有限差分正演模拟方法及其设备
CN115857001A (zh) 一种伪深度域声波同位网格逆时偏移成像方法以及设备
Mironov et al. Strong motion record processing of the Baikal rift zone

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