基于非分裂完全匹配层吸收边界的正演模拟方法与流程

1.本发明涉及地球物理勘查技术领域，特别是涉及到一种基于非分裂卷积完全匹配层弹性波数值模拟方法。


背景技术：

2.波动方程法是应对日益复杂的地质构造地震波传播数值模拟的最重要方法。由于无法构建无限区域的模型，使用波动方程法进行地震波正演数值模拟时，一般需要用有限的区域模拟无限区域波传播过程，即添加人工边界。为了消除地震波传播到人工边界引起的虚假反射，就需要设置吸收边界条件，使得地震波可以无反射的穿过人工边界。因此，发展形成有效而稳定的吸收边界条件一直是地震波数值模拟的一项重要研究内容。完全匹配层吸收边界是在模型区域外人为增加了一层吸收介质，当波传播到该介质中时，波的能量将被介质以指数衰减的速度迅速吸收，直至消失，而不会反弹回模型区域。理论上，通过合理的设置该吸收层的吸收系数，无论什么频率的波以任意角度从模型区域进入吸收介质层时，都不产生反射。与传统的吸收边界条件相比，完全匹配层具有更优秀的吸收效果。目前应用较为广泛的完全匹配层吸收边界为分裂完全匹配层(spml)吸收边界。但在大角度入射的情况下，spml吸收边界的吸收效果并不理想，虚假反射将会更加明显，对薄层介质、震源靠近模型离散的网格边缘或者大偏移距模型的数值模拟效果不理想，此外，spml边界无法吸收掠射波，并且正演时对内存需求较大。在分裂完全匹配层(spml)吸收边界的基础上发展起来的卷积完全匹配层(convolution perfectly matched layer，cpml)吸收边界能够吸收大角度入射波和掠射波，但目前构建的卷积完全匹配层吸收边界多适用于一阶波动方程，难以直接应用于二阶波动方程。
3.为解决以上技术问题，本发明提出一种基于非分裂完全匹配层吸收边界正演模拟方法。


技术实现要素：

4.本发明的目的是提供一种适用于二阶弹性波方程的非分裂卷积完全匹配层吸收边界条件的数值模拟方法。
5.本发明的目的可通过如下技术措施来实现：基于非分裂完全匹配层吸收边界的正演模拟方法，该基于非分裂完全匹配层吸收边界的正演模拟方法包括：步骤1，建立二阶弹性波动方程；步骤2，利用迭代法代替卷积运算的思想构建适用于二阶弹性波方程的非分裂完全匹配层吸收边界；步骤3，实现非分裂完全匹配层吸收边界二阶精度弹性波正演模拟。
6.本发明的目的还可通过如下技术措施来实现：
7.在步骤1中，通过柯西方程、几何方程以及纳维尔方程推导出二阶弹性波波动方程：
[0008][0009]
其中λ和μ为介质的拉梅常数；ρ表示介质密度；u表示地震波场；u
x
和uz分别为水平方向和垂直方向的位移波场，t表示地震波传播时间，x，z分别表示水平方向和垂直方向。
[0010]
在步骤2中，采用复坐标变换，忽略衰减因子空变影响计算出二阶弹性波方程的空间二阶偏导算子：
[0011][0012]
其中：δ(t)和h(t)分别代表狄拉克函数以及单位阶跃函数，t为时间，d
x
、dz为衰减透射波的衰减参数；α
x
、αz为滤波的频率参数；为x方向在复频域的拉伸变量，x与z方向在复频域的拉伸变量，k
x
，kz为衰减渐消波的尺度因子。
[0013]
在步骤2中，利用迭代法代替卷积运算的思想，将二阶偏导算子中的时域卷积运算转换成迭代运算计算出时间域基于非分裂完全匹配层吸收边界的二阶弹性波方程。
[0014]
在步骤3中，由连续方程的时间二阶、空间十阶有限差分计算基于非分裂完全匹配层吸收边界二阶弹性波方程的离散方程。
[0015]
在步骤3中，对均匀介质模型和复杂层状模型分别开展分裂完全匹配层和非分裂完全匹配层吸收边界二阶精度弹性波方程正演模拟，通过成图验证非分裂完全匹配层吸收边界具有更好的吸收效果、更快的计算速度和对计算机内存较低的需求。
[0016]
本发明的基于非分裂完全匹配层吸收边界的正演模拟方法，由连续方程的时间二阶、空间十阶有限差分计算基于非分裂完全匹配层吸收边界二阶弹性波方程的离散方程，基于卷积完全匹配层(cpml)吸收边界利用迭代运算取代卷积运算的思想推导出非分裂完全匹配层吸收边界(ncpml，new convolution perfectly matched layer，又称新卷积完全匹配层)，非分裂完全匹配层吸收边界能够吸收大角度的入射波和掠射波，能够直接对二阶弹性波方程做数值模拟，减小内存使用空间，提高计算效率，更适用于实际生产需求。
附图说明
[0017]
图1是本发明的基于非分裂完全匹配层吸收边界的正演模拟方法的一具体实施例流程图；
[0018]
图2是本发明的具体实施例中均匀介质模型纵波源激发分裂完全匹配层(spml)和非分裂完全匹配层(ncpml)两种不同吸收边界条件下二阶弹性波正演的水平分量波场图；(a)为spml边界700ms时刻的波场快照；(b)为ncpml边界700ms时刻的波场快照；
[0019]
图3是本发明的具体实施例中不同吸收边界条件下x＝900m处700ms时刻水平分量波形对比图；
[0020]
图4是本发明的具体实施例中层状速度模型的示意图；
[0021]
图5是本发明的具体实施例中层状模型t＝350ms时两种不同吸收边界条件下水平分量波场快照对比图；(a)层状模型spml吸收边界条件下水平分量波场快照；(b)层状模型
ncpml吸收边界条件下水平分量波场快照；
[0022]
图6是本发明的具体实施例中第150道spml吸收边界和ncpml吸收边界条件下的水平分量波形曲线与理论水平分量波形曲线对比图；
[0023]
图7是本发明的具体实施例中图6在深度为400m至800m处局部放大水平分量波形曲线对比图；
[0024]
图8是本发明的具体实施例中图6spml吸收边界和ncpml吸收边界波场水平分量的数值解与理论值的误差曲线对比图；
[0025]
图9是本发明的具体实施例中ncpml吸收边界的层状模型对应的地震记录的水平分量和垂直分量图,(a)水平分量地震记录，(b)垂直分量。
具体实施方式
[0026]
为使本发明的上述和其他目的、特征和优点能更明显易懂，下文特举出较佳实施例，并配合附图所示，作详细说明如下。
[0027]
图1所示为本发明的基于非分裂完全匹配层弹性波数值模拟方法的流程图。
[0028]
第一步：建立二阶弹性波波动方程
[0029]
通过柯西方程、几何方程以及纳维尔方程推导出二阶弹性波波动方程：
[0030][0031]
其中λ和μ为介质的拉梅常数；ρ表示介质密度；u
x
和uz分别为水平方向和垂直方向的位移波场；
[0032]
第二步：利用迭代法代替卷积运算的思想构建适用于二阶弹性波方程的非分裂完全匹配层吸收边界
[0033]
cpml吸收边界x方向的复拉伸变量可以表示为：
[0034][0035]
其中，
x
≥1,d
x
≥0,α
x
≥0，虚数部分d
x
＞0表示吸收层中衰减因子，用来衰减透射波，α
x
，d
x
是新引入的衰减因子，s
x
表示x方向复拉伸变量；在非吸收层中，k
x
＝1且d
x
＝0。
[0036]
采用以上复坐标变换公式(2)，将二阶弹性波方程表示为：
[0037][0038]
其中和分别表示纵、横波速度；表示波场水平分量，为波场垂直分量；s
x
表示x方向复拉伸变量，sz为z方向复拉伸变量；f-1
为fourier反变换。
[0039]
忽略衰减因子空变影响计算出方程(3)中空间二阶偏导算子：
[0040][0041]
其中：δ(t)和h(t)分别代表狄拉克函数以及单位阶跃函数，t为时间，d
x
、fz为衰减透射波的衰减参数；α
x
、αz为滤波的频率参数；为x方向在复频域的拉伸变量，为x与z方向在复频域的拉伸变量，k
x
，kz为衰减渐消波的尺度因子。
[0042]
利用迭代法代替卷积运算的思想，将时域卷积运算转换成迭代运算并计算得到时间域基于ncpml的二阶弹性波方程：
[0043]
1)当d
x
/k
x
α
x
≠dz/kz αz时，基于ncpml时间域的控制方程为：
[0044][0045]
2)当d
x
/k
x
α
x
＝dz/kz αz时，基于ncpml时间域的控制方程为：
[0046][0047]
方程(5)(6)中：是将卷积运算转换成迭代运算的迭代项；b
x
、a
x
为转换过程中产生的两个系数，为转换过程中产生的两个系数，类似。w为波场垂直分量，u为波场水平分量，t为采样时间。v
p
表示纵波速度、vs表示横波速度。其余符号含义与公式(4)相同。
[0048]
第三步：实现非分裂完全匹配层吸收边界二阶精度弹性波正演模拟为实现正演模拟，对连续方程进行时间二阶、空间十阶差分计算，差分方程格式如下：
[0049]
时间二阶精度差分格式：
[0050][0051]
空间十阶精度差分格式：
[0052]
[0053][0054]
其中a0、zn、均为差分系数，δx、δz为网格间距，u(x,z)为波场，t表示采样时间，δt表示时间采样间隔。
[0055]
离散后的弹性波方程为：
[0056]
1)当d
x
/k
x
α
x
≠dz/kz αz时，基于ncpml时间域的控制方程的离散形式为:
[0057][0058][0059]
2)当d
x
/k
x
α
x
＝dz/kz αz时，基于ncpml时间域的控制方程离散形式为：
[0060][0061]
方程(10)和方程(11)中，为n 1时刻的水平分量，为n 1时刻垂直分量，δt表示时间采样间隔，其它参数同公式(5)和公式(6)。
[0062]
正演模拟使用的子波主频为20hz，时间采样间隔为1ms，采样点数为1000。均匀模型震源位于网格点(100，100)，层状模型震源位于(1000m，0m)处。
[0063]
通过数值模拟对比ncpml边界和spml边界在吸收效果、运行速度以及内存占用方面的不同，并用复杂的层状模型测试ncpml边界对复杂模型的适定性。首先通过层状模型对比分析ncpml吸收边界与spml吸收边界的吸收效果以及计算效率；然后通过层状模型验证所推导的二阶弹性波ncpml吸收边界的稳定性。
[0064]
图2是本发明的一具体实施例中均匀介质模型纵波源激发的水平分量波场图，模型大小为200*200，vp＝3000m/s，vs＝2000m/s，网格间距为10m。(a)为spml边界700ms时刻的波场快照；(b)为ncpml边界700ms时刻的波场快照；对比图2(a)和图2(b)，直观上两种吸收边界对人工边界产生的反射波都具有良好的吸收效果，在边界处没有明显的虚假反射。定量分析两种吸收边界的吸收效果，图3是图2中虚线处(在x＝900m的垂线上)用两种pml吸收边界数值模拟得出的波场值与理论值的对比图，其中理论值是通过增加计算区域得到的，其计算结果等同于对边界处的透射波完全吸收。通过图2和图3可以看出与spml吸收边界相比，ncpml吸收边界更接近理论值，吸收效果更好。
[0065]
图4为存在多个速度界面的复杂层状速度模型。在层状模型数值模拟中，图5(a)为层状模型spml吸收边界t＝350ms时水平分量波场快照，(b)为层状模型ncpml吸收边界t＝350ms时水平分量波场快照；图6是第150道spml吸收边界和ncpml吸收边界的水平分量波形曲线与理论水平分量波形曲线对比；由图5和图6可以看出ncpml吸收边界与spml吸收边界的吸收效果都较好。为进一步分析两者的区别，对图6所示曲线进行局部放大显示，并分别计算两者与理论值的误差。图7是图6在深度为400m至800m处水平分量波形曲线；图8是用图6中spml吸收边界和ncpml吸收边界水平分量波场的数值解与理论值的误差曲线；通过图7和图8可以看出ncpml吸收边界的数值更加接近理论值，其吸收效果更好。
[0066]
表1 spml边界与ncpml边界计算效率与占用内存数据对比
[0067][0068]
表1为均匀介质模型spml边界与ncpml边界计算效率与占用内存的数据对比，通过表1的数据对比可知ncpml吸收边界相较于spml吸收边界计算速度更快，所占内存更小，更适用于生产实践。图9是基于ncpml吸收边界的层状模型对应的地震记录，图9(a)水平分量，图9(b)垂直分量；由图9可以看出，得到的地震记录都没有明显的虚假反射波，证明了本方法对于复杂模型的数值模拟稳定性良好。