一种基于拟空间域声波方程的有限差分计算方法与流程

本发明属于地震勘探领域，具体地涉及一种基于拟空间域声波方程的有限差分计算方法。

背景技术：

地震勘探是目前探明地下地质构造最有效的地球物理方法，其主要包括地震信息采集、数据处理和成果解释三大环节。而地震数据处理是连接信息采集与成果解释环节的重要桥梁。长期以来，精确的地下构造偏移成像一直是地震数据处理中的研究热点。

在偏移成像领域，基于双程波动方程理论进行逆时偏移成像是公认的精确方法。该技术方法在当前的实现中普遍应用到重要的数值计算方法——有限差分法。众所周知，传统的波动方程有限差分计算方法本身存在着难以克服的固有问题——空间域的矩形网格剖分必然造成地震速度界面的畸变。图1给出了一个包含倾斜速度界面模型的常规有限差分法网格剖分示意图。图中黑色实线为倾斜速度界面，而黑色空心点和黑色实心点分别表示界面两侧的速度。显然，矩形网格剖分的结果使得原本光滑的倾斜速度界面畸变为明显的阶梯状折线，当地震波传播至这样的界面时会形成一系列的假散射；此外常规波动方程波场延拓产生的层间反射波会在偏移剖面上形成低频干扰和偏移假象，从而影响偏移剖面的精度。

为解决界面畸变的问题，一些学者采用空间变步长网格的思路，在介质变化剧烈的区域采用精细网格，在介质变化平缓的区域采用较粗网格剖分，但该方法仍未摆脱矩形网格的局限性；褚春雷等根据有限元的思想采用基于非规则三角网格剖分方法实现有限差分模拟，该方法可对弯曲的模型界面实现精细地描述，但计算量与传统矩形网格有限差分相比明显增加。为了减少层间反射波的影响，baysal等在假设地下介质波阻抗恒定的情况下推导出无反射声波方程，其压制层间反射的效果仅在垂直入射时有较好效果；何兵寿等发展了任意广角声波方程逆时偏移技术，但该方法在浅层的成像效果较差；yoon等将poynting矢量成像条件引入到逆时偏移中实现不同方向波场的互相关成像，但在复杂构造区域应用中，poynting矢量计算误差较大。

可以看出，目前的方法在解决界面畸变和压制层间反射波方面均存在局限性。因此需要一种数值模拟方法，不但能够有效解决常规波动方程逆时偏移中的速度界面畸变及假散射问题，同时还可显著减少波场延拓中的层间反射波。

技术实现要素：

本发明要解决的技术问题在于提供一种基于拟空间域声波方程的波场延拓有限差分数值计算方法。本发明首先提出了一种拟空间域声波方程，其主要思想是将常规声波方程变换到拟空间域后，空间域的等间隔的“距离”网格步长转化为可非等间隔的“传播时间”步长，由此人们可在严格定义速度模型的基础上精确计算速度界面两侧的“传播时间”。然后在此基础上给出了该方程及其完全匹配层(pml)边界条件的2n(n为正整数)阶精度有限差分表达式，基于此可实现逆时偏移过程中地震波的有限差分波场延拓。

本发明采取以下技术方案：

一种基于拟空间域声波方程的有限差分计算方法，它包括(1)对于实际地下连续的地质构造模型，根据单位网格长度将其按照矩形网格进行离散化，利用离散化网格点上的介质速度求出模型中各网格点处的拟空间域采样间隔，从而将空间域的声波方程转化为拟空间域声波方程；(2)在每个拟空间域差分计算网格点上，求取拟空间域声波方程中各一阶导数项和的拟空间域差分系数计算方程组中所需的一系列拟空间域传播时间，并由此计算各拟空间域差分计算网格点处的差分系数，进而得出各一阶导数项拟空间域2n阶精度差分表达式；(3)在中心波场区域，将各一阶导数项拟空间域2n阶精度差分表达式代入拟空间域声波方程中，得到整个拟空间域声波方程2n阶精度差分表达式；(4)在边界区域，定义拟空间域衰减因子和将衰减因子引入拟空间域声波方程，得到拟空间域声波方程完全匹配层边界条件，然后将各拟空间域阻尼项和写成差分形式并与各拟空间域一阶导数项差分格式同时代入边界条件表达式，得到拟空间域声波方程完全匹配层(pml)边界条件2n阶精度有限差分格式；(5)根据观测系统参数，给定震源子波函数、炮集记录以及波场初始条件，基于拟空间域声波方程2n阶精度差分表达式和拟空间域声波方程完全匹配层边界条件的2n阶精度有限差分格式，沿正时方向或逆时方向递推计算各时刻的正时延拓波场值或逆时延拓波场值。

进一步，所述的步骤(1)为设空间域连续模型按照矩形网格离散化之后的单位网格长度为△ξ，ξ代表x或z，其中x和z分别表示空间的水平方向和垂直方向；离散化之后的网格模型中点坐标为(i,j)，i和j分别表示x和z方向上离散网格点坐标；

利用离散化之后各网格点处对应的纵波速度cp求出地震波在长度为△ξ的网格上的拟空间域采样间隔为△τξ，τξ代表τx或τz，其中τx和τz分别表示拟空间域的水平方向和垂直方向，这里空间网格△ξ与拟空间域采样间隔△τξ之间满足关系△ξ＝cp△τξ；由此可根据式(1)将应力波场p和速度分量波场vx和vz对空间域的导数变换为对拟空间域的导数；

式中cp表示介质的纵波速度，p表示应力波场，τξ代表τx或τz，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)点之间的拟空间域采样间隔；

根据式(1)将常规空间域声波方程变换为拟空间域声波方程，其表达式如下：

上式中ρ表示介质的密度，t表示时间。在式(2)中，波动方程中不再有速度参数，即将常规空间域声波方程变换到拟空间域后，空间域的等间隔的“距离”网格步长转化为可非等间隔的“传播时间”步长，这里的“传播时间”步长即拟空间域采样间隔；

当速度界面与空间域相邻网格点之间单位网格线相交时，按照与速度界面的实际相交位置将单位网格线划分成对应不同速度的距离线段，在每个线段上根据对应的速度与线段长度分别计算传播时间，各线段“传播时间”之和即为该单位网格长度对应的拟空间域采样间隔。

进一步，所述的步骤(2)为定义与空间域网格模型中坐标为(i,j)的点相对应的拟空间域坐标为(τi,τj)；为了实现基于拟空间域声波方程的有限差分计算，需要对式(2)中各一阶导数项进行差分表示，得到和的拟空间域2n阶精度差分表达式。

进一步，的拟空间域2n阶精度差分表达式的获得方法为设应力波场p对变量τx具有2n阶导数，其中n为正整数，则应力波场p对变量τx在τx＝τi+1/2处一阶导数的2n阶精度差分表达式为：

式中或表示拟空间域τx方向上以τx＝τi+1/2为中心点的不同网格点的坐标，m＝1,2,…,n-1,n，其中τi+1/2表示τi+1与τi之间的拟空间域半网格点位置；表示各阶差分系数，m＝1,2,…,2n-1,2n，可通过求解式(4)所示的方程组求解；

式中和表示拟空间域传播时间，m＝1,2,…,n-1,n，其可分别表示为

其中表示中心点(τi+1/2,τj)与变量p所处网格点(τi+m,τj)之间的拟空间域传播时间，表示中心点(τi+1/2,τj)与变量p所处网格点(τi-(m-1),τj)之间的拟空间域传播时间；

按上述类似的方式分别得到p在τz＝τj+1/2处一阶导数的拟空间域2n阶精度差分系数和差分表达式，vx在τx＝τi处一阶导数的拟空间域2n阶精度差分系数和差分表达式，以及vz在τz＝τj处一阶导数的拟空间域2n阶精度差分系数和差分表达式，其中所有m均等于1,2,…,2n-1,或2n。

进一步，所述的步骤(3)为将p,vx,vz关于时间一阶导数的二阶精度差分格式以及步骤(2)得到的p,vx,vz关于τx和τz的一阶导数的拟空间域2n阶精度差分表达代入式(2)得到

其中k表示离散的时间点时刻，时间t＝k△t，△t表示离散时间步长；表示拟空间域中(τi,τj)点处在k时刻时应力波场p的值，同理可得以及的含义；式(7)即为拟空间域声波方程的2n阶精度有限差分表达式。

进一步，所述的步骤(4)为在边界区域，定义拟空间域衰减因子和将衰减因子引入拟空间域声波方程，得到拟空间域声波方程的完全匹配层边界条件：

其中px和pz表示应力p在τx和τz方向上的分量。和表示τx和τz方向上的衰减因子，其表达式为：

其中τm为拟空间域pml层中的点到中心波场外沿的垂直拟空间域距离，r为拟空间域pml层的理论反射系数，其取值为10^-5～10^-7，τl为拟空间域pml层厚度中最外层的点到中心波场外沿的垂直拟空间域传播时间，在需要进行边界吸收的区域令或

将拟空间域阻尼项和写成差分形式并与各拟空间域一阶导数项差分格式同时代入边界条件表达式，得到拟空间域声波方程完全匹配层边界条件2n阶精度有限差分格式，如下式所示：

进一步，所述的步骤(5)为给定观测系统参数、震源子波函数以及炮集记录，通过沿正时方向在每个k时刻将震源子波函数加载到式(7)中的应力波场p上，即可通过式(7)所示的拟空间域声波方程有限差分表达式递推求出全部时刻的p,vx,vz的正时延拓波场值；而通过沿逆时方向在每个k时刻将炮集记录加载到式(7)中的应力波场p上，即可通过式(7)所示的拟空间域声波方程有限差分表达式递推求出全部时刻的p,vx,vz的逆时延拓波场值；在边界区域，基于(10)式实现拟空间域声波方程人工边界反射的吸收衰减。

本发明与现有技术相比的有益效果：

本发明提出了拟空间域声波方程及其有限差分数值计算方法。在拟空间域，速度界面两侧的“传播时间”是按其实际的速度与传播距离分别计算，而与该网格线对应的时间采样间隔为不同“传播时间”段之和。理论上讲，在拟空间域不再存在速度界面的畸变问题，甚至还会减弱相邻网格点之间模型参量的突变程度，从而可在偏移计算中期望降低假散射和界面反射。通过模型试验计算可得，与常规声波方程相比，拟空间域声波方程有限差分方法可以很好地解决常规声波方程在逆时偏移剖面中速度界面形态畸变的问题；此外，基于拟空间域声波方程进行波场延拓可明显弱化界面假散射和层间反射波，从而进一步提高偏移成像的质量。

附图说明

图1为常规有限差分法网格剖分示意图；

图2为含倾斜界面的两层速度模型(a)原光滑界面模型(b)10m网格间距的模型；

图3为合成炮集记录示例(第76炮)；

图4为倾斜速度界面模型空间域网格剖分示意图；

图5为逆时偏移正时波场波前快照(第76炮，0.9s时刻)(a)基于常规声波方程(b)基于拟空间域声波方程；

图6为逆时偏移地震剖面(a)基于常规声波方程(b)基于拟空间域声波方程；

图7为逆时偏移地震剖面局部放大图示(a)基于常规声波方程(b)基于拟空间域声波方程。

具体实施方式

下面通过实施例结合附图来对本发明的技术方案作进一步解释，但本发明的保护范围不受实施例任何形式上的限制。

本发明提出了一种基于拟空间域声波方程的波场有限差分数值计算方法。其主要思想是将常规空间域变换到拟空间域后，空间域的等间隔的“距离”网格步长转化为可非等间隔的拟空间域“传播时间”步长，由此人们可在严格定义速度模型的基础上精确计算速度界面两侧的“传播时间”。在此基础上，给出了该方程及其完全匹配层(pml)边界条件的2n(n为正整数)阶精度有限差分表达式，然后基于此实现拟空间域声波波场的有限差分数值计算。

本发明的主要实施过程分为五个步骤：(1)根据单位网格长度将连续的地下介质模型按照矩形网格进行离散化，利用离散化矩形网格点上的介质速度求出各网格点处的拟空间域采样间隔，从而将空间域的声波方程转化为拟空间域声波方程；(2)在每个拟空间域差分计算网格点上，求取拟空间域声波方程中各一阶导数项和的拟空间域差分系数计算方程组中所需的一系列拟空间域传播时间，并由此计算各拟空间域差分计算网格点处的差分系数，进而得出各一阶导数项拟空间域2n阶精度差分表达式；(3)将各一阶导数差分表达式代入拟空间域声波方程中，得到整个拟空间域声波方程2n阶精度差分表达式；(4)在边界区域，定义拟空间域衰减因子和将衰减因子引入拟空间域声波方程，得到拟空间域声波方程完全匹配层边界条件，然后将各拟空间域阻尼项和写成差分形式并与各拟空间域一阶导数项差分格式同时代入边界条件表达式，得到拟空间域声波方程完全匹配层(pml)边界条件2n阶精度有限差分格式；(5)根据观测系统参数，给定震源子波函数、炮集记录以及波场初始条件，基于拟空间域声波方程2n阶精度差分表达式和拟空间域声波方程完全匹配层边界条件的2n阶精度有限差分格式，沿正时方向或逆时方向递推计算各时刻的正时延拓波场值或逆时延拓波场值。

实施例1

本实验的主要目的是验证拟空间域声波方程有限差分计算方法在逆时偏移中解决速度界面畸变以及压制界面假散射和层间反射波方面的有效性。

实验采用一个含有倾斜界面的两层速度模型，如图2(a)所示，其横向和纵向长度分别为4000m和2000m，界面上下两侧纵波速度分别为2500m/s和3500m/s，密度为2000kg/m³。将此倾斜界面模型以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可以逐一求出地震波在模型中各个长度为△ξ的网格上的拟空间域采样间隔为△τξ，τξ代表τx或τz，其中τx和τz分别表示拟空间域的水平方向和垂直方向；这里空间网格△ξ与拟空间域采样间隔△τξ之间满足关系△ξ＝cp△τξ。由此可根据式(1)将应力波场p和速度分量波场vx和vz对空间域的导数变换为对拟空间域的导数。根据式(1)可将常规空间域声波方程变换为拟空间域声波方程，其表达式如式(2)所示。

这样经过拟空间域坐标变换后，与空间域点(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)可以看出，波动方程中不再有速度参数，即将波动方程变换到拟空间域后，空间域的等间隔的“距离”网格步长转化为可非等间隔的“传播时间”步长，由此人们在严格定义速度模型的基础上来精确计算速度界面两侧的“传播时间”。例如，图4是包含一条倾斜界面的速度模型(图中黑色实线为速度界面)。可以看出，将速度界面模型按照常规矩形网格划分后，倾斜界面畸变为明显的阶梯状折线。但在拟空间域，速度界面两侧的“传播时间”是按其实际的速度与传播距离分别计算，而与该网格线对应的时间采样间隔为不同“传播时间”段之和。在图4中p和q点为速度模型按照矩形网格剖分之后两相邻空间网格点，黑色实线所示的速度界面与网格线pq相交于b点，此时p和q点之间的拟空间域采样间隔为△τ＝△τpb+△τbq，其中△τpb表示线段pb之间的旅行时，△τbq表示线段bq之间的旅行时。

理论上讲，在拟空间域不再存在速度界面的畸变问题，甚至还会减弱相邻网格点之间模型参量的突变程度，从而可在偏移计算中期望降低假散射和界面反射。

(2)定义与空间域网格模型中坐标为(i,j)的点相对应的拟空间域坐标为(τi,τj)。为了实现基于拟空间域声波方程的有限差分计算，需要对式(2)中各一阶导数项进行差分表示，得到和的拟空间域2n阶精度差分表达式。

以为例，设应力波场p对变量τx具有2n阶导数，其中n为正整数，则应力波场p对变量τx在τx＝τi+1/2处一阶导数的2n阶精度差分表达式为：

式中或表示拟空间域τx方向上以τx＝τi+1/2为中心点的不同网格点的坐标，m＝1,2,…,n-1,n，其中τi+1/2表示τi+1与τi之间的拟空间域半网格点位置；表示各阶差分系数，m＝1,2,…,2n-1,2n，可通过求解式(4)所示的方程组求解。

式中和表示拟空间域传播时间，m＝1,2,…,n-1,n，其可分别表示为

其中表示中心点(τi+1/2,τj)与变量p所处网格点(τi+m,τj)之间的拟空间域传播时间，表示中心点(τi+1/2,τj)与变量p所处网格点(τi-(m-1),τj)之间的拟空间域传播时间。由于拟空间域传播时间与网格速度和空间网格大小有关，因此即使在差分阶数相同的条件下，不同的速度和空间网格大小对应的差分系数也不同(表1和表2分别为空间网格点速度为3000m/s、空间网格步长为5m和空间网格点速度为4000m/s、空间网格步长为10m时的各阶精度差分系数)。

表1不同精度的拟空间域差分系数(速度v＝3000m/s，空间网格步长△ξ＝5m)

表2不同精度的拟空间域差分系数(速度v＝4000m/s，空间网格步长△ξ＝10m)

按上述类似的方式分别得到p在τz＝τj+1/2处一阶导数的拟空间域2n阶精度差分系数和差分表达式，vx在τx＝τi处一阶导数的拟空间域2n阶精度差分系数和差分表达式，以及vz在τz＝τj处一阶导数的拟空间域2n阶精度差分系数和差分表达式，其中所有m均等于1,2,…,2n-1,或2n。

(3)将p,vx,vz关于时间一阶导数的二阶精度差分格式以及步骤(2)得到的p,vx,vz关于τx和τz的一阶导数的拟空间域2n阶精度差分格式代入式(2)可得到

其中k表示离散的时间点时刻，时间t＝k△t，△t表示离散时间步长；表示拟空间域中(τi,τj)点处在k时刻时应力p的波场值，同理可得以及等的含义。式(7)即为拟空间域声波方程的2n阶精度有限差分表达式。

(4)在边界区域，定义拟空间域衰减因子和将衰减因子引入拟空间域声波方程，得到拟空间域声波方程的完全匹配层边界条件：

其中px和pz表示应力p在τx和τz方向上的分量。和表示τx和τz方向上的衰减因子，其表达式为：

其中τm为拟空间域pml层中的点到中心波场外沿的垂直拟空间域传播时间，r为拟空间域pml层的理论反射系数，其取值为10^-5～10^-7，τl为拟空间域pml层中最外层的点到中心波场外沿的垂直拟空间域传播时间，在需要进行边界吸收的区域令或

将拟空间域阻尼项和写成差分形式并与各拟空间域一阶导数项差分格式同时代入边界条件表达式，得到拟空间域声波方程完全匹配层边界条件的2n阶精度有限差分格式，如下式所示：

(5)实验中建立道固定、炮移动的观测系统，炮点位于500m至3480m之间，炮间隔为20m，共150炮；每炮接收道数为401道，各道位于0m至4000m之间，道间隔为10m；炮点与接收点深度均为10m，震源子波函数采用35hz的ricker子波，炮集记录如图3所示。通过沿正时方向在每个k时刻将35hz的ricker子波加载到式(7)中的应力波场p上，即可通过式(7)所示的拟空间域声波方程有限差分表达式递推求出全部时刻的p,vx,vz的正时延拓波场值；而通过沿逆时方向在每个k时刻将图3所示所示的炮集记录加载到式(7)中的应力波场p上，即可通过式(7)所示的拟空间域声波方程有限差分表达式递推求出全部时刻的p,vx,vz的逆时延拓波场值；在边界区域，基于(10)式实现拟空间域声波方程人工边界反射的吸收衰减。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)中箭头处界面反射波可以看出，基于拟空间域声波方程的有限差分计算方法对界面反射波(尤其是近垂直入射的反射波)具有明显压制效果。由此可见，本发明所述方法在压制界面假散射和层间反射波的有效性。