CN115600373A - 黏滞各向异性介质qP波模拟方法、系统、设备及应用 - Google Patents

Abstract

本发明属于地震波场数值模拟数据识别技术领域，公开了黏滞各向异性介质qP波模拟方法、系统、设备及应用。所述方法包括：通过傅里叶变换将时间‑空间域VTI介质波动方程转换到频率‑空间域，推导出黏滞性VTI介质中qP波波动方程，消除横波干扰；采用有限差分方法将地震模型离散化，剖分为多个的网格，用网格点替代速度模型；通过构建频率‑空间域25点的波动方程的优化差分格式进行数值模拟计算。本发明具有很高的数值模拟精度。通过准纵波频率‑空间域波动方程25点优化差分算法的数值模拟可以研究地震波的传播特征，可以指导波场的识别和纵波波场的分离，同时也可以为反演、偏移成像等后续工作提供理论基础。

Description

黏滞各向异性介质qP波模拟方法、系统、设备及应用

技术领域

本发明属于地震波场数值模拟数据识别技术领域，尤其涉及一种黏滞各向异性介质qP波模拟方法、系统、设备及应用。

背景技术

目前，随着我国经济的快速持续发展，石油、天油气消费的对外依存度已超过50％，严重影响我国经济社会发展的安全。而且，随着我国国民经济的进一步持续发展，对外依存度还将进一步提高。同时，我国油气平均探明率只有29％，远低于国际64％的平均水平。因此加大我国油气资源的勘探开发力度成为当务之急和战略方向。但伴随着我国油气的勘探开发，剩余油气资源的勘探开发难度越来越大。目前我国的油气勘探开发中存在的主要问题有：

1)大规模油气资源已经非常匮乏，目前的勘探对象主要是低渗透、隐蔽型、裂缝型和低幅度构造油气藏，对这类油气藏进行地震勘探，在技术理论与处理技术等方面，均面临新的困难和问题；

2)复杂地质条件与陆上复杂地表条件下形成的低信噪比地震资料难以准确成像，对这些类型储层的精细勘探提出了巨大的挑战；

3)这些类型储层的复杂性使得当前的处理技术在地震数据成像和油气识别等问题上存在较大的误差，难以满足油气勘探开发的需要。

目前，中国的石油、天燃气地球物理勘探绝大多数都在构造地形复杂地区进行。复杂地质条件与陆上复杂地表条件下低信噪比的地震资料难以准确成像是其中一个非常重要的问题，地震勘探的记录数据信噪比低、面貌复杂、有效信息难以分辨，因此需要以数值模拟来辅助波场中各种波的识别及资料的处理和解释。地震数值模拟是通过模拟地震波在地下实际介质中的传播，研究地震波的传播规律与波场特征，是开展地震资料采集、处理和解释等研究工作的有效辅助方法，对正确认识地震波的传播规律、进行实际地震资料解释和处理有重要的意义。

地震波场数值模拟根据实现域的不同主要分为时间-空间域和频率-空间域两种。时间-空间域数值模拟是通过对时间域波动方程进行离散，然后按照时间步进的方法模拟地震波的传播。频率-空间域数值模拟是通过对频率域波动方程进行离散并求解各个频率的波场，将不同频率的波场加在一起进行傅里叶逆变换得到模拟的地震波场。与时间-空间域数值模拟方法相比，在频率域求解波动方程的优点在于：1、在频率域求解波动方程不会累加误差。这主要是因为在频率域进行求解时，不同频率之间互不干扰，计算模型中全部网格节点的不同频率切片，将误差分配到所有的网格点，因此地震波场不会累加误差。2、在频率-空间域进行数值模拟各频率切片之间的计算可以实现并行运算，大大提高了计算效率。3、将波动方程通过傅里叶变换转换为频率域后更容易添加衰减项，地震波的衰减更容易实现。因此频率-空间域数值模拟有着巨大的发展潜力和广泛的应用前景，为后续的深入研究和勘探提供理论依据。

地震波场数值模拟方法根据实现域的不同主要分为时间-空间域和频率-空间域两种。其中，频率-空间域虽然发展较晚，但由于频率-空间域数值模拟方法的误差积累较小，适合多炮计算，具有较大的发展潜力。

频率-空间域数值模拟最早可以追溯到20世纪70年代。1972年Lysmer和Drake首次实现了频率-空间域数值模拟。1984年Marfurt通过对比分析时间-空间域数值模拟和频率-空间域数值模拟的特点，指出频率-空间域数值模拟能更好的模拟地震波的传播和地震波的衰减效应。1990年Pratt和Worthington推导出频率-空间域声波方程差分格式和频率-空间域弹性波方程的5点有限差分格式。1996年Liao采用常规5点差分格式实现了频率-空间域粘滞声波的数值模拟和反演计算的研究。1996年Jo和Shin等人采用优化9 点差分格式实现了频率-空间域声波数值模拟计算，有效的压制了数值模拟过程中产生的频散现象。1998年Shin在Jo的基础上，将频率-空间域波动方程的差分格式扩展到25点有限差分格式，提高了数值模拟的精度。2002年Min等人提出并运用25点差分格式对弹性波进行了频率-空间域数值模拟，模拟结果表明 25点有限差分格式比较契合复杂介质中的地震波数值模拟，可以有效的缩减计算量，提高模拟精度。2004 年Hustede等人对比交错网格法和混合网格法发现混合网格法在内存需求和计算效率方面优于交错网格法。 2005年吴国忱等人实现了VTI介质和TTI介质中25点优化差分算子数值模拟。2009年任浩然等人采用二维声波方程最优化9点差分数值模拟，减小了计算存储空间，提高了数值模拟的计算速度，并研究讨论了边界条件和大型稀疏矩阵存储技术对数值模拟结果的影响。2011年李桂花等人实现了粘弹性VTI介质中频率-空间域25点最优化差分方法数值模拟，得到了qP波的传播特征。2014年李振春实现了起伏地表条件下的频率-空间域声波数值模拟，有效的压制了由于矩形网格离散产生的散射，使用变网格方法提高了计算效率。2022年吴国忱提出使用交错网格和混合网格进行频率-空间域非均质声波方程数值模拟，并验证了该方法的准确性和稳定性。

通过上述分析，现有技术存在的问题及缺陷为：

(1)现有技术的频率-空间域正演模拟方法采用常规的5点差分网格或9点差分网格，差分精度相对较低；

(2)现有技术的数值模拟方法是对声波的模拟，声波是把地下介质简化为各种同性介质情况下对地震波场的一种简化，但实际地下介质是各向异性的，而且还是粘弹性的，地震波场更为复杂，难以识别，只是模拟声波不能显示出野外实际资料中纵波的传播特征，qP波保留了各向异性特征中P波的特性，而地下绝大多数介质属于VTI介质，通过黏滞VTI介质波动方程推导出的qP波方程数值模拟得到的波场更符合实际地下介质中P波的传播规律，模拟的波场更有利于识别复介质中P波的识别和分离，为后期的处理方法的研究打好基础。

频率-空间域25点优化差分算法，可以提高数值模拟的精度。

但是，现有技术黏滞性VTI介质频率、空间域qР波25点优化差分数值算法在边界条件模拟上，模拟精度有待提高，而且对边界反射的各种波吸收效果差。

发明内容

为克服相关技术中存在的问题，本发明公开实施例提供了一种黏滞各向异性介质准qP波模拟方法、系统、设备及应用。具体涉及一种黏滞各向异性介质qP波数值模拟方法。

所述技术方案如下：一种黏滞各向异性介质qP波数值模拟方法包括以下步骤：

所述黏滞各向异性介质qP波数值模拟方法包括以下步骤：

基于时间-空间域VTI介质波动方程，通过傅里叶变换将时间-空间域VTI介质波动方程转换到频率- 空间域，结合Thomsen参数和Alkhalifah声波假设，推导出粘弹性VTI介质中qP波波动方程，消除横波干扰；

采用有限差分方法将地震模型离散化，剖分为多个的网格，用网格点速度替代速度模型；

构建空间域25点的优化差分格式，采用高斯牛顿迭代法计算差分系数，对频率域每个频率片并行进行频率-空间域的25点有限差分运算，对每个网格点都进行差分，结合震源项G，建立矩阵方程组；求解矩阵方程组得对应频率的空间各点的波场值即某频率切片，傅氏变换后变换到时间域得地震炮集记录。

在一个实施例中，采用数值模拟的计算方程如公式(1)：

式(1)中，

F为 qP波的波场；ω为圆频率；V_p0为qP波在垂直方向上传播的速度；χ＝1+2ε，η＝ε-δ，其中ε、γ、δ为衡量VTI介质中各向异性程度的Thomsen参数，没有单位，是地质模型输入信息，是已知的，ε表征纵波的各向异性强弱程度，γ表征横波的各向异性强弱程度，ε和γ同时增加或减少或者同时为零，δ为影响qP波传播速度的参数。

在一个实施例中，所述黏滞各向异性介质qP波数值模拟方法具体包括以下步骤：

S1，读取输入数据参数，包括初始地质模型文件，震源子波的类型，震源子波的主频，震源子波的衰减系数；

S2，确定模型剖分网格空间间隔Δx，Δz，及吸收边界的厚度；

S3，根据时间采样率Δt、道长Lt，根据地震数据在地表激发和接收规则，检波器的深度、位置参数，确定检波点数量，定义地面接收观测系统；

S4，根据震源子波类型和子波主频及震源子波的衰减系数生成多种震源子波；对确定生成的子波进行快速傅里叶变化，把震源子波变换到频率域，并将傅里叶变换后的实部和虚部采样点个数的前N/2+1 个点的值作为震源的实部和虚部放置在震源点位置作为震源项；

S5，在进行频率-空间域差分前先用Gauss-Newton法求解25点差分算子的最优加权系数；

S6，定义归一化速度为差分方程相速度与波动方程相速度的比值

归一化相速度表征频散程度，f的值接近1；

S7，对频率域震源项涉及的每个频率片并行进行频率-空间域的25点有限差分运算，对每个网格点都进行差分，结合震源项G，建立矩阵方程组；求解矩阵方程组得某频率对应的频率切片；

S8，采用在计算区域设置人工区域，在人工区域采用CPML法和特征分析法相组合的方法，吸收边界反射的影响；

S9，行运算出每个频率片的各个节点处的波场值，进行傅里叶变换，变换到时间域得到不同时刻时间域的波场值，输出波场时间切片结果和地面地震记录结果。

本发明的另一目的在于提供一种黏滞各向异性介质qP波数值模拟系统包括：

qP波波动方程获取模块，用于基于时间-空间域VTI介质波动方程，通过傅里叶变换将时间-空间域 VTI介质波动方程转换到频率-空间域，结合Thomsen参数和Alkhalifah声波假设，推导出粘弹性VTI介质中qP波波动方程，消除横波干扰；

网格剖分模块，用于采用有限差分方法将地震模型离散化，剖分为多个的网格，用网格点替代速度模型；

数值模拟计算模块，用于通过构建频率-空间域25点的波动方程的优化差分格式进行数值模拟计算。

本发明的另一目的在于提供一种计算机设备，所述计算机设备包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行所述黏滞各向异性介质qP波数值模拟方法。

本发明的另一目的在于提供一种所述黏滞各向异性介质qP波数值模拟方法在模拟地震波衰减的地震记录、地震波不衰减的地震记录、各向异性介质的波场快照和地震记录、各向同性介质的波场快照和模拟记录上的应用。

结合上述的所有技术方案，本发明所具备的优点及积极效果为：

第一、针对上述现有技术存在的技术问题以及解决该问题的难度，紧密结合本发明的所要保护的技术方案以及研发过程中结果和数据等，详细、深刻地分析本发明技术方案如何解决的技术问题，解决问题之后带来的一些具备创造性的技术效果。具体描述如下：

本发明是在黏滞性VTI介质频率-空间域qР波25点优化差分数值模拟的基础上，进一步优化了算法，边界条件，使得模拟精度更高，效果更好。

本发明给出了频率-空间域qP波波动方程，对qP波波动方程进行频率-空间域25点优化差分算法的有限差分数值模拟，分析并讨论了CPML和特征分析法组合边界对边界反射的吸收效果，模拟的qP波比声波来说，不仅更符合实际地下介质是各向异性的也考虑了地下介质的黏滞性性质，更符合地震波在地下传播能量会衰减的特性，模拟的地震波场精度更高，效果更好。

现有的方法中，有限差分法有在时间-空间域数值模拟的，也有在频率域-空间域进行数值模拟的，时间 -空间域数值模拟在每一个步进过程中会累积计算误差，而在频率-空间域数值模拟，具有不会累加误差的优点。

本发明是一种高精度的数值模拟方法，该数值模拟算法首先通过建立具有各向异性性质的地球物理模型，并且设置吸收衰减系数，模拟去掉横波的更符合实际地震波在地下介质中P波的传播规律，有利于进行P波的理论研究和对比分析野外地震数据。

本发明在模型范围内设置炮点坐标及各检波器坐标，在炮点给定震源激发函数，从VTI介质弹性波波动方程出发，推导出频率-空间域qP波波动方程，并构建25点高精度有限差分格式，利用高斯-牛顿法计算VTI介质差分算子的优化加权系数，能够有效地提高数值模拟的精度。

本发明为了解决人工边界造成的边界反射，通过对数值模型设置组合边界及特征分析法边界和CPML 组合边界能有效的吸收边界反射的各种波，这种组合的人工吸收边界条件对边界反射的吸收效果非常显著。

本发明因频率域数值模拟，实现了各频率片之间的并行计算，大大提高了计算效率。

第二、把技术方案看作一个整体或者从产品的角度，本发明所要保护的技术方案具备的技术效果和优点，具体描述如下：

本发明具体涉及VTI介质频率-空间域qP波波动方程、25点有限差分格式、特征分析法和CPML边界的组合边界技术，尤其涉及一种VTI介质频率-空间域25点优化插值方法，可用于高精度的构造复杂地区地震波的波场数值模拟，并获得地震记录及波场快照，用于研究地震波的传播规律与波场特征。

相比于现有技术，本发明的优点还包括：

本发明是一种波动方程qP波数值模拟方法。因一般情况下，地震勘探都是纵波勘探，而本模拟方法模拟的qP波，消除了横波干扰，可以有效地指导野外地震记录波场的识别，可以更好地被用于纵波地震波的成像及其他方法的研究；建立的VTI介质模型，其VTI的性质更符合实际地下沉积地层多为横向各项同向(即VTI)型介质性质，另外模型的黏滞性性质模拟的地震波更符合实际地震波在地下传播振幅和相位会衰减变化的情况；采用的频率域有限差分算法相比时间域的有限差分不仅不会累积计算误差，而且每个频率都可以采用并行算法并行运算提高计算的效率；有限差分采用25点优化差分算法，采用了牛顿-赛德尔迭代法计算优化差分系数，25点差分增强了该算法的精度，具有高精度的优点；并采用特征分析法和复完全匹配层技术相结合的组合边界对人工边界进行吸收，以压制边界反射效应。该数值模拟方法在提高数值模拟精度同时，又可以进一步提高准纵波模拟的计算效率。将本发明方法具体应用到实际的地质模型当中，与实际的地震记录相比，本发明具有很高的数值模拟精度。通过准纵波频率-空间域波动方程25点优化差分算法的数值模拟可以研究地震波的传播特征，可以指导波场的识别和纵波波场的分离，同时也可以为反演、偏移成像等后续工作提供理论基础。

附图说明

此处的附图被并入说明书中并构成本说明书的一部分，示出了符合本公开的实施例，并与说明书一起用于解释本公开的原理。

图1是本发明实施例提供的黏滞各向异性介质qP波数值模拟方法流程图；

图2是本发明实施例提供的黏滞各向异性介质qP波数值模拟方法原理图；

图3是本发明实施例提供的关于

近似的优化差分算子模型示意图；

图4是本发明实施例提供的关于

近似的优化差分算子模型示意图；

图5是本发明实施例提供的关于

近似的优化差分算子模型示意图；

图6是本发明实施例提供的关于ω⁴F近似的优化差分算子示意图；

图7是本发明实施例提供的CPML边界示意图；

图8(a)是本发明实施例提供的未加边界时50Hz单频波快照图；

图8(b)是本发明实施例提供的150ms时间域波前快照图；

图9(a)是本发明实施例提供的施加组合边界条件后的50Hz单频波快照图；

图9(b)是本发明实施例提供的150ms时间域波前快照图；

图10是本发明实施例提供的水平层状(多层)模型图；

图11(a)是本发明实施例提供的水平层状(多层)模型70ms波场快照图；

图11(b)是本发明实施例提供的水平层状(多层)模型100ms波场快照图；

图11(c)是本发明实施例提供的水平层状(多层)模型120ms波场快照图；

图11(d)是本发明实施例提供的水平层状(多层)模型150ms波场快照图；

图12是本发明实施例提供的洼陷模型图；

图13(a)是本发明实施例提供的洼陷模型70ms波场快照图；

图13(b)是本发明实施例提供的洼陷模型100ms波场快照图；

图13(c)是本发明实施例提供的洼陷模型120ms波场快照图；

图13(d)是本发明实施例提供的洼陷模型150ms波场快照图；

图14是本发明实施例提供的断层模型图；

图15(a)是本发明实施例提供的断层模型70ms波场快照图；

图15(b)是本发明实施例提供的断层模型100ms波场快照图；

图15(c)是本发明实施例提供的断层模型120ms波场快照图；

图15(d)是本发明实施例提供的断层模型150ms波场快照图；

图16是本发明实施例提供的井间地震共炮点道集的Z分量剖面图；

图17(a)是本发明实施例提供的优化差分算法模拟结果图；

图17(b)是本发明实施例提供的交错网格有限差分算法模拟结果图；

图18是本发明实施例提供的水平层状介质模型图；

图19(a)是本发明实施例提供的水平层状介质模型50Hz单频波快照对比中弹性介质图；

图19(b)是本发明实施例提供的水平层状介质模型50Hz单频波快照对比中粘弹性介质图；

图20(a)是本发明实施例提供的弹性和粘弹性介质qP波井间地震记录的对比中弹性介质图；

图20(b)是本发明实施例提供的弹性和粘弹性介质qP波井间地震记录的对比中粘弹性介质图；

图21(a)是本发明实施例提供的弹性和粘弹性远离震源检波器接收(首道)记录的频率振幅图的对比中弹性介质图；

图21(b)是本发明实施例提供的弹性和粘弹性远离震源检波器接收(首道)记录的频率振幅图的对比中粘弹性介质图；

图22(a)是本发明实施例提供的均匀各向异性介质频率域50Hz单频波快照对比图中常规9点差分算子图；

图22(b)是本发明实施例提供的均匀各向异性介质频率域50Hz单频波快照对比图中25点优化差分算子图；

图23(a)是本发明实施例提供的均匀各向异性介质时间域150ms波场快照对比图中常规9点差分算子图；

图23(b)是本发明实施例提供的均匀各向异性介质时间域150ms波场快照对比图中25点优化差分算子图；

图24(a)是本发明实施例提供的均匀各向异性介质地震记录的对比图中常规9点差分算子图；

图24(b)是本发明实施例提供的均匀各向异性介质地震记录的对比图中25点优化差分算子图。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图对本发明的具体实施方式做详细的说明。在下面的描述中阐述了很多具体细节以便于充分理解本发明。但是本发明能够以很多不同于在此描述的其他方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似改进，因此本发明不受下面公开的具体实施的限制。

一、解释说明实施例：

本发明实施例提供一种黏滞各向异性介质qP波数值模拟方法包括：

所述黏滞各向异性介质qP波数值模拟方法包括以下步骤：

基于时间-空间域VTI介质波动方程，通过傅里叶变换将时间-空间域VTI介质波动方程转换到频率- 空间域，结合Thomsen参数和Alkhalifah声波假设，推导出粘弹性VTI介质中qP波波动方程，消除横波干扰；

采用有限差分方法将地震模型离散化，剖分为多个的网格，用网格点速度替代速度模型；

构建空间域25点的优化差分格式，采用高斯牛顿迭代法计算差分系数，对频率域每个频率片并行进行频率-空间域的25点有限差分运算，对每个网格点都进行差分，结合震源项G，建立矩阵方程组；求解矩阵方程组得对应频率的空间各点的波场值即某频率切片，傅氏变换后变换到时间域得地震炮集记录

具体地，如图1所示，本发明实施例提供的黏滞各向异性介质qP波数值模拟方法包括以下步骤：

S101，基于时间-空间域VTI介质波动方程，通过傅里叶变换将时间-空间域VTI介质波动方程转换到频率-空间域，结合Thomsen参数和Alkhalifah声波假设思想，推导出粘弹性VTI介质中qP波波动方程，消除横波干扰；

S102，采用有限差分方法将地震模型离散化，剖分为一个个的网格，用网格点替代速度模型；

S103，通过构建频率-空间域25点的波动方程的优化差分格式进行数值模拟计算。

在本发明实施例中，提供了一种高效率-高精度的频率-空间域VTI介质准纵波波动方程数值模拟方法，分析了精度较高的频率-空间域准纵波25点优化有限差分格式，迭代计算了差分系数，并采用吸收效果较好的特征分析法和复最佳匹配层的组合边界，设计建立黏滞性VTI介质模型，进行不同构造形态模型的数值模拟来计算获得地面地震记录及其波场快照，分析地震波的传播规律及其波场特征，指导波场的识别和纵波波场的分离，同时也可以为反演、偏移成像等后续工作提供理论基础。

实施例1

本发明实施例提供一种频率-空间域的准纵波数值模拟方法，为了解决各向异性介质复杂构造地区波场复杂，波场复杂不易识别且成像难等问题，本数值模拟方法可以进行基于VTI介质的准纵波波场的模拟，更符合实际介质中纵波的波场。本发明在频率-空间域进行差分计算，具有时间域所不具有的优点，不仅不会累计计算误差，而且可以进行并行计算，提高计算效率，而且采用25点优化差分算法大大地提高了计算精度。同时提出特征分析法吸收边界和复完全匹配层组成组合边界对人工边界反射进行压制。通过实施例验证本发明的正确性，并将本发明应用到某油田地区的实际模型当中，与实际地震数据进行了对比，对该油田地区地震波的传播机制有更好的解释。本发明需要输入地质模型，震源子波类型，吸收边界厚度，观测系统的设置参数，就可以实现qP波波场的数值模拟。

本发明采用数值模拟的计算方程如公式(1)：

式(1)中，

F为qP波的波场；ω为圆频率；V_p0为qP波在垂直方向上传播的速度；χ＝1+2ε，η＝ε-δ，其中ε、γ、δ为衡量VTI介质中各向异性程度的Thomsen参数，没有单位，是地质模型输入信息，是已知的，ε表征纵波的各向异性强弱程度，γ表征横波的各向异性强弱程度，ε和γ同时增加或减少或者同时为零，δ为影响qP波传播速度的参数。

用25点差分格式进行频率-空间域VTI介质qP波数值模拟可以提高数值模拟精度，降低频散。

近似的25点优化差分格式示意图如图3所示。每行用5个网格点构造2个二阶差分算子，其中，用第1项、第3项、第5项(画圈)构造间距为2△x的差分算子，第2项、第3项、第4项构造间距为△x的差分算子，每行的两个差分算子用加权系数C和D平均，这样共得到5个差分算子，再根据与中间行的距离不同分别用加权系数B₁、B₂、B₃平均，中间行的差分算子用加权系数B₁进行加权平均，和中间行距离为△z 的两个差分算子用B₂进行加权平均，和中间行距离为2△z的两个差分算子用B₃进行加权平均。最后得到了

近似的25点优化差分格式表达式为公式(2)。同理，

近似的25点优化差分格式示意图如图4 所示，表达式为公式(3)。

的近似优化差分格式如图5所示。将25点优化差分算子代入到

的差分算子中构造间距为△x、△z和2△x、2△z的2个差分算子，每个差分算子都是由9项构成的四阶差分算子，再用加权系数e 和f对这两个四阶差分算子进行加权平均，其表达式为式(4)。

对于ω⁴F项，将质量加速度算子分配到以包含所求网格点在内的共25个网格点上，采用加权系数a₁对所求网格点进行加权，按照距离所求中心网格点的远近分别采用加权系数a₂，a₃，a₄，a₅，a₆对其余的网格点进行加权平均，与中心网格点距离为△x、△z的网格点用加权系数a₂进行加权平均，与中心网格点距离为

的网格点用加权系数a₃进行加权平均，与中心网格点距离为2△x、2△z的网格点用加权系数a₄进行加权平均，与中心网格点距离为

的网格点用加权系数a₅进行加权平均，与中心网格点距离为

的网格点用加权系数a₆进行加权平均，用来近似质量加速度项所用的网格点示意图如图6所示。

将25点差分算子代入F_i,j，可得：

其中，

在二维情况下，设模型的X方向上有nx个空间采样点，Z方向上有nz个空间采样点，F为对应的地震波场值。根据上述求得的F_i,j项的差分格式，对每个网格点F_i,j都可建立一个这样的方程，结合震源项G，可建立矩阵方程组的表达式为：

A_nx×nzF＝G (7)

其中：

F＝(F_0,0,F_1,0,…,F_nx-1,0,F_0,1,…,F_nx-1,1,…,F_nx-1,nz-1)^T

G＝(G_0,0,G_1,0,…,G_nx-1,0,G_0,1,…,G_nx-1,1,…,G_nx-1,nz-1)^T

用迭代法求解式(7)就可以得每个计算频率对应的空间网格点上的波场值F。

实施例2

如图2所示，本发明实施例提供的一种频率空间域qP波高精度地震波场数值模拟的方法包含以下步骤:

步骤1，读取输入数据参数，包括初始地质模型文件，震源子波的类型，震源子波的主频，震源子波的衰减系数。

步骤2，确定模型剖分网格空间间隔Δx，Δz，及吸收边界的厚度。

步骤3，根据时间采样率Δt、道长Lt，根据地震数据在地表激发和接收规则，检波器的深度、位置等参数，确定检波点数量，定义好地面接收观测系统。

步骤4，根据震源子波类型和子波主频及震源子波的衰减系数生成震源子波，可以生成Gabor子波、 Berlage波、穆勒子波：雷克子波、高斯子波、带通子波、最小相位子波7种类型子波。对确定生成的子波进行快速傅里叶变化，把震源子波变换到频率域，并将傅里叶变换后的实部和虚部采样点个数的前 N/2+1个点的值作为震源的实部和虚部放置在震源点位置作为震源项。

步骤5，在进行频率-空间域差分之前需要首先用Gauss-Newton法求解25点差分算子的最优加权系数得

a₁＝0.41462，a₂＝0.1178，a₃＝0.018237

a₄＝-0.0018609，a₅＝0.002508，a₆＝-0.00046843

B₁＝0.6098，B₂＝0.15388，B₃＝0.0011674

C＝0.66482，D＝0.39359，e＝0.57896，f＝0.39843

步骤6，定义归一化速度为差分方程相速度与波动方程相速度的比值

归一化相速度表征频散程度，f的值越接近1，说明频散程度越低。

步骤7，对频率域震源项涉及的每个频率片并行进行频率-空间域的25点有限差分运算，具体差分公式见式(6)，对每个网格点都按照(6)式进行差分，结合震源项G，可建立矩阵方程组的表达式(7)。求解(7)式得某频率对应的频率切片。

步骤8，为了压制边界反射，本发明实施例采用了在计算区域设置一个人工区域，在人工区域采用 CPML法和特征分析法相组合的方法，当地震波传播到这个吸收边界时，地震波经过了两种方法的联合吸收衰减，很大程度上减少了人工边界反射的干扰，从而有效的吸收了边界反射的影响。其中CPML边界的复频伸展函数为：

式(8)中：κ_x、κ_z为增加的辅助衰减系数，且κ_x≥1，κ_z≥¹，a_x≥⁰，a_z≥0。如果设定κ_x＝0，κ_z＝0，a_x＝1，a_z＝1，方程(8)就会变成传统的PML吸收边界条件。

图7为在计算区域的周围施加CPML匹配层边界示意图，其中，

标号④的区域为非匹配层区域，标号第一匹配层区域①、第二匹配层区域②、第三匹配层区域③为匹配层区域；

在①区域内，σ_x≠0，σ_z≠0；在②区域内，σ_x＝0，σ_z≠0；在③区域内，σ_x≠0，σ_z＝0。

结合衰减函数σ_x、σ_z和VTI介质弹性波动方程，得到二维VTI介质qP波波动方程为：

另外，特征分析法同时加载在匹配层区域上，对应区域上的Qp波波动方程见式(10)：

式(10)中，

步骤9，并行运算出每个频率片的各个节点处的波场值，进行傅里叶变换，变换到时间域得到不同时刻时间域的波场值，输出波场时间切片结果和地面地震记录结果。

实施例3

基于实施例2记载的频率空间域qP波高精度地震波场数值模拟的方法，为了让本发明的目的、技术和优势更加清晰，下面将通过原理分析以及结合附图作进一步的说明。

本发明实施例根据实际地层中比较常见的几种地层结构设计了均匀介质模型，水平层状模型，断层模型，凹陷模型，并对野外采集的实际地震记录，通过本发明实施例建立模型模拟的记录与野外记录进行对比。

本发明实施例在建立一个简单的均匀介质模型进行数值正演模拟来验证边界条件的吸收效果，模型参数设置为：模型大小700m×700m，qP波速度v_p0＝2500m/s，网格间距10m，采样间隔2ms，采样点数512，震源位于(350m，350m)，采用50Hz雷克子波，边界厚度100m，ρ为1.7，ε为0.255，δ为0.32，在匹配层区域采用施加CPML吸收边界和特征分析法边界组合边界来压制边界反射。

图8(a)是本发明实施例提供的未加边界时50Hz单频波快照图；

图8(b)是本发明实施例提供的150ms时间域波前快照图；

由图8(a)、图8(b)可以看出，未施加边界条件时产生了强烈的边界反射，由于入射波和产生的边界反射相互影响，单频波快照中波形抖动剧烈，时间域波前快照中无法明确地显示出波前。

图9(a)是本发明实施例提供的施加组合边界条件后的50Hz单频波快照图；

图9(b)是本发明实施例提供的150ms时间域波前快照图。

由图9(a)、图9(b)可以看出，当采用本发明实施例所提的组合边界时，单频波快照更清晰光滑，时间域波前快照上基本看不到边界反射回来的能量。

(1)水平层状模型

建立的水平层状模型示意图如图10所示，模型大小为700m×700m，震源位于模型中心(350m，350m)，采样间隔为2ms，采样点数为512点，组合吸收边界厚度为100m，各向异性参数：ρ为1.7，ε为0.255，δ为0.32，横波速度为0，纵波速度为500m/s、2200m/s、800m/s、2700m/s、2800m/s，边界条件设定为 CPML边界和特征分析法组合边界条件。选择网格间距为5m，采用30Hz雷克子波作为震源主频作为所设定数值模型的参数通过模型进行验证VTI介质频率-空间域qP波有限差分数值模拟。

其中，图11(a)是本发明实施例提供的水平层状(多层)模型70ms波场快照图；图11(b)是本发明实施例提供的水平层状(多层)模型100ms波场快照图；图11(c)是本发明实施例提供的水平层状(多层)模型120ms波场快照图；图11(d)是本发明实施例提供的水平层状(多层)模型150ms波场快照图；

从图11(a)-图11(d)的波场快照中可以看出，波前在一圈一圈的前进，由于模型上部存在低速区域，地震波在低速带中传播速度一般远低于其下伏地层的速度。所以模型上部分的波前传播的比下部分慢很多。在150ms的波场快照中可以看出，波前传播到边界时没有产生边界反射。

(2)洼陷模型

建立洼陷模型如图12所示，洼陷模型共有三个反射面，第一个和第三个反射面是水平的，第二个反射面是是由两个不同深度的水平面和两个倾斜面组成的。模型大小为700m×700m，采样间隔为2ms，采样点数为512，边界厚度为100m，各向异性参数：ρ为1.7，ε为0.255，δ为0.32，边界条件设定为CPML 边界和特征分析法组合边界条件。通过对洼陷模型对频率-空间域qP波25点有限差分数值模拟进行验证，如图13(a)是本发明实施例提供的洼陷模型70ms波场快照图；图13(b)是本发明实施例提供的洼陷模型100ms波场快照图；图13(c)是本发明实施例提供的洼陷模型120ms波场快照图；图13(d)是本发明实施例提供的洼陷模型150ms波场快照图；

从图13(a)-图13(d)波场快照中可以看出没有出现明显的波前抖动现象，同时也可以看到地震波遇到第三层反射层面后产生的反射波，而且边界的吸收效果很好，没有边界反射。也可以看到地震波遇到洼陷部分反射出来的波。

(3)断层模型

断层模型示意图如图14所示，模型大小为700m×700m，模型的其他参数设置为：采样间隔为2ms，采样点数为512，边界厚度为100m，各向异性参数：ρ为1.7，ε为0.255，δ为0.32，横波速度为0，边界条件设定为CPML边界和特征分析法组合边界条件。模型内介质的纵波速度V_p0从上到下依次是 2200m/s，2000m/s，2800m/s，2700m/s。如图15(a)是本发明实施例提供的断层模型70ms波场快照图；图15(b)是本发明实施例提供的断层模型100ms波场快照图；图15(c)是本发明实施例提供的断层模型120ms波场快照图；图15(d)是本发明实施例提供的断层模型150ms波场快照图所示。由于模型内存在垂直断层和倾斜的地层，从图15(a)-图15(d)中可以看出，从150ms的波场快照中可以清晰的看出两条反射波，一条向左上角，一条水平向右，可以清晰的判断出断层和倾斜面的位置。

(4)野外采集的实际记录的模拟

野外采集的井间地震记录，井间距299m，炮点深度1402m，检波器深度从986m到1811m，检波器间隔 3m，共276个检波点。地震记录时间采样间隔为0.25ms，记录长度1s。图16为原始采集的井间地震Z分量记录。

图17(a)是本发明实施例提供的优化差分算法模拟结果图；

图17(b)是本发明实施例提供的交错网格有限差分算法模拟结果图；从图17(a)、图17(b)中剖面上可以看出该记录噪声背景比较大，反射波场被噪音所压制，q^p初至1很清晰，即使接收井内与炮点深度对应的检波点所接收到的初至波也隐约可见，除了纵波的初至外，在250ms～400ms之间有几组能量很强的同相轴，这些是横波的信息，经分析对比，确定图中3为慢横波q^SV的直达波，4为慢横波q^SV的下行反射波，5为慢横波q^SV的上行反射波，由于其能量很强，几乎压制了q^p初至外的其他所有的纵波信息。

表1根据测井的声波曲线和密度测井曲线求得的各层的纵波速度和密度

利用测井曲线获得如表1所示地层的速度和密度信息，并结合初至波或直达波信息，计算初至旅行时或直达波旅行时，纵波初至波与实际记录的q^p初至数值相等时，即确定了慢纵波(即0°纵波)和快纵波 (即90°纵波)的比值。根据实际资料的分析和计算，在数值模拟中采用了慢纵波和快纵波速度比0.715，从而依次求得各网格点的快慢纵波速度及密度，建立了各向异性的地质模型。

利用本发明频率空间域25点优化差分算法模拟的井间地震qP记录的剖面，其模拟参数：时间采样间隔0.5ms，空间采样间隔3m。震源子波采用80Hz的雷克子波。利用时间域交错网格模拟的Z分量剖面，其模拟参数分别为：时间采样间隔0.25ms，空间采样间隔3m，震源震源子波采用80Hz的雷克子波。由于频率空间域对矩阵是整体求解，占用很大的内存量，受机器性能限制，在空间采样间隔相同的情况下，时间采样间隔仅能达到0.5ms采样。从图(16)对比可以看出这两幅图在P波直达波、P上行反射波和下行反射波的时间对应关系都是一样的，但(a)图中，由于仅有qP波的信息，所以波形相对简单，易于辨认波场。

实施例4

本发明实施例提供的黏滞各向异性介质qP波数值模拟方法，还包含在弹性特殊情况。Qp值趋于无穷大时，趋近于弹性各向异性介质，通过调整各向异性系数，各向异性介质也既可以变为弹性各项同向介质。

实施例5

本发明实施例提供的黏滞各向异性介质qP波数值模拟方法不仅可以模拟地震波衰减的地震记录；也可以模拟地震波不衰减的地震记录。

实施例6

本发明实施例提供的黏滞各向异性介质qP波数值模拟方法不仅可以模拟各向异性介质的波场快照和地震记录，也可以模拟各向同性介质的波场快照和模拟记录。

本发明实施例除了本案例里的模型，可以任意给出不同结构的模型进行数值模拟。

实施例7

本发明实施例提供的黏滞各向异性介质qP波数值模拟方法采用了25点优化差分算法，若计算精度不需要太高，也可以采用9点常规差分算法的模拟算法，传统的9点差分就是把25点外圈的16点去掉，只采用中间的9个点来进行差分，来降低计算量。

实施例8

本发明实施例提供一种黏滞各向异性介质qP波数值模拟系统包括：

qP波波动方程获取模块，用于基于时间-空间域VTI介质波动方程，通过傅里叶变换将时间-空间域 VTI介质波动方程转换到频率-空间域，结合Thomsen参数和Alkhalifah声波假设，推导出粘弹性VTI介质中qP波波动方程，消除横波干扰；

网格剖分模块，用于采用有限差分方法将地震模型离散化，剖分为多个的网格，用网格点替代速度模型；

数值模拟计算模块，用于通过构建频率-空间域25点的波动方程的优化差分格式进行数值模拟计算。

在上述实施例中，对各个实施例的描述都各有侧重，某个实施例中没有详述或记载的部分，可以参见其它实施例的相关描述。

上述装置/单元之间的信息交互、执行过程等内容，由于与本发明方法实施例基于同一构思，其具体功能及带来的技术效果，具体可参见方法实施例部分，此处不再赘述。

所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，仅以上述各功能单元、模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能单元、模块完成，即将所述装置的内部结构划分成不同的功能单元或模块，以完成以上描述的全部或者部分功能。实施例中的各功能单元、模块可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中，上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。另外，各功能单元、模块的具体名称也只是为了便于相互区分，并不用于限制本发明的保护范围。上述系统中单元、模块的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。

二、应用实施例：

应用例1

本发明实施例提供的黏滞各向异性介质qP波数值模拟方法涉及VTI介质频率-空间域qP波波动方程、 25点有限差分格式、特征分析法和CPML边界的组合边界技术，尤其涉及一种VTI介质频率-空间域25 点优化插值方法，可用于高精度的构造复杂地区地震波的波场数值模拟，并获得地震记录及波场快照，用于研究地震波的传播规律与波场特征。

应用例2

本发明应用实施例还提供了一种计算机设备，该计算机设备包括：至少一个处理器、存储器以及存储在所述存储器中并可在所述至少一个处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述任意各个方法实施例中的步骤。

应用例3

本发明应用实施例还提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时可实现上述各个方法实施例中的步骤。

应用例4

本发明应用实施例还提供了一种信息数据处理终端，所述信息数据处理终端用于实现于电子装置上执行时，提供用户输入接口以实施如上述各方法实施例中的步骤，所述信息数据处理终端不限于手机、电脑、交换机。

应用例5

本发明应用实施例还提供了一种服务器，所述服务器用于实现于电子装置上执行时，提供用户输入接口以实施如上述各方法实施例中的步骤。

应用例6

本发明应用实施例提供了一种计算机程序产品，当计算机程序产品在电子设备上运行时，使得电子设备执行时可实现上述各个方法实施例中的步骤。

所述集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明实现上述实施例方法中的全部或部分流程，可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一计算机可读存储介质中，该计算机程序在被处理器执行时，可实现上述各个方法实施例的步骤。其中，所述计算机程序包括计算机程序代码，所述计算机程序代码可以为源代码形式、对象代码形式、可执行文件或某些中间形式等。所述计算机可读介质至少可以包括：能够将计算机程序代码携带到拍照装置/终端设备的任何实体或装置、记录介质、计算机存储器、只读存储器(Read-OnlyMemory，ROM)、随机存取存储器(RandomAccessMemory，RAM)、电载波信号、电信信号以及软件分发介质。例如U盘、移动硬盘、磁碟或者光盘等。

三、实施例相关效果的证据：

实验

本发明实施例提供的模拟方法虽然是一种VTI各向异性粘弹性介质的地震数值模拟方法，但弹性和各项同向是一种简化，所以本发明可以模拟的波场包括粘弹性介质、弹性介质、VTI各向异性介质、各项同向介质，当采用不用的模型参数，可以模型出不同的介质模型的波场快照和地震记录。为了说明其不同的效果，利用图18所示水平层状模型结构水意图分别对弹性介质和粘弹性介质进行了井间地震波前快照和地震记录的数值模拟试算。模型共有5层，水平层状介质模型的X方向范围从0m-1000m，Z方向范围从 0m-1000m，空间采样间隔△x＝△z＝10m，模型网格101×101。时间采样间隔2.0ms，震源采用主频为 50Hz的雷克子波，震源位置为(0m，500m)，接收波波器布置在如图18所示模型右侧黑粗线。对于弹性介质模型参数见表1，粘弹性介质模型参数见表2。图19(a)50Hz单频波快照对比中水平层状介质模型弹性介质；图19(b)50Hz单频波快照对比中粘弹性介质；从图19(a)、图19(b)上可以看出，其单频波快照在形态上是一样的，只是数值上不同，从图的灰度颜色深浅上可以看出在单频波快照上，粘弹性介质情况下的振幅值比弹性介质振幅值小，灰度颜色变浅。图20(a)为模拟生成的弹性介质qP波井间地震记录，图20(b)为模拟生成的粘弹性介质井间地震记录，从图18中可以看出，同相轴在单位时间内波形变宽，说明频率降低。

图21(a)是本发明实施例提供的弹性和粘弹性远离震源检波器接收(首道)记录的频率振幅图的对比中弹性介质图；即模拟的弹性VTI和粘弹性VTI介质一个检波器接收到的一道记录的频率振幅对比图；

图21(b)是本发明实施例提供的弹性和粘弹性远离震源检波器接收(首道)记录的频率振幅图的对比中粘弹性介质图；

图21(a)-图21(b)中上面曲线是整个地震记录所有道的频率振幅平均值，下面曲线是该道的频率振幅值。从图中可以看出，对于同一道数据，在粘弹性介质模拟的记录的振幅明显比弹性介质中模拟的记录小的多，在粘弹性介质中，高频成分的震源明显的降低，低频成分的震幅值降低的相对比较小，整个记录主频向低频移动。

表1对图16中水平层状介质模型弹性介质参数表

表2对图16中水平层状介质模型粘弹性介质参数表

在本发明实施例中，为了说明满足精度要求的情况下为了减少计算工作量可以采用9点的插值方法，去掉25点外圈的16个点，只采用中间的3*3网格点进行插值，为了区别两者的效果，应用本发明建立一个均匀VTI介质模型，介质的垂向速度为υ_p0＝2500m/s，VTI介质的各向异性参数ε＝0.189，δ＝0.204。模型网格为101×101，空间间隔△x＝△z＝10m，时间采样间隔2ms，采用主频为50Hz的雷克子波作为震源子波，震源位于模型中心处。

图22(a)是本发明实施例提供的均匀各向异性介质频率域50Hz单频波快照对比图中常规9点差分算子图；图22(b)是本发明实施例提供的均匀各向异性介质频率域50Hz单频波快照对比图中25点优化差分算子图；

图23(a)是本发明实施例提供的均匀各向异性介质时间域150ms波场快照对比图中常规9点差分算子图；图23(b)是本发明实施例提供的均匀各向异性介质时间域150ms波场快照对比图中25点优化差分算子图；

图24(a)是本发明实施例提供的均匀各向异性介质地震记录的对比图中常规9点差分算子图；图24 (b)是本发明实施例提供的均匀各向异性介质地震记录的对比图中25点优化差分算子图；

从上图可以看出25点优化差分对常规的9点差分精度要高。不过若精度要求不是特别高，可以通过缩小空间网格间隔来提高精度，降低频散，用替代方案采用9点差分来进行数值模拟。

以上所述，仅为本发明较优的具体的实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，都应涵盖在本发明的保护范围之内。

Claims (10)

1.一种黏滞各向异性介质qP波数值模拟方法，其特征在于，该方法基于时间-空间域VTI介质波动方程，通过傅里叶变换将时间-空间域VTI介质波动方程转换到频率-空间域，结合Thomsen参数和Alkhalifah声波假设，推导出粘弹性VTI介质中qP波波动方程，消除横波干扰；

采用有限差分方法将地震模型离散化，剖分为多个的网格，用网格点速度替代速度模型；构建空间域25点的优化差分格式，采用高斯牛顿迭代法计算差分系数，对频率域每个频率片并行进行频率-空间域的25点有限差分运算，对每个网格点都进行差分，结合震源项G，建立矩阵方程组；求解矩阵方程组得对应频率的空间各点的波场值，即某频率切片，通过傅氏变换后变换到时间域得地震炮集记录。

2.根据权利要求1所述的黏滞各向异性介质qP波数值模拟方法，其特征在于，采用数值模拟的计算方程为：

式(1)中，

F为qP波的波场；ω为圆频率；V_p0为qP波在垂直方向上传播的速度；χ＝1+2ε，η＝ε-δ，其中ε、γ、δ为衡量VTI介质中各向异性程度的Thomsen已知参数，ε表征纵波的各向异性强弱程度，γ表征横波的各向异性强弱程度，ε和γ同时增加或减少或者同时为零，δ为影响qP波传播速度的参数。

3.根据权利要求1所述的黏滞各向异性介质qP波数值模拟方法，其特征在于，所述黏滞各向异性介质qP波数值模拟方法具体包括以下步骤：

S1，读取输入数据参数，该数据参数包括初始地质模型文件，震源子波的类型，震源子波的主频，震源子波的衰减系数；

S2，确定模型剖分网格空间间隔Δx，空间间隔Δz，及吸收边界的厚度；

S3，根据时间采样率Δt、道长Lt，根据地震数据在地表激发和接收规则，检波器的深度、位置参数，确定检波点数量，定义地面接收观测系统；

S4，根据震源子波类型和子波主频及震源子波的衰减系数生成多种震源子波；对确定生成的子波进行快速傅里叶变化，把震源子波变换到频率域，并将傅里叶变换后的实部和虚部采样点个数的前N/2+1个点的值作为震源的实部和虚部放置在震源点位置作为震源项；

S5，在进行频率-空间域差分前先用Gauss-Newton法求解25点差分算子的最优加权系数；

S6，定义归一化速度为差分方程相速度与波动方程相速度的比值

归一化相速度表征频散程度，f的值接近1；

S7，对频率域震源项涉及的每个频率片并行进行频率-空间域的25点有限差分运算，对每个网格点都进行差分，结合震源项G，建立矩阵方程组；求解矩阵方程组得某频率对应的频率切片；

S8，在计算区域设置人工区域，在人工区域采用CPML法和特征分析法相组合的方法，吸收边界反射的影响；

S9，运算出每个频率片的各个节点处的波场值，进行傅里叶变换，变换到时间域得到不同时刻时间域的波场值，输出波场时间切片结果和地面地震记录结果。

4.根据权利要求3所述的黏滞各向异性介质qP波数值模拟方法，其特征在于，在步骤S4中，震源子波包括：Gabor子波、Berlage波、穆勒子波、雷克子波、高斯子波、带通子波、最小相位子波七种类型子波。

5.根据权利要求3所述的黏滞各向异性介质qP波数值模拟方法，其特征在于，在步骤S5中，Gauss-Newton法求解得到空间网格差分的25点差分算子的最优加权系数：

a₁＝0.41462，a₂＝0.1178，a₃＝0.018237，a₄＝-0.0018609，a₅＝0.002508，

a₆＝-0.00046843；

B₁＝0.6098，B₂＝0.15388，B₃＝0.0011674；

C＝0.66482，D＝0.39359，e＝0.57896，f＝0.39843。

6.根据权利要求3所述的黏滞各向异性介质qP波数值模拟方法，其特征在于，在步骤S7中，对频率域震源项涉及的每个频率片并行进行频率-空间域的25点有限差分运算，对每个网格点都进行差分，结合震源项G，建立矩阵方程组；求解矩阵方程组得某频率对应的频率切片具体包括：

用25点差分格式进行频率-空间域VTI介质qP波数值模拟可以提高数值模拟精度，25点优化差分格式每行用5个网格点构造2个二阶差分算子；

其中，用第1项、第3项、第5项构造空间间距为2Δx的差分算子；

其中，Δx为模型剖分X方向的空间间隔；

第2项、第3项、第4项构造间距为Δx的差分算子，每行的两个差分算子用加权系数C和D平均，这样共得到5个差分算子，再根据与中间行的距离不同分别用加权差分系数B₁、B₂、B₃平均，中间行的差分算子用加权系数B₁进行加权平均，和中间行距离为Δz的两个差分算子用B₂进行加权平均，和中间行距离为2Δz的两个差分算子用B₃进行加权平均；

最后得到了

近似的25点优化差分格式表达式为公式(2)；

同理，

近似的25点优化差分格式表达式为公式(3)；

式(2)、式(3)中，F为(i，j)点的差分之后的qP波波场值，F下标不同代表(i，j)及其周围25个点不同点的qP波波场值，利用这25点的波场值，差分得F点对X方向的二阶偏导数

和F点对Y方向

的二阶偏导数；

将25点优化差分算子代入到

的差分算子中构造间距为Δx、Δz(Δz为模型剖分Z方向的空间间隔)和2Δx、2Δz的2个差分算子，每个差分算子都是由9项构成的四阶差分算子，再用加权系数e和f对这两个四阶差分算子进行加权平均，表达式为式(4)；

对于ω⁴F项，F为所求的网格点(i，j)处差分求得的qP波的波场值；ω为圆频率，将差分算子F_i,j分配到包含所求网格点(i，j)在内的共25个网格点上，采用加权系数a₁对所求网格点进行加权，按照距离所求中心网格点的远近分别采用加权系数a₂，a₃，a₄，a₅，a₆对其余的网格点进行加权平均，与中心网格点距离为Δx、Δz的网格点用加权系数a₂进行加权平均，与中心网格点距离为

的网格点用加权系数a₃进行加权平均，与中心网格点距离为2Δx、2Δz的网格点用加权系数a₄进行加权平均，与中心网格点距离为

的网格点用加权系数a₅进行加权平均，与中心网格点距离为

的网格点用加权系数a₆进行加权平均，用来近似网格点(i，j)的qP波波场值F_i,j；

将25点差分算子代入F_i,j，得：

其中，

在二维情况下，设模型的X方向上有nx个空间采样点，Z方向上有nz个空间采样点，F为对应的地震波场值；根据上述求得的F_i,j项的差分格式，对每个网格点F_i,j都能建立，结合震源项G，建立矩阵方程组的表达式为：

A_nx×nzF＝G (7)

其中：

F＝(F_0,0,F_1,0,…,F_nx-1,0,F_0,1,…,F_nx-1,1,…,F_nx-1,nz-1)^T

G＝(G_0,0,G_1,0,…,G_nx-1,0,G_0,1,…,G_nx-1,1,…,G_nx-1,nz-1)^T

用迭代法求解式(7)得每个计算频率对应的空间网格点上的qP波波场值F。

7.根据权利要求1所述的黏滞各向异性介质qP波数值模拟方法，其特征在于，在步骤S8中，CPML边界的复频伸展函数e_x,e_z为：

式(8)中：κ_x、κ_z为增加的辅助衰减系数，且κ_x≥1，κ_z≥1，实数部分a_x,a_z为尺度因子a_x≥0，a_z≥0；如果设定κ_x＝0，κ_z＝0，a_x＝1，a_z＝1，σ_x、σ_z分别为X和Z方向的衰减函数，式(8)就会变成传统的PML吸收边界条件；

在计算区域的周围施加CPML匹配层边界中，计算区域包括非匹配层区域、第一匹配层区域、第二匹配层区域、第三匹配层区域；

在第一匹配层区域内，σ_x≠0，σ_z≠0；在第二匹配层区域内，σ_x＝0，σ_z≠0；在第三匹配层区域内，σ_x≠0，σ_z＝0；

结合衰减函数σ_x、σ_z和VTI介质弹性波动方程，得到二维VTI介质qP波波动方程为：

特征分析法同时加载在匹配层区域上，对应区域上的Qp波波动

方程见式(10)：

式(10)中，

8.一种实现如权利要求1-7任意一项所述黏滞各向异性介质qP波数值模拟方法的系统，其特征在于，该黏滞各向异性介质qP波数值模拟系统包括：

qP波波动方程获取模块，用于基于时间-空间域VTI介质波动方程，通过傅里叶变换将时间-空间域VTI介质波动方程转换到频率-空间域，结合Thomsen参数和Alkhalifah声波假设，推导出粘弹性VTI介质中qP波波动方程，消除横波干扰；

网格剖分模块，用于采用有限差分方法将地震模型离散化，剖分为多个的网格，用网格点替代速度模型；

数值模拟计算模块，用于通过构建频率-空间域25点的波动方程的优化差分格式进行数值模拟计算。

9.一种计算机设备，其特征在于，所述计算机设备包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行权利要求1～7任意一项所述黏滞各向异性介质qP波数值模拟方法。

10.一种权利要求1～7任意一项所述黏滞各向异性介质qP波数值模拟方法在模拟地震波衰减的地震记录、地震波不衰减的地震记录、各向异性介质的波场快照和地震记录、各向同性介质的波场快照和模拟记录上的应用。

Images