CN107479092B - 一种基于方向导数的频率域高阶声波方程正演模拟方法 - Google Patents

Abstract

本发明公开了一种基于方向导数的频率域高阶声波方程正演模拟方法，属于地震勘探技术领域，旨在提供一种模拟精度更高的频率域二维标量声波方程正演模拟方法。方法包括：根据频率域二维标量声波方程，利用方向导数建立其包含多个加权系数的四阶17点有限差分方程；进行归一化相速度频散分析，通过优化算法求取最优化加权系数；构建带有吸收边界条件的有限差分方程；利用该四阶17点有限差分方程进行地震波场数值模拟，得到地震波正演记录。本发明能够最大程度地压制频散，提高地震波场数值模拟的精度，还能够适应纵横向网格尺寸不等的情况。本发明主要用于地震勘探技术领域，为地震波场模拟与分析、地震反演成像以及地质建模等提供基础数据和技术支持。

Description

一种基于方向导数的频率域高阶声波方程正演模拟方法

技术领域

本发明属于地震勘探技术领域，涉及一种频率域高阶声波正演模拟方法。

背景技术

地震反演就是一种可以获得地下构造形态和物性分布的地球物理反演方法，分为基于褶积模型的反演和基于波动方程的反演。后者根据全部波场信息重建波动方程系数项(即物性参数)，能提供比较可靠的储层地质参数，在油藏描述、岩性勘探和开发地震研究中起着重要的作用，是地震勘探反问题的一个热点研究方向。

全波形反演是一种有效的波动方程反演方法，即利用叠前地震波场的运动学和动力学信息重建地层结构、获取物性参数，具有揭示复杂地质背景下的地质构造与储层物性细节信息的能力。全波形反演可以分为时间域反演、频率域反演、拉普拉斯域反演以及混合域反演。频率域全波形反演的存在以下优势：首先，频率域全波形反演只需使用几个频率的地震数据就能够实现高精度的建模，并且更容易引入吸收衰减因子等参数，所以在实际资料应用中更有意义；其次，由于频率域波场是相互解耦的，可以根据需要选用部分或全部频段的数据同时进行反演，更容易实现从低频到高频的多尺度反演策略，从而解决全波形反演的强非线性问题；此外，频率域波场正演，由于单频波场和震源的线性关系，在解决多源问题时存在天然的优势，比时间域有更高的计算效率，且不存在累积误差。无论是时间域全波形反演还是频率域全波形反演抑或是拉普拉斯域全波形反演，都离不开地震波场正演模拟。

频率域正演是频率域全波形反演的基础，正演方法的选择决定着反演的精度和效率。改进有限差分方案能够有效地提高频率域正演的计算精度和效率。Pratt andWorthington(1990)提出了一种经典的二维频率域声波方程5点有限差分方法，能够适用于不等的方向采样间隔；但是这种差分格式的数值频散严重，最小波长内需要有13个网格才能满足相速度误差低于1％。为了减少频率域正演模拟中最小波长需要的网格点数，旋转坐标系被引入到了差分格式。基于旋转坐标系，Jo等(1996)提出了优化9点有限差分格式，在相速度误差控制在1％内时，最小波长所需网格数减小到4个，且明显压制了数值频散，提高了计算的精度。基于优化9点差分格式，Shin和Sohn(1998)引入了0°、26.6°、45°、63.4°四种旋转坐标系来近似拉普拉斯算子，从而构造出了25点有限差分格式，该差分格式使得每个波长内只要2.5个网格点就能将相速度误差控制在1％内，但是阻抗矩阵的带宽却增加了一倍。刘璐等(2013)从减小系数矩阵带宽的角度出发得到了优化15点有限差分格式，在减小带宽的同时保持了相当的计算精度(每个波长内只需2.97个网格)，并提高了计算效率。由于以上这些差分格式都是二阶精度，不能满足高精度成像的需求。为此，曹书红和陈景波(2012)在常规四阶9点差分格式的基础上引入了45°方向的旋转坐标系，提出了高精度的四阶优化17点有限差分格式，将最小波长内所需网格数降低到了2.56个。

在实际使用中，当地质模型比较复杂的时候，通常会出现纵横采样间隔不相等的情况，这时上述提出的基于旋转坐标系差分方法不再适用。因此，针对上述旋转坐标系差分方法的局限性，Chen(2012)基于平均导数方法(Average-derivative method:ADM)提出了ADM二阶优化9点有限差分格式，将最小波长内所需网格数控制在了4个，并将其推广到了三维介质和粘滞声波方程。此外，Chen(2013)基于方向导数技术(Directional-derivativemethod:DDM)提出了一般的优化9点有限差分格式。这两种差分格式都能够适用于不同的方向空间采样间隔，增强了9点有限差分方法的灵活性并扩展了应用范围。申请号为201310522864.5的发明专利申请就公开了一种频率域优化混合交错网格有限差分正演模拟方法，其步骤为：①给出时间域二维声波方程；②消除人工边界反射，得带完全匹配层边界条件的时间域二维声波方程；③对方程左右两边时间变量进行傅立叶变换得频率域声波方程；④对匹配层边界条件频率域声波方程按常规交错网格进行有限差分离散，得有限差分离散格式；⑤对匹配层边界条件频率域声波方程按旋转交错网格进行有限差分离散得有限差分离散格式；⑥将常规交错网格和旋转交错网格优化混合，用两套网格系统中的加权平均，质量加速度项为中心点与其周围8点的加权平均；⑦在相速度误差最小的准则下，求取最优化系数。本发明加权系数使有限差分离散引起的频散误差最小，大大的提高了频率域正演模拟的精度。考虑到高精度成像的需要，张衡等(2014)将平均导数方法推广了四阶有限差分，提出了ADM四阶25点有限差分格式，并且最小波长内所需网格数为2.78。唐祥德等(2015)将平均导数法应用于常规的17点差分，提出的ADM四阶17点有限差分格式，使得最小波长内所需网格数降到了2.4，基于平均导数法的四阶有限差分格式，不仅提高了算法的灵活性，而且还进一步提高了计算精度。同样的，李美娟等(2016)将平均导数法15差分格式相结合，导出了一种能适应纵横向采样间隔不等情况的ADM优化15有限差分格式；而且最小波长内所需网格数降低到了2.85个，其在提高计算效率的同时兼顾了模拟精度。

综上所述，即使上述技术一步一步地对模拟精度进行了提高，但是目前的频率域标量声波方程正演技术仍不能满足高精度成像的需求，且部分现有的高阶正演模拟方法存在应用条件限制的问题，使得频率域高阶声波方程正演的模拟精度仍然较低。

发明内容

本发明的目的在于：提供一种基于方向导数的频率域高阶声波方程正演模拟方法，将方向导数技术与旋转坐标系相结合，克服常规四阶优化17点有限差分格式不能适用于纵横向网格尺寸不等的限制，提高数值模拟的精度。

本发明采用的技术方案如下：

一种基于方向导数的频率域高阶声波方程正演模拟方法，综合利用方向导数技术和旋转坐标系得到能够适应纵横向网格尺寸不等情况的高阶有限差分格式，最大程度地压制频散，提高地震波场数值模拟的精度，包括以下步骤：

步骤1：根据频率域二维标量声波方程利用方向导数技术建立包含多个加权系数的四阶17点有限差分方程：

其中，P_m,n＝P(m△x,n△z)表示离散网格点(m,n)处的压力波场，△x、△z分别表示速度模型在X轴方向、Z轴方向的采样间隔，下标m、n分别表示X轴方向、Z轴方向的网格坐标，v_m,n表示速度模型离散网格点(m,n)处的速度，ω_j为计算角频率，下标j为角频率离散点号，a、b、c、d、e、f均为加权系数，且b+4c+4d+4e+4f＝1，差分方程左边的第一项为原始正交坐标系下拉普拉斯算子的四阶差分项、第二项为旋转坐标系下用方向导数获得拉普拉斯算子的四阶差分项、第三项为质量加速度项；

步骤2：进行归一化相速度频散分析，通过优化算法求取最优化加权系数；

步骤3：构建带有吸收边界条件的有限差分方程；

步骤4：利用四阶17点有限差分方程进行地震波场数值模拟，得到地震波正演记录。

作为优选，在步骤1中，利用方向导数技术建立包含多个加权系数的四阶17点有限差分方程的步骤为：

步骤1.1：将频率域二维标量声波方程的拉普拉斯算子项进行泛化，表示为原始正交坐标系下的拉普拉斯算子和旋转坐标系下的用方向导数获得的拉普拉斯算子的加权组合：

其中，分别为原始正交坐标系下的拉普拉斯算子和旋转坐标系下的用方向导数获得的拉普拉斯算子；

步骤1.2：将四阶有限差分格式分别应用于正交坐标系下的拉普拉斯算子和旋转坐标系下的用方向导数获得的拉普拉斯算子，得到拉普拉斯算子项的四阶有限差分格式：

步骤1.3：将频率域二维标量声波方程的质量加速度项表示为对应于拉普拉斯算子差分网格点的线性组合：

步骤1.4：将步骤1.2得到的拉普拉斯算子的四阶有限差分格式、步骤1.3得到的质量加速度的线性组合公式带入频率域二维标量声波方程得到包含多个加权系数的四阶17点有限差分方程。

作为优选，步骤1.1中：旋转坐标系下用方向导数获得拉普拉斯算子的方法为：利用方向导数技术，推导出拉普拉斯算子：

其中，△r＝(△x²+△z²)^1/2，l₁,l₂为旋转坐标系的两个方向，当△x≠△z时，方向l₁与l₂不再是相互正交的。

作为优选，在步骤2中，求取在频散最小条件下的最优化加权系数的步骤为：

步骤2.1：基于公式(1)，导出归一化相速度：

将平面波带入公式(1)，化简得到如下离散的归一化相速度公式，

其中，

A＝cos(△xk_x)，B＝cos(2△xk_x)，C＝cos(△zk_z)，D＝cos(2△zk_z)，V_ph＝ω/k为相速度，

P₀为平面波振幅，i为虚数单位，k_x、k_z分别为X轴方向、Z轴方向的波数，a、b、c、d、e、f均为加权系数，k_x＝k sinθ，k_z＝k cosθ，θ为平面波传播方向于Z轴正方向的夹角，为平面波波数，△＝max(△x,△z)，即取最大的采样间隔，G表示每个波长内需要的网格数；当△x≠△z时，波数k的取值由较大的采样间隔决定，并且在求取最优化加权系数的时候需要分为△x≥△z和△x<△z两种情况分别进行讨论；由于公式(1)中的拉普拉斯算子项和质量加速度项的有限差分格式在△x≥△z和△x<△z两种情况下具有几何对称性，即在求取最优化加权系数的时候只需要考虑其中一种情况；

当△x≥△z时，令网格比r＝△x/△z，则化简上述等式(2)得到△x≥△z情况下的归一化相速度公式：

其中，

步骤2.2：建立频散分析目标函数，即频散误差最小化目标函数：

其中将公式(3)带入公式(4)中得到具体的频散误差最小化目标函数；

步骤2.3：采用模式搜索法求解目标函数，得到目标函数最小条件下不同网格比r对应的最优加权系数a、b、c、d、e、f，具体参见图8。

作为优选，在步骤3中，构建带有吸收边界条件的有限差分方程，通常采用完全匹配层(Perfectly Matched Layer，PML)边界，添加了吸收边界条件的有限差分方程为：

其中，

e_x,e_z分别为X轴和Z轴方向的拉伸函数，d_x,d_z分别为X轴和Z轴方向的衰减因子。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

本发明中，采用方向导数技术和旋转坐标系组合，扩展了常规了四阶优化17点有限差分格式的适用范围，即通过最优化算法获得不同纵横向网格比对应的最优化加权系数，使其能够适应不等纵横向网格尺寸的情况；通过模式搜索法得到的最优化的加权系数，在最大程度上压制了数值频散；再加上四阶有限差分，有效地提高了数值模拟的精度，为高精度成像奠定了坚实的基础；最小波场内所需网格数降到了2.4，在相同精度要求下有效地降低了内存需求和计算量，优于其他类似方法；采用最优化方法求取加权系数，极大地压制了数值频散，有效地提高计算精度；从而解决了目前的频率域标量声波方程正演技术不能满足高精度成像的需求，以及部分现有的高阶正演方法存在应用条件限制的问题

附图说明

图1为本发明的方法流程图；

图2为17点有限差分格式网格分布图；

图3为简单的层状介质速度模型及观测系统分布图；

图4为正演结果剖面图；

图5为二维地堑速度模型及观测系统分布图；

图6为正演结果剖面图；

图7为采用时间域正演方法得到的正演结果剖面图；

图8为△x≥△z时对应不同网格比的最优化加权系数。

具体实施方式

为了使本发明的目的、方法流程和特征更加清楚，下面将结合附图对本发明实施例和应用进行详细说明。具体实施例和应用中所提供的描述信息仅为示例，用于解释本发明，并能作为本发明的限定。这里所描述的实施例的各种延伸和组合对于本领域的技术人员是显而易见的，在不脱离本发明的实质和范围的情况下，本发明定义的一般原则可以应用到其他实施例和应用中。

下面将结合附图对本发明进行详细说明。

图1是本发明方法的流程图：

步骤1：根据频率域二维标量声波方程利用方向导数技术建立包含多个加权系数的四阶17点有限差分方程：

其中，P_m,n＝P(m△x,n△z)表示离散网格点(m,n)处的压力波场，△x、△z分别表示速度模型在X轴方向、Z轴方向的采样间隔，下标m、n分别表示X轴方向、Z轴方向的网格坐标，v_m,n表示速度模型离散网格点(m,n)处的速度，ω_j为计算角频率，下标j为角频率离散点号，a、b、c、d、e、f均为加权系数，且b+4c+4d+4e+4f＝1，差分方程左边的第一项为原始正交坐标系下拉普拉斯算子的四阶差分项、第二项为旋转坐标系下用方向导数获得拉普拉斯算子的四阶差分项、第三项为质量加速度项；具体的17点有限差分网格模板参见图2，图中黑色圆点代表差分格式中的网格结点。

上述步骤1中的基于方向导数的四阶17点有限差分方程，能够适用于△x≠△z的情况，也可以针对于△x＝△z的情况，该四阶17点有限差分方程的建立包括以下步骤：

步骤1.1：将频率域二维标量声波方程的拉普拉斯算子项进行泛化，表示为原始正交坐标系下的拉普拉斯算子和旋转坐标系下用方向导数获得拉普拉斯算子的加权组合：

其中，分别为原始正交坐标系下的拉普拉斯算子和旋转坐标系下用方向导数获得拉普拉斯算子；

在步骤1.1中：旋转坐标系下用方向导数获得拉普拉斯算子的方法为：利用方向导数技术，推导出拉普拉斯算子：

其中，△r＝(△x²+△z²)^1/2，l₁,l₂为旋转坐标系的两个方向，当△x≠△z时，方向l₁与l₂不再是相互正交的。

步骤1.2：将四阶有限差分格式分别应用于正交坐标系下的拉普拉斯算子和旋转坐标系下用方向导数获得拉普拉斯算子，得到拉普拉斯算子项的四阶有限差分格式：

步骤1.3：将频率域二维标量声波方程的质量加速度项表示为对应于拉普拉斯算子差分网格点的线性组合：

步骤1.4：将步骤1.2得到的拉普拉斯算子项的四阶有限差分格式、步骤1.3得到的质量加速度的线性组合公式带入频率域二维标量声波方程得到包含多个加权系数的四阶17点有限差分方程。

步骤2：进行归一化相速度频散分析，通过优化算法求取在频散最小条件下的最优化加权系数，具体包括以下步骤：

步骤2.1：基于公式(1)，导出归一化相速度：

将平面波带入公式(1)，化简得到如下离散的归一化相速度公式，

其中，

A＝cos(△xk_x)，B＝cos(2△xk_x)，C＝cos(△zk_z)，D＝cos(2△zk_z)，V_ph＝ω/k为相速度，

P₀为平面波振幅，i为虚数单位，k_x、k_z分别为X轴方向、Z轴方向的波数，a、b、c、d、e、f均为加权系数，k_x＝k sinθ，k_z＝k cosθ，θ为平面波传播方向于Z轴正方向的夹角，为平面波波数，△＝max(△x,△z)，即取最大的采样间隔，G表示每个波长内需要的网格数；当△x≠△z时，波数k的取值由较大的采样间隔决定，并且在求取最优化加权系数的时候需要分为△x≥△z和△x<△z两种情况分别进行讨论；由于公式(1)中的拉普拉斯算子项和质量加速度项的有限差分格式在△x≥△z和△x<△z两种情况下具有几何对称性，即在求取最优化加权系数的时候只需要考虑其中一种情况；

当△x≥△z时，令网格比r＝△x/△z，则化简上述等式(2)得到△x≥△z情况下的归一化相速度公式：

其中，

步骤2.2：建立频散分析目标函数，即频散误差最小化目标函数：

其中将公式(3)带入公式(4)中得到具体的频散误差最小化目标函数；

步骤2.3：采用模式搜索法求解目标函数，得到目标函数最小条件下不同网格比r对应的最优加权系数a、b、c、d、e、f，具体参见图8。

步骤3：采用完全匹配层(Perfectly Matched Layer，PML)边界，构建带有吸收边界条件的有限差分方程，添加了吸收边界条件的有限差分方程为：

其中，

e_x,e_z分别为X轴和Z轴方向的拉伸函数，d_x,d_z分别为X轴和Z轴方向的衰减因子。

步骤4：利用四阶17点有限差分方程进行地震波场数值模拟包括：输入速度模型、加载震源子波、设置正演参数、进行频率循环计算所有单频波场、提取检波器接收波场并输出正演记录。

上述步骤4中的内容均为现有技术中频率于正演的常规方法流程，故在本发明申请中不再赘述。

为了更清楚地展示上述方法及装置的正演效果，包括同相轴位置、边界吸收以及频散压制情况，下面结合具体实施例对本发明进行进一步说明。

实施例1

采用上述的模拟方法对一个简单的两层介质速度模型进行正演模拟试验。图3为速度模型结构和观测系统分布，速度模型速度分别为1000/s和2000/s，图中红色倒三角形代表震源，黑色的正三角形代表检波器。速度模型网格间距△x＝△z＝5m，纵横向网格点数为101×101。观测系统类型为中间激发两端接收，设置一个单震源观测系统，其震源位置为(x_s,z_s)＝(250m,10m)，排列中检波器对称分布于震源两侧，道间距为5m，检波器埋深10m，总共11道接收。震源采用主频为20Hz的雷克子波，时间采样间隔为1ms，记录时间长度为1s。由于该速度模型的纵横向网格尺寸相等，所以选择图8中r＝1时的加权系数构建稀疏阻抗矩阵。图4为其正演结果剖面。正演结果中的同相轴符合地震波传播规律且位置正确，边界反射吸收效果好，数值频散压制效果好。试验表明，由于采用了四阶精度的差分格式，本发明方法具有较高的数值模拟精度，且最优化加权系数保证了较好的频散压制。

实施例2

采用上述的模拟方法对一个二维的地堑速度模型进行正演模拟试验。图5为速度模型结构与观测系统分布，速度模型速度从上到下分别为1500m/s，2000m/s和3000m/s，图中红色倒三角形代表震源，黑色的正三角形代表检波器。速度模型网格间距△x＝10m，△z＝5m，纵横向网格点数为201×201。与上述实施例1采用相同类型的观测系统和震源参数，震源位置为(x_s,z_s)＝(1000m,5m)，道间距为10m，检波器埋深5m，总共201道接收。时间采样间隔1ms，记录时间长度为1.2s。由于该速度模型的纵横向网格尺寸不相等，△x/△z＝2，所以选择图8中的r＝2时的加权系数构建稀疏阻抗矩阵。图6为本发明方法实施例的正演结果，图7为时间域方法的正演结果。通过对比发现两个正演结果一致，说明本发明方法是正确的，能适用于复杂的地质模型，对纵横向网格尺寸不等的情况具有较强的鲁棒性。

上述即为本发明的实施列，本发明不限于上述实施方式。本发明可以以任何有效的形式实现，包括硬件、软件、固件和它们的任意组合。本发明实施列的元素和模块可以在物理上、功能上以及逻辑上以任何有效的方式实现。任何相关领域的技术人员，在本发明的启发下做出的修改和组合，凡是采用了相同的机制和方法，均落入本发明的保护范围之内。

Claims (6)

1.一种基于方向导数的频率域高阶声波方程正演模拟方法，其特征在于，包括以下步骤：

步骤1：根据频率域二维标量声波方程利用方向导数技术建立包含多个加权系数的四阶17点有限差分方程：

其中，P_m,n＝P(m△x,n△z)表示离散网格点(m,n)处的压力波场，△x、△z分别表示速度模型在X轴方向、Z轴方向的采样间隔，下标m、n分别表示X轴方向、Z轴方向的网格坐标，v_m,n表示速度模型离散网格点(m,n)处的速度，ω_j为计算角频率，下标j为角频率离散点号，a、b、c、d、e、f均为加权系数，且b+4c+4d+4e+4f＝1，差分方程左边的第一项为原始正交坐标系下拉普拉斯算子的四阶差分项、第二项为旋转坐标系下用方向导数获得拉普拉斯算子的四阶差分项、第三项为质量加速度项；

步骤2：进行归一化相速度频散分析，通过优化算法求取最优化加权系数；

步骤3：构建带有吸收边界条件的有限差分方程；

步骤4：利用四阶17点有限差分方程进行地震波场数值模拟，得到地震波正演记录。

2.如权利要求1所述的一种基于方向导数的频率域高阶声波方程正演模拟方法，其特征在于，步骤1中，利用方向导数技术建立包含多个加权系数的四阶17点有限差分方程的步骤包括：

步骤1.1：将频率域二维标量声波方程的拉普拉斯算子项进行泛化，表示为原始正交坐标系下的拉普拉斯算子和旋转坐标系下的用方向导数获得的拉普拉斯算子的加权组合；

步骤1.2：将四阶有限差分格式分别应用于正交坐标系下的拉普拉斯算子和旋转坐标系下的用方向导数获得的拉普拉斯算子，得到拉普拉斯算子项的四阶有限差分格式；

步骤1.3：将频率域二维标量声波方程的质量加速度项表示为对应于拉普拉斯算子差分网格点的线性组合；

步骤1.4：将步骤1.2得到的拉普拉斯算子项的四阶有限差分格式、步骤1.3得到的质量加速度的线性组合公式带入频率域二维标量声波方程，得到包含多个加权系数的四阶17点有限差分方程。

3.如权利要求2所述的一种基于方向导数的频率域高阶声波方程正演模拟方法，其特征在于，步骤1.1中：旋转坐标系下用方向导数获得拉普拉斯算子的方法为：利用方向导数技术，推导出拉普拉斯算子。

4.如权利要求1所述的一种基于方向导数的频率域高阶声波方程正演模拟方法，其特征在于，步骤2中，求取最优化加权系数的步骤为：

步骤2.1：基于公式(1)，导出归一化相速度：

将平面波带入公式(1)，化简得到如下离散的归一化相速度公式，

其中，

A＝cos(△xk_x)，B＝cos(2△xk_x)，C＝cos(△zk_z)，D＝cos(2△zk_z)，V_ph＝ω/k为相速度，P₀为平面波振幅，i为虚数单位，k_x、k_z分别为X轴方向、Z轴方向的波数，a、b、c、d、e、f均为加权系数，k_x＝ksinθ，k_z＝kcosθ，θ为平面波传播方向于Z轴正方向的夹角，为平面波波数，△＝max(△x,△z)，即取最大的采样间隔，G表示每个波长内需要的网格数；当△x≠△z时，波数k的取值由较大的采样间隔决定，并且在求取最优化加权系数的时候需要分为△x≥△z和△x<△z两种情况分别进行讨论；由于公式(1)中的拉普拉斯算子项和质量加速度项的有限差分格式在△x≥△z和△x<△z两种情况下具有几何对称性，即在求取最优化加权系数的时候只需要考虑其中一种情况；

当△x≥△z时，令网格比r＝△x/△z，则化简上述等式(2)得到△x≥△z情况下的归一化相速度公式：

其中，

步骤2.2：建立频散分析目标函数，即频散误差最小化目标函数：

其中将公式(3)带入公式(4)中得到具体的频散误差最小化目标函数；

步骤2.3：采用优化算法求解目标函数，得到目标函数最小条件下不同网格比r对应的最优加权系数a、b、c、d、e、f；优化算法包括：模式搜索算、全局搜索、遗传算法、模拟退火优化算法。

5.如权利要求1所述的一种基于方向导数的频率域高阶声波方程正演模拟方法，其特征在于，步骤3中，构建带有吸收边界条件的有限差分方程，用于吸收正演过程中来自边界的反射波场，添加了吸收边界条件的有限差分方程为：

其中，

e_x,e_z分别为X轴和Z轴方向的拉伸函数，d_x,d_z分别为X轴和Z轴方向的衰减因子。

6.如权利要求1所述的一种基于方向导数的频率域高阶声波方程正演模拟方法，其特征在于，步骤4中，利用四阶17点有限差分方程进行地震波场数值模拟包括：输入速度模型、加载震源子波、设置正演参数、进行频率循环计算所有单频波场、提取检波器接收波场并输出正演记录。