CN108108331B - 一种基于拟空间域弹性波方程的有限差分计算方法 - Google Patents

一种基于拟空间域弹性波方程的有限差分计算方法 Download PDF

Info

Publication number
CN108108331B
CN108108331B CN201711326455.2A CN201711326455A CN108108331B CN 108108331 B CN108108331 B CN 108108331B CN 201711326455 A CN201711326455 A CN 201711326455A CN 108108331 B CN108108331 B CN 108108331B
Authority
CN
China
Prior art keywords
spatial domain
quasi
wave
difference
point
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
CN201711326455.2A
Other languages
English (en)
Other versions
CN108108331A (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.)
Ocean University of China
National Deep Sea Center
Original Assignee
Ocean University of China
National Deep Sea Center
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 Ocean University of China, National Deep Sea Center filed Critical Ocean University of China
Priority to CN201711326455.2A priority Critical patent/CN108108331B/zh
Publication of CN108108331A publication Critical patent/CN108108331A/zh
Application granted granted Critical
Publication of CN108108331B publication Critical patent/CN108108331B/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
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Operations Research (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)

Abstract

一种基于拟空间域弹性波方程的有限差分计算方法，属于地震勘探领域，其主要思想是将常规弹性波方程变换为拟空间域弹性波方程，从而使空间域的等间隔的“距离”网格步长转化为可非等间隔的“传播时间”步长，由此人们可在严格定义速度模型的基础上精确计算速度界面两侧的“传播时间”。然后在此基础上给出了该方程及其完全匹配层边界条件的2N(N为正整数)阶精度有限差分表达式，基于此可实现逆时偏移过程中地震波的有限差分波场延拓。本发明方法可以很好地解决常规弹性波方程在逆时偏移剖面中速度界面形态畸变的问题；此外，基于拟空间域弹性波方程进行波场延拓可明显弱化界面假散射和层间反射波，从而进一步提高偏移成像的质量。

Description

一种基于拟空间域弹性波方程的有限差分计算方法
技术领域
本发明属于地震勘探领域，具体地涉及一种基于拟空间域弹性波方程的有限差分计算方法。
背景技术
地震勘探是目前探明地下地质构造最有效的地球物理方法，其主要包括地震信息采集、数据处理和成果解释三大环节。而地震数据处理是连接信息采集与成果解释环节的重要桥梁。长期以来，精确的地下构造偏移成像一直是地震数据处理中的研究热点。
在偏移成像领域，基于双程波动方程理论进行逆时偏移成像是公认的精确方法。该技术方法在当前的实现中普遍应用到重要的数值计算方法——有限差分法。众所周知，传统的波动方程有限差分计算方法本身存在着难以克服的固有问题——空间域的矩形网格剖分必然造成地震速度界面的畸变。图1给出了一个包含倾斜速度界面模型的常规有限差分法网格剖分示意图。图中黑色实线为倾斜速度界面，而黑色空心点和黑色实心点分别表示界面两侧的速度。显然，矩形网格剖分的结果使得原本光滑的倾斜速度界面畸变为明显的阶梯状折线，当地震波传播至这样的界面时会形成一系列的假散射；此外常规波动方程波场延拓产生的层间反射波会在偏移剖面上形成低频干扰和偏移假象，从而影响偏移剖面的精度。
为解决界面畸变的问题，一些学者采用空间变步长网格的思路，在介质变化剧烈的区域采用精细网格，在介质变化平缓的区域采用较粗网格剖分，但该方法仍未摆脱矩形网格的局限性；褚春雷等根据有限元的思想采用基于非规则三角网格剖分方法实现有限差分模拟，该方法可对弯曲的模型界面实现精细地描述，但计算量与传统矩形网格有限差分相比明显增加。为了减少层间反射波的影响，Baysal等在假设地下介质波阻抗恒定的情况下推导出无反射声波方程，其压制层间反射的效果仅在垂直入射时有较好效果；何兵寿等发展了任意广角声波方程逆时偏移技术，但该方法在浅层的成像效果较差；Yoon等将Poynting矢量成像条件引入到逆时偏移中实现不同方向波场的互相关成像，但在复杂构造区域应用中，Poynting矢量计算误差较大。
可以看出，目前的方法在解决界面畸变和压制层间反射波方面均存在局限性，而且在弹性波方程领域研究尚少。因此需要一种数值模拟方法，不但能够有效解决常规波动方程逆时偏移中的速度界面畸变及假散射问题，同时还可显著减少波场延拓中的层间反射波。
发明内容
本发明要解决的技术问题在于提供一种基于拟空间域弹性波方程的波场延拓有限差分数值计算方法。本发明首先提出了一种拟空间域弹性波方程，其主要思想是将常规变换到拟空间域后，空间域的等间隔的“距离”网格步长转化为可非等间隔的“传播时间”步长，由此人们可在严格定义速度模型的基础上精确计算速度界面两侧的“传播时间”。然后在此基础上给出了该方程及其完全匹配层(PML)边界条件的2N(N为正整数)阶精度有限差分表达式，基于此可实现逆时偏移过程中地震波的有限差分波场延拓。
本发明采取以下技术方案：
一种基于拟空间域弹性波方程的有限差分计算方法，它包括(1)对于实际地下连续的地质构造模型，根据单位网格长度将其按照矩形网格进行离散化，利用离散化网格点上的介质纵波和横波速度求出模型中各网格点处的纵波和横波拟空间域采样间隔，从而将空间域的弹性方程转化为拟空间域弹性波方程；(2)在每个拟空间域差分计算网格点上，求取拟空间域弹性波方程各一阶导数项 的拟空间域差分系数计算方程组中所需的一系列纵波或横波拟空间域传播时间，并由此计算各拟空间域差分计算网格点处的差分系数，进而得出各一阶导数项拟空间域2N阶精度差分表达式；(3)将步骤(2)得到的各一阶导数的2N阶精度差分表达式代入拟空间域弹性波方程中，得到整个拟空间域弹性波方程2N阶精度差分表达式；(4)在边界区域，定义拟空间域衰减因子和将衰减因子引入拟空间域弹性波方程，得到拟空间域弹性波方程完全匹配层边界条件，然后将各拟空间域阻尼项 及写成差分形式并与各拟空间域一阶导数项差分格式同时代入边界条件表达式，得到拟空间域弹性波方程完全匹配层(PML)边界条件2N阶精度有限差分格式；(5)根据观测系统参数，给定震源子波函数、炮集记录以及波场初始条件，基于拟空间域弹性波方程2N阶精度差分表达式和拟空间域弹性波方程完全匹配层边界条件的2N阶精度有限差分格式，沿正时方向或逆时方向递推计算各时刻的应力和速度分量的正时延拓波场值或逆时延拓波场值。
进一步，所述的步骤(1)为设空间域连续模型按照矩形网格离散化之后的单位网格长度为Δξ，ξ代表x或z，其中x和z分别表示空间的水平方向和垂直方向；离散化之后的网格模型中点坐标为(i,j)，i和j分别表示x和z方向上离散网格点坐标；
利用模型离散化之后各网格点对应的纵波速度cp求出地震波在长度为Δξ的网格上的纵波拟空间域采样间隔为Δτ(p)ξ，下标P表示纵波；τ(p)ξ代表τ(p)x或τ(p)z，其中τ(p)x和τ(p)z分别表示纵波拟空间域的水平方向和垂直方向；利用离散化之后各网格点对应的横波速度cs求出地震波在长度为Δξ的网格上的横波拟空间域采样间隔为Δτ(s)ξ，下标S表示横波；τ(s)ξ代表τ(s)x或τ(s)z，其中τ(s)x和τ(s)z分别表示横波拟空间域的水平方向和垂直方向；这里空间网格Δξ与纵横波拟空间域采样间隔Δτ(p)ξ和Δτ(s)ξ之间满足关系
由此可根据式(1)将应力波场σ，σ代表σp,σxxs,σzzs或σxz，以及速度分量波场vx和vz对空间域的导数变换为对拟空间域的导数，如式(2)所示；
式中cp表示介质的纵波速度，cs表示介质的横波速度，σ表示应力波场，σ代表σp,σxxs,σzzs或σxz，vx和vz分别表示水平和垂直的速度分量波场；
经过拟空间域坐标变换后，对于纵波拟空间域采样间隔，与空间域点(i,j)对应的拟空间域点的周围有四个纵波拟空间域采样间隔，这里统一记作其中l代表x或z；“-”代表坐标序号小于i或j的一侧；“+”代表坐标序号大于i或j的一侧；和分别代表(i-1,j)与(i,j)点之间的纵波拟空间域采样间隔和(i,j-1)与(i,j)点之间的纵波域拟空间域采样间隔；和分别代表(i+1,j)与(i,j)点之间的纵波拟空间域采样间隔和(i,j+1)与(i,j)点之间的纵波拟空间域采样间隔；对于横波拟空间域采样间隔，与空间域点(i,j)对应的拟空间域点的周围有四个横波拟空间域采样间隔，这里统一记作其中l代表x或z；“-”代表坐标序号小于i或j的一侧；“+”代表坐标序号大于i或j的一侧；和分别代表(i-1,j)与(i,j)点之间的横波拟空间域采样间隔和(i,j-1)与(i,j)点之间的横波拟空间域采样间隔；和分别代表(i+1,j)与(i,j)点之间的横波拟空间域采样间隔和(i,j+1)与(i,j)点之间的横波拟空间域采样间隔；
根据式(2)将常规空间域弹性波方程变换为拟空间域弹性波方程，其表达式如下：
上式中ρ表示介质的密度，t表示时间；在式(3)中，波动方程中不再有速度参数，即将常规空间域弹性波方程变换到拟空间域后，空间域的等间隔的“距离”网格步长转化为可非等间隔的“传播时间”步长，这里的“传播时间”步长即拟空间域采样间隔；
当速度界面与空间域相邻网格点之间单位网格线相交时，按照与速度界面的实际相交位置将单位网格线划分成对应不同纵波或横波速度的距离线段，在每个线段上根据对应的速度与线段长度分别计算传播时间，各线段“传播时间”之和即为该单位网格长度对应的纵波或横波拟空间域采样间隔。
进一步，所述的步骤(2)为定义与空间域网格模型中坐标为(i,j)的点相对应的拟空间域坐标为(τi,τj)；为了实现基于拟空间域弹性波方程的有限差分计算，需要对式(3)中各一阶拟空间域导数项进行差分表示，得到 的拟空间域2N阶精度差分表达式。
进一步，的拟空间域2N阶精度差分表达式的获得方法为设应力波场σp对变量τ(p)x具有2N阶导数，其中N为正整数，则应力波场σp对变量τ(p)x在τ(p)x＝τi+1/2处一阶导数的2N阶精度差分表达式为：
式中或表示纵波拟空间域τ(p)x方向上以τ(p)x＝τi+1/2点为中心点的不同网格点的坐标，m＝1,2,…,N-1,N，其中τi+1/2表示τi+1与τi之间的拟空间域半网格点位置；表示各阶差分系数，m＝1,2,…,2N-1,2N，通过求解式(5)所示的方程组可得到的值；
式中和表示拟空间域传播时间，m＝1,2,…,N-1,N，其可分别表示为
其中表示中心点(τi+1/2,τj)与变量σp所处网格点(τi+m,τj)之间的纵波拟空间域传播时间，表示中心点(τi+1/2,τj)与变量σp所处网格点(τi-(m-1),τj)之间的纵波拟空间域传播时间；
按上述类似的方式分别得到在(τi+1/2,τj)处的拟空间域2N阶精度差分表达式和差分系数m＝1,2,…,2N-1,2N；在(τi,τj+1/2)处的拟空间域2N阶精度差分表达式和差分系数m＝1,2,…,2N-1,2N； 在(τi,τj)处的拟空间域2N阶精度差分表达式和差分系数m＝1,2,…,2N-1,2N；在(τi+1/2,τj+1/2)处的拟空间域2N阶精度差分表达式和差分系数m＝1,2,…,2N-1,2N。
进一步，所述步骤(3)为将各应力波场σ，σ代表σp,σxxs,σzzs或σxz，以及水平和垂直的速度分量波场vx和vz关于时间一阶导数的时间域二阶精度差分表达式以及步骤(2)得到的各应力波场σ以及水平和垂直的速度分量波场vx和vz关于τ(p)x、τ(p)z、τ(s)x以及τ(s)z的一阶导数 的拟空间域2N阶精度差分表达式代入式(3)可得到
其中k表示离散的时间点时刻，t＝kΔt，Δt表示离散时间步长；表示拟空间域中(τi,τj)点处在k时刻时应力σp的波场值，同理可得其他应力分量波场 以及速度分量波场的含义；式(8)即为拟空间域弹性波方程的有限差分表达式。
进一步，所述的步骤(4)为在边界区域，定义拟空间域衰减因子和将衰减因子引入拟空间域弹性波方程，得到拟空间域弹性波方程的完全匹配层边界条件：
其中σpx和σpz表示应力σp在τx和τz方向上的分量，σxzx和σxzz表示应力σxz在τx和τz方向上的分量；vxx和vxz表示速度分量vx在τx和τz方向上的分量，vzx和vzz表示速度分量vz在τx和τz方向上的分量。和表示τx和τz方向上的衰减因子，其表达式为：
其中τm为拟空间域PML层中的点到中心波场外沿的垂直纵波拟空间域距离，R为拟空间域PML层的理论反射系数，其取值为10-5～10-7，τL为拟空间域PML层中最外层的点到中心波场外沿的垂直纵波拟空间域传播时间，在需要进行边界吸收的区域令或
将拟空间域阻尼项及写成差分形式并与各拟空间域一阶导数项差分格式同时代入边界条件表达式，得到拟空间弹性波方程完全匹配层边界条件2N阶精度有限差分格式：
式(11)为拟空间域弹性波方程关于速度分量的完全匹配层边界条件2N阶精度有限差分格式，式(12)为拟空间域弹性波方程中关于应力的完全匹配层边界条件2N阶精度有限差分格式。
进一步，所述的步骤(5)为给定观测系统参数、震源子波函数以及炮集记录，通过沿正时方向在每个k时刻将震源子波函数加载到相应的波场中，即可通过式(8)所示的拟空间域弹性波方程有限差分表达式递推求出全部时刻的应力波场σp,σxxs,σzzs，σxz以及速度分量波场vx和vz的正时延拓波场值；而通过沿逆时方向在每个k时刻将炮集记录加载到相应波场中，即可通过式(8)所示的拟空间域弹性波方程有限差分表达式递推求出全部时刻的应力波场σp,σxxs,σzzs，σxz以及速度分量波场vx和vz的逆时延拓波场值；在边界区域，基于(11)和(12)式实现拟空间域弹性波方程人工边界反射的吸收衰减。
本发明与现有技术相比的有益效果：
本发明提出了拟空间域弹性波方程及其有限差分数值计算方法。在拟空间域，速度界面两侧的“传播时间”是按其实际的速度与传播距离分别计算，而与该网格线对应的时间采样间隔为不同“传播时间”段之和。理论上讲，在拟空间域不再存在速度界面的畸变问题，甚至还会减弱相邻网格点之间模型参量的突变程度，从而可在偏移计算中期望降低假散射和界面反射。通过模型试验计算可得，与常规弹性波方程相比，拟空间域弹性波方程有限差分方法可以很好地解决常规弹性波方程在逆时偏移剖面中速度界面形态畸变的问题；此外，基于拟空间域弹性波方程进行波场延拓可明显弱化界面假散射和层间反射波，从而进一步提高偏移成像的质量。
附图说明
图1为常规有限差分法网格剖分示意图；
图2为含倾斜界面的两层速度模型(a)原光滑界面模型(b)10m网格间距的模型；
图3为合成炮集记录示例(第76炮)(a)水平速度分量(b)垂直速度分量；
图4为倾斜速度界面模型空间域网格剖分示意图；
图5为逆时偏移正时波场波前快照(第76炮，0.9s时刻)(a)基于常规弹性波方程(b)基于拟空间域弹性波方程；
图6为逆时偏移地震剖面(a)基于常规弹性波方程(b)基于拟空间域弹性波方程；
图7为逆时偏移地震剖面局部放大图示(a)基于常规弹性波方程(b)基于拟空间域弹性波方程。
具体实施方式
下面通过实施例结合附图来对本发明的技术方案作进一步解释，但本发明的保护范围不受实施例任何形式上的限制。
本发明提出了一种基于拟空间域弹性波方程的波场延拓有限差分数值计算方法。其主要思想是将常规弹性波方程变换到拟空间域后，空间域的等间隔的“距离”网格步长转化为可非等间隔的“传播时间”步长，由此人们可在严格定义速度模型的基础上精确计算速度界面两侧的“传播时间”。在此基础上，给出了该方程及其完全匹配层(PML)边界条件的2N(N为正整数)阶精度有限差分表达式，然后基于此可实现拟空间域弹性波波场的有限差分数值计算。
本发明的主要实施过程分为五个步骤：(1)对于实际地下连续的地质构造模型，根据单位网格长度将其按照矩形网格进行离散化，利用离散化网格点上的介质纵波和横波速度求出模型中各网格点处的纵波和横波拟空间域采样间隔，从而将空间域的弹性方程转化为拟空间域弹性波方程；(2)在每个拟空间域差分计算网格点上，求取拟空间域弹性波方程各一阶导数项 的拟空间域差分系数计算方程组中所需的一系列纵波和横波拟空间域传播时间，并由此计算各拟空间域差分计算网格点处的差分系数，进而得出各一阶导数项拟空间域2N阶精度差分表达式；(3)将步骤(2)得到的各一阶导数的2N阶精度差分表达式代入拟空间域弹性波方程中，得到整个拟空间域弹性波方程2N阶精度差分表达式；(4)在边界区域，定义拟空间域衰减因子和将衰减因子引入拟空间域弹性波方程，得到拟空间域弹性波方程完全匹配层边界条件，然后将各拟空间域阻尼项 及写成差分形式并与各拟空间域一阶导数项差分格式同时代入边界条件表达式，得到拟空间域弹性波方程完全匹配层(PML)边界条件2N阶精度有限差分格式；(5)根据观测系统参数，给定震源子波函数、炮集记录以及波场初始条件，基于拟空间域弹性波方程2N阶精度差分表达式和拟空间域弹性波方程完全匹配层边界条件的2N阶精度有限差分格式，沿正时方向或逆时方向递推计算各时刻的应力和速度分量的正时延拓波场值或逆时延拓波场值。
实施例1
本实验的主要目的是验证拟空间域弹性波方程有限差分计算方法在逆时偏移中解决速度界面畸变以及压制界面假散射和层间反射波方面的有效性。
实验采用一个含有倾斜界面的两层速度模型，如图2(a)所示，其横向和纵向长度分别为4000m和2000m，界面上下两侧纵波速度分别为2500m/s和3500m/s，密度为2000kg/m3。采用纵波震源激发，而模型横波速度由纵波速度和泊松比值计算得到(本实验所用泊松比为0.25，密度为2000kg/m3)。将此倾斜界面模型以10m的纵横向网格间距剖分后所得网格模型如图2(b)所示，可见原光滑的速度界面已变为明显的阶梯状界面(如图中白色箭头所指位置)。实验中建立道固定、炮移动的观测系统，炮点位于500m至3480m之间，炮间隔为20m，共150炮；每炮接收道数为401道，各道位于0m至4000m之间，道间隔为10m；炮点与接收点深度均为10m。
显然，要检验一个偏移方法在解决速度界面畸变方面的有效性，需要保证采集的地震记录是精确的。本实验所需的炮集记录是采用传统的有限差分法波动方程正演模拟获得。理论上讲，只有采用足够小的网格间距才能保证所获得的炮集记录相对精确。为此，本文先采用纵横向网格间距均为1m网格模型进行正演模拟(注：即使网格点数仅增加一倍也将造成计算量的巨增，因此在实际处理中无论正演模拟还是偏移反演，一般均难以采用如此小的网格间距)。模拟采用主频为35Hz的Ricker子波，差分精度为时间二阶空间十六阶，共获得150炮的地震记录，其中第76炮的记录如图3所示。
现基于如图2(b)所示的10m网格间距模型，采用差分精度为时间二阶空间/拟空间十六阶的拟空间域弹性波方程有限差分算法进行逆时偏移计算(注：拟空间域的逆时偏移需要加入速度界面信息计算两网格点之间的旅行时步长)。
根据逆时偏移的实现过程，其主要包括两个部分，首先基于拟空间域弹性波方程有限差分计算得到地震波传播的正时波场和逆时波场，然后依据互相关成像条件对波场计算求得地下偏移成像点的成像值。其中第一部分所涉及地计算地震波传播波场的过程具体步骤如下：(1)根据图2(b)所示的网格模型，空间域离散化之后的单位网格长度为Δξ＝10m，ξ可以代表x或z，其中x和z分别表示空间的水平方向和垂直方向；离散化之后的网格模型中点坐标为(i,j)，i和j分别表示x和z方向上离散网格点坐标；
利用模型离散化之后各网格点对应的纵波速度cp求出地震波在长度为Δξ的网格上的纵波拟空间域采样间隔为Δτ(p)ξ，下标P表示纵波；τ(p)ξ代表τ(p)x或τ(p)z，其中τ(p)x和τ(p)z分别表示纵波拟空间域的水平方向和垂直方向；利用离散化之后各网格点对应的横波速度cs求出地震波在长度为Δξ的网格上的横波拟空间域采样间隔为Δτ(s)ξ，下标S表示横波；τ(s)ξ代表τ(s)x或τ(s)z，其中τ(s)x和τ(s)z分别表示横波拟空间域的水平方向和垂直方向；这里空间网格Δξ与纵横波拟空间域采样间隔Δτ(p)ξ和Δτ(s)ξ之间满足关系
由此可根据式(1)将应力σ，σ代表σp,σxxs,σzzs或σxz，以及和速度分量vx和vz对空间域的导数变换为对拟空间域的导数，如式(2)所示；
式中cp表示介质的纵波速度，cs表示介质的横波速度，σ表示应力波场，σ代表σp,σxxs,σzzs或σxz，vx和vz分别表示水平和垂直的速度分量波场；
经过拟空间域坐标变换后，对于纵波拟空间域采样间隔，与空间域点(i,j)对应的拟空间域点的周围有四个纵波拟空间域采样间隔，这里统一记作其中l代表x或z；“-”代表坐标序号小于i或j的一侧；“+”代表坐标序号大于i或j的一侧；和分别代表(i-1,j)与(i,j)点之间的纵波拟空间域采样间隔和(i,j-1)与(i,j)点之间的纵波拟空间域采样间隔；和分别代表(i+1,j)与(i,j)点之间的纵波拟空间域采样间隔和(i,j+1)与(i,j)点之间的纵波拟空间域采样间隔；对于横波拟空间域采样间隔，与空间域点(i,j)对应的拟空间域点的周围有四个横波拟空间域采样间隔，这里统一记作其中l代表x或z；“-”代表坐标序号小于i或j的一侧；“+”代表坐标序号大于i或j的一侧；和分别代表(i-1,j)与(i,j)点之间的横波拟空间域采样间隔和(i,j-1)与(i,j)点之间的横波拟空间域采样间隔；和分别代表(i+1,j)与(i,j)点之间的横波拟空间域采样间隔和(i,j+1)与(i,j)点之间的横波拟空间域采样间隔；
根据式(2)可将常规空间域弹性波方程变换为拟空间域弹性波方程，其表达式如下：
上式中ρ表示介质的密度，t表示时间；
以图4为例说明纵波拟空间域采样间隔计算的过程，在图4中P和Q点为速度模型按照矩形网格剖分之后两相邻空间网格点，黑色实线所示的速度界面与网格线PQ相交于B点，此时P和Q点之间的纵波拟空间域采样间隔为Δτ(p)ξ＝ΔτPB+ΔτBQ，其中ΔτPB表示线段PB之间的纵波旅行时，ΔτBQ表示线段BQ之间的纵波旅行时。同理可推导出横波拟空间域采样间隔的计算过程。
(2)定义与空间域网格模型中坐标为(i,j)的点相对应的拟空间域坐标为(τi,τj)。为了实现基于拟空间域弹性波方程的有限差分计算，需要对式(3)中各一阶拟空间域导数项进行差分表示，得到 的拟空间域2N阶精度差分表达式；
以为例，设应力σp对变量τ(p)x具有2N阶导数，其中N为正整数，则应力波场σp对变量τ(p)x在τ(p)x＝τi+1/2处一阶导数的2N阶精度差分表达式为：
式中或表示纵波拟空间域τ(p)x方向上以τ(p)x＝τi+1/2点为中心点的不同网格点的坐标，m＝1,2,…,N-1,N，其中τi+1/2表示τi+1与τi之间的拟空间域半网格点位置；表示各阶差分系数，m＝1,2,…,2N-1,2N，通过求解式(5)所示的方程组可得到的值。
式中和表示拟空间域传播时间，m＝1,2,…,N-1,N，其可分别表示为
其中表示中心点(τi+1/2,τj)与变量σp所处网格点(τi+m,τj)之间的纵波拟空间域传播时间，表示中心点(τi+1/2,τj)与变量σp所处网格点(τi-(m-1),τj)之间的纵波拟空间域传播时间。
可以看出，差分系数的解取决于拟空间域传播时间和(m＝1,2,…,N-1,N)。当空间网格速度不变时，则此时拟空间域为均匀步长，有(k＝1,2,…,N-1,N)；当空间网格速度变化时，存在关系即此时拟空间域为变步长，满足(k＝1,2,…,N-1,N)。由于拟空间域传播时间与网格速度和空间网格大小有关，因此即使在差分阶数相同的条件下，不同的速度和空间网格大小对应的差分系数也不同(表1和表2分别为空间网格点速度为3000m/s、空间网格步长为5m和空间网格点速度为4000m/s、空间网格步长为10m时的各阶精度差分系数)。
表1 不同精度的拟空间域差分系数(速度v＝3000m/s，空间网格步长Δξ＝5m)
表2 不同精度的拟空间域差分系数(速度v＝4000m/s，空间网格步长Δξ＝10m)
同理，可以上述类似的方式分别得到在(τi+1/2,τj)处的拟空间域2N阶精度差分表达式和差分系数m＝1,2,…,2N-1,2N；在(τi,τj+1/2)处的拟空间域2N阶精度差分表达式和差分系数m＝1,2,…,2N-1,2N； 在(τi,τj)处的拟空间域2N阶精度差分表达式和差分系数 m＝1,2,…,2N-1,2N；在(τi+1/2,τj+1/2)处的拟空间域2N阶精度差分表达式和差分系数m＝1,2,…,2N-1,2N；
(3)将各应力波场σ，σ代表σp,σxxs,σzzs或σxz，以及水平和垂直的速度分量波场vx和vz关于时间一阶导数的时间域二阶精度差分表达式以及步骤(2)得到的各应力波场σ以及水平和垂直的速度分量波场vx和vz关于τ(p)x、τ(p)z、τ(s)x以及τ(s)z的一阶导数 的拟空间域2N阶精度差分表达式代入式(3)可得到
其中k表示离散的时间点时刻，t＝kΔt(Δt表示离散时间步长)；表示拟空间域中(τi,τj)点处在k时刻时应力σp的波场值，同理可得其他应力分量波场 以及速度分量波场的含义；式(8)即为拟空间域弹性波方程的有限差分表达式。
(4)在边界区域，定义拟空间域衰减因子和将衰减因子引入拟空间域弹性波方程，得到拟空间域弹性波方程的完全匹配层边界条件，如式(9)所示，式中和表示τx和τz方向上的衰减因子，其表达式为：
其中τm为拟空间域PML层中的点到中心波场外沿的垂直纵波拟空间域传播时间，R为拟空间域PML层的理论反射系数，其取值为10-5～10-7，τL为拟空间域PML层中最外层的点到中心波场外沿的垂直纵波拟空间域传播时间，在需要进行边界吸收的区域令或
将拟空间域阻尼项及写成差分形式并与各拟空间域一阶导数项差分格式同时代入边界条件表达式，得到拟空间弹性波方程完全匹配层边界条件2N阶精度有限差分格式，如式(11)和(12)所示，其中式(11)为拟空间域弹性波方程关于速度分量的完全匹配层边界条件2N阶精度有限差分格式，式(12)为拟空间域弹性波方程中关于应力的完全匹配层边界条件2N阶精度有限差分格式。
(5)根据前面给定的观测系统参数，这里利用35Hz的Ricker子波和图3所示的炮集记录，通过沿正时间方向在每个k时刻将震源子波函数加载到相应的波场中，即可通过式(8)所示的拟空间域弹性波方程有限差分表达式递推求出全部时刻的应力波场σp,σxxs,σzzs，σxz以及速度分量波场vx和vz的正时延拓波场值；而通过沿逆时间方向在每个k时刻将炮集记录加载到相应波场中，即可通过式(8)所示的拟空间域弹性波方程有限差分表达式递推求出全部时刻的应力波场σp,σxxs,σzzs，σxz以及速度分量波场vx和vz的逆时延拓波场值；在边界区域，基于式(11)和式(12)实现拟空间域弹性波方程人工边界反射的吸收衰减。0.9s时刻、第76炮对应的应力波场P的正时延拓波场如图5(a)所示。
以上即为基于拟空间域弹性波方程的有限差分差分方法实现逆时偏移中正时延拓波场和逆时延拓波场计算的具体过程。基于由拟空间域弹性波方程有限差分计算方法所得到的地震波波场，再应用归一化互相关成像条件，可以得到逆时偏移剖面，如图6(b)所示。为更直观的对比本发明方法与常规有限差分方法得到的偏移剖面中倾斜界面的形态，基于如图2(b)所示的10m网格间距模型，采用差分精度为时间二阶空间十六阶的常规弹性波方程有限差分算法进行逆时偏移，得到逆时偏移剖面如图6(a)所示。对图6(a,b)中白色椭圆区域内界面同相轴进行放大显示，如图7(a,b)所示。可以看出，常规弹性波方程逆时偏移剖面中倾斜界面形态(图7(a)中灰色虚线)相比于真实界面形态(图7(a)灰色实线)出现明显畸变，而拟空间域弹性波方程逆时偏移剖面中倾斜界面的形态则与真实界面形态(图7(b)中灰色实线)基本吻合。由此可见，本发明所述方法可有效地解决逆时偏移中速度界面的畸变问题。
为验证拟空间域弹性波方程在压制界面假散射和层间反射波方面的有效性，这里基于如图2(b)所示的10m网格间距模型，采用差分精度为时间二阶空间十六阶的常规弹性波方程有限差分算法得到的0.9s时刻、第76炮的正时延拓波场的波前快照如图5(a)所示。从图5(a,b)可以看出，常规弹性波方程的波场中存在明显的界面假散射(如图5(a)中椭圆区域所示)，而基于拟空间域弹性波方程的波场中并无明显的界面假散射(如图5(b)中椭圆内区域所示)。通过对比图5(a)和图5(b)中箭头处界面反射波可以看出，基于拟空间域弹性波方程的有限差分计算方法对界面反射波(尤其是近垂直入射的反射波)具有明显压制效果。由此可见，本发明所述方法在压制界面假散射和层间反射波的有效性。

Claims (6)

1.一种基于拟空间域弹性波方程的有限差分计算方法，其特征在于它包括(1)对于实际地下连续的地质构造模型，根据单位网格长度将其按照矩形网格进行离散化，利用离散化网格点上的介质纵波和横波速度求出模型中各网格点处的纵波和横波拟空间域采样间隔，从而将空间域的弹性方程转化为拟空间域弹性波方程；(2)在每个拟空间域差分计算网格点上，求取拟空间域弹性波方程各一阶导数项 的拟空间域差分系数计算方程组中所需的一系列纵波或横波拟空间域传播时间，并由此计算各拟空间域差分计算网格点处的差分系数，进而得出各一阶导数项拟空间域2N阶精度差分表达式；(3)将步骤(2)得到的各一阶导数的2N阶精度差分表达式代入拟空间域弹性波方程中，得到整个拟空间域弹性波方程2N阶精度差分表达式；(4)在边界区域，定义拟空间域衰减因子和将衰减因子引入拟空间域弹性波方程，得到拟空间域弹性波方程完全匹配层边界条件，完全匹配层简写为PML，然后将各拟空间域阻尼项及写成差分形式并与各拟空间域一阶导数项差分格式同时代入边界条件表达式，得到拟空间域弹性波方程完全匹配层边界条件2N阶精度有限差分格式；(5)根据观测系统参数，给定震源子波函数、炮集记录以及波场初始条件，基于拟空间域弹性波方程2N阶精度差分表达式和拟空间域弹性波方程完全匹配层边界条件的2N阶精度有限差分格式，沿正时方向或逆时方向递推计算各时刻的应力和速度分量的正时延拓波场值或逆时延拓波场值；
所述的步骤(1)为设空间域连续模型按照矩形网格离散化之后的单位网格长度为△ξ，ξ代表x或z，其中x和z分别表示空间的水平方向和垂直方向；离散化之后的网格模型中点坐标为(i,j)，i和j分别表示x和z方向上离散网格点坐标；
利用模型离散化之后各网格点对应的纵波速度cp求出地震波在长度为△ξ的网格上的纵波拟空间域采样间隔为△τ(p)ξ，下标P表示纵波；τ(p)ξ代表τ(p)x或τ(p)z，其中τ(p)x和τ(p)z分别表示纵波拟空间域的水平方向和垂直方向；利用离散化之后各网格点对应的横波速度cs求出地震波在长度为△ξ的网格上的横波拟空间域采样间隔为△τ(s)ξ，下标S表示横波；τ(s)ξ代表τ(s)x或τ(s)z，其中τ(s)x和τ(s)z分别表示横波拟空间域的水平方向和垂直方向；这里空间网格△ξ与纵横波拟空间域采样间隔△τ(p)ξ和△τ(s)ξ之间满足关系
由此可根据式(1)将应力波场σ，σ代表σp,σxxs,σzzs或σxz，以及速度分量波场vx和vz对空间域的导数变换为对拟空间域的导数，如式(2)所示；
式中cp表示介质的纵波速度，cs表示介质的横波速度，σ表示应力波场，vx和vz分别表示水平和垂直速度分量波场；
经过拟空间域坐标变换后，对于纵波拟空间域采样间隔，与空间域点(i,j)对应的拟空间域点的周围有四个纵波拟空间域采样间隔，这里统一记作其中l代表i或j；“-”代表坐标序号小于i或j的一侧；“+”代表坐标序号大于i或j的一侧；和分别代表(i-1,j)与(i,j)点之间的纵波拟空间采样间隔和(i,j-1)与(i,j)点之间的纵波拟空间采样间隔；和分别代表(i+1,j)与(i,j)点之间的纵波拟空间采样间隔和(i,j+1)与(i,j)点之间的纵波拟空间采样间隔；对于横波拟空间域采样间隔，与空间域点(i,j)对应的拟空间域点的周围有四个横波拟空间域采样间隔，这里统一记作其中l代表i或j；“-”代表坐标序号小于i或j的一侧；“+”代表坐标序号大于i或j的一侧；和分别代表(i-1,j)与(i,j)点之间的横波拟空间采样间隔和(i,j-1)与(i,j)点之间的横波拟空间采样间隔；和分别代表(i+1,j)与(i,j)点之间的横波拟空间采样间隔和(i,j+1)与(i,j)点之间的横波拟空间采样间隔；
根据式(2)将常规空间域弹性波方程变换为拟空间域弹性波方程，其表达式如下：
上式中ρ表示介质的密度，t表示时间；在式(3)中，波动方程中不再有速度参数，即将常规空间域弹性波方程变换到拟空间域后，空间域的等间隔的“距离”网格步长转化为可非等间隔的“传播时间”步长，这里的“传播时间”步长即拟空间域采样间隔；
当速度界面与空间域相邻网格点之间单位网格线相交时，按照与速度界面的实际相交位置将单位网格线划分成对应不同纵波或横波速度的距离线段，在每个线段上根据对应的速度与线段长度分别计算传播时间，各线段“传播时间”之和即为该单位网格长度对应的纵波或横波拟空间域采样间隔。
2.据权利要求1所述的方法，其特征在于所述的步骤(2)为定义与空间域网格模型中坐标为(i,j)的点相对应的拟空间域坐标为(τi,τj)；为了实现基于拟空间域弹性波方程的有限差分计算，需要对式(3)中各一阶拟空间域导数项进行差分表示，得到 的拟空间域2N阶精度差分表达式。
3.根据权利要求2所述的方法，其特征在于的拟空间域2N阶精度差分表达式的获得方法为设应力波场σp对变量τ(p)x具有2N阶导数，其中N为正整数，则应力波场σp对变量τ(p)x在τ(p)x＝τi+1/2处一阶导数的2N阶精度差分表达式为：
式中或表示纵波拟空间域τ(p)x方向上以τ(p)x＝τi+1/2点为中心点的不同网格点的坐标，m＝1,2,…,N-1,N，其中τi+1/2表示τi+1与τi之间的拟空间域半网格点位置；表示各阶差分系数，m＝1,2,…,2N-1,2N，通过求解式(5)所示的方程组可得到的值；
式中和表示拟空间域传播时间，m＝1,2,…,N-1,N，其可分别表示为
其中表示中心点(τi+1/2,τj)与变量σp所处网格点(τi+m,τj)之间的纵波拟空间域传播时间，表示中心点(τi+1/2,τj)与变量σp所处网格点(τi-(m-1),τj)之间的纵波拟空间域传播时间；
按上述的方式分别得到在(τi+1/2,τj)处的拟空间域2N阶精度差分表达式和差分系数m＝1,2,…,2N-1,2N；在(τi,τj+1/2)处的拟空间域2N阶精度差分表达式和差分系数m＝1,2,…,2N-1,2N； 在(τi,τj)处的拟空间域2N阶精度差分表达式和差分系数m＝1,2,…,2N-1,2N；在(τi+1/2,τj+1/2)处的拟空间域2N阶精度差分表达式和差分系数m＝1,2,…,2N-1,2N。
4.根据权利要求1所述的方法，其特征在于所述步骤(3)为将各应力波场σ，σ代表σp,σxxs,σzzs或σxz，以及水平和垂直速度分量波场vx和vz关于时间一阶导数 的时间域二阶精度差分表达式以及步骤(2)得到的各应力波场σ以及水平和垂直的速度分量波场vx和vz关于τ(p)x、τ(p)z、τ(s)x以及τ(s)z的一阶导数 的拟空间域2N阶精度差分表达式代入式(3)可得到
其中k表示离散的时间点时刻，ρ表示介质的密度，t＝k△t，△t表示离散时间步长；表示拟空间域中(τi,τj)点处在k时刻时应力σp的波场值，同理可得其他应力分量波场 以及速度分量波场 的含义；式(8)即为拟空间域弹性波方程的有限差分表达式。
5.根据权利要求4所述的方法，其特征在于所述的步骤(4)为在边界区域，定义拟空间域衰减因子和将衰减因子引入拟空间域弹性波方程，得到拟空间域弹性波方程的完全匹配层边界条件：
其中σpx和σpz表示应力σp在τx和τz方向上的分量，σxzx和σxzz表示应力σxz在τx和τz方向上的分量；vxx和vxz表示速度分量vx在τx和τz方向上的分量，vzx和vzz表示速度分量vz在τx和τz方向上的分量，和表示τx和τz方向上的衰减因子，其表达式为：
其中τm为拟空间域PML层中的点到中心波场外沿的垂直纵波拟空间域传播时间，R为拟空间域PML层的理论反射系数，其取值为10-5～10-7，τL为拟空间域PML层中最外层的点到中心波场外沿的垂直纵波拟空间域传播时间，在需要进行边界吸收的区域令或
将拟空间域阻尼项及写成差分形式并与各拟空间域一阶导数项差分格式同时代入边界条件表达式，得到拟空间弹性波方程完全匹配层边界条件2N阶精度有限差分格式：
式(11)为拟空间域弹性波方程关于速度分量的完全匹配层边界条件2N阶精度有限差分格式，式(12)为拟空间域弹性波方程中关于应力的完全匹配层边界条件2N阶精度有限差分格式。
6.根据权利要求5所述的方法，其特征在于所述的步骤(5)为给定观测系统参数、震源子波函数以及炮集记录，通过沿正时方向在每个k时刻将震源子波函数加载到相应的波场中，即可通过式(8)所示的拟空间域弹性波方程有限差分表达式递推求出全部时刻的应力波场σp,σxxs,σzzs，σxz以及速度分量波场vx和vz的正时延拓波场值；而通过沿逆时方向在每个k时刻将炮集记录加载到相应波场中，即可通过式(8)所示的拟空间域弹性波方程有限差分表达式递推求出全部时刻的应力波场σp,σxxs,σzzs，σxz以及速度分量波场vx和vz的逆时延拓波场值；在边界区域，基于(11)和(12)式实现拟空间域弹性波方程人工边界反射的吸收衰减。
CN201711326455.2A 2017-12-13 2017-12-13 一种基于拟空间域弹性波方程的有限差分计算方法 Active CN108108331B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711326455.2A CN108108331B (zh) 2017-12-13 2017-12-13 一种基于拟空间域弹性波方程的有限差分计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711326455.2A CN108108331B (zh) 2017-12-13 2017-12-13 一种基于拟空间域弹性波方程的有限差分计算方法

Publications (2)

Publication Number Publication Date
CN108108331A CN108108331A (zh) 2018-06-01
CN108108331B true CN108108331B (zh) 2018-11-02

Family

ID=62215741

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711326455.2A Active CN108108331B (zh) 2017-12-13 2017-12-13 一种基于拟空间域弹性波方程的有限差分计算方法

Country Status (1)

Country Link
CN (1) CN108108331B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108983285B (zh) * 2018-07-19 2019-12-13 中国石油大学(北京) 一种基于矩张量的多种震源波场模拟方法及装置
CN109782338B (zh) * 2019-01-03 2020-06-30 深圳信息职业技术学院 一种频率域地震正演模拟的高精度差分数值法
CN113139266B (zh) * 2020-01-18 2024-05-17 中国科学院地质与地球物理研究所 纵、横波数值模拟方法及系统
CN111257930B (zh) * 2020-02-18 2022-03-29 中国石油大学(华东) 一种黏弹各向异性双相介质区域变网格求解算子
CN112526605B (zh) * 2020-12-24 2022-09-02 广州海洋地质调查局 一种采用地震数值模拟勘探天然气水合物的方法
CN115081267B (zh) * 2022-05-18 2023-08-29 内蒙古农业大学 基于优化的时空变步长有限差分地震波数值模拟方法
CN115453621B (zh) * 2022-09-14 2024-03-22 成都理工大学 基于一阶速度-应力方程的纵横波解耦分离假象去除方法
CN116822297B (zh) * 2023-06-30 2024-01-16 哈尔滨工程大学 一种应用于弹性波正演的三阶Higdon阻尼吸收边界方法

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
NO322089B1 (no) * 2003-04-09 2006-08-14 Norsar V Daglig Leder Fremgangsmate for simulering av lokale prestakk dypmigrerte seismiske bilder
CN101021568A (zh) * 2007-02-07 2007-08-22 匡斌 基于最大能量旅行时计算的三维积分叠前深度偏移方法
US9348044B2 (en) * 2012-04-19 2016-05-24 Cgg Services Sa Vectorization of fast fourier transform for elastic wave propogation for use in seismic underwater exploration of geographical areas of interest
CN103149585B (zh) * 2013-01-30 2016-02-17 中国石油天然气集团公司 一种弹性偏移地震波场构建方法及装置
WO2016108896A1 (en) * 2014-12-31 2016-07-07 Landmark Graphics Corporation Seismic elastic wave simulation for tilted transversely isotropic media using adaptive lebedev staggered grid
CN107315192B (zh) * 2016-04-26 2019-07-05 中国石油化工股份有限公司 基于二维各向同性介质的弹性波波场数值的模拟方法
CN106597535A (zh) * 2016-12-01 2017-04-26 桂林理工大学 一种提高弹性波逆时偏移计算率和空间分辨率的方法

Also Published As

Publication number Publication date
CN108108331A (zh) 2018-06-01

Similar Documents

Publication Publication Date Title
CN108108331B (zh) 一种基于拟空间域弹性波方程的有限差分计算方法
CN108051855B (zh) 一种基于拟空间域声波方程的有限差分计算方法
Evans et al. Teleseismic velocity tomography using the ACH method: theory and application to continental-scale studies
US8560242B2 (en) Pseudo-analytical method for the solution of wave equations
Day FINITE-ELEMENT ANALYSIS OF SEISMIC SCATTERING PROBLEMS.
Furumura et al. Subduction zone guided waves and the heterogeneity structure of the subducted plate: Intensity anomalies in northern Japan
US9075163B2 (en) Interferometric seismic data processing
US20070233437A1 (en) Method of Evaluating the Interaction Between a Wavefield and a Solid Body
CN106932819A (zh) 基于各向异性马尔科夫随机域的叠前地震参数反演方法
CN104820242B (zh) 一种面向叠前反演的道集振幅分频补偿方法
CN110007340A (zh) 基于角度域直接包络反演的盐丘速度密度估计方法
Keers et al. Viscoacoustic crosswell imaging using asymptotic waveforms
Yin et al. Improving horizontal resolution of high-frequency surface-wave methods using travel-time tomography
CN115600373A (zh) 黏滞各向异性介质qP波模拟方法、系统、设备及应用
CN110780341B (zh) 一种各向异性地震成像方法
CN108363097A (zh) 一种地震资料偏移成像方法
Cormier et al. Calculation of strong ground motion due to an extended earthquake source in a laterally varying structure
Zhao et al. A new scheme for fast calculation of seismic traveltimes and ray paths in heterogeneous media
CN110161561A (zh) 一种油气储层中的可控层位分阶层间多次波模拟方法
Kuhn A numerical study of Lamb's problem
Wei et al. A high-resolution microseismic source location method based on contrast source algorithm
Shin et al. Laplace-domain full waveform inversion using irregular finite elements for complex foothill environments
CN110609325B (zh) 弹性波场数值模拟方法及系统
Stead Finite differences and a coupled analytic technique with applications to explosions and earthquakes
Jansen et al. Acoustic tomography in solids using a bent ray SIRT algorithm

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