基于对称波动方程的TTI介质全波形反演方法及装置

本技术涉及地球物理勘探，尤其是涉及一种基于对称波动方程的tti介质全波形反演方法及装置。

背景技术：

1、全波形反演是一种能够利用地震波全部波形信息（包括走时、振幅、相位及多波至等）定量重建地下介质弹性参数（速度、密度和各向异性参数）的勘探地球物理方法，为油气藏识别与评价提供高精度的定量模型。

2、全波形反演通过最小化观测地震记录与理论正演记录之间的失配，能够恢复地下介质的高分辨率参数模型，是地震勘探领域的一项关键技术。在地下构造相对简单的区域，各向同性或垂直横向各向同性（vti）假设能够较好地描述波传播特征，并已在工业界获得广泛应用。然而，对于复杂沉积盆地、倾斜层状介质、断褶带及构造变形显著区域，地层往往存在褶皱、断层等地质形态，导致介质的对称轴不再竖直，而是具有空间变化的倾角，此时实际介质更符合倾斜横向各向同性（tti）模型。若在这些区域仍采用vti介质假设进行建模与反演，不仅会导致介质参数表征失真，还可能引起梯度计算偏差，进而降低波场模拟和全波形反演的精度与稳定性。

3、针对tti介质，现有全波形反演方法主要从传统bond变换形式的tti介质弹性系数矩阵出发构建波动方程。该类方法虽然物理上严格，但在工程实现中存在两个突出问题：其一，tti介质弹性系数与速度、密度、thomsen参数及倾角之间呈现复杂的非线性关系，梯度推导通常需要反复进行链式求导，显著提高了计算代价与实现复杂度；其二，传统tti介质波动方程的算子一般不具对称结构，导致伴随状态方程在数学形式和自由表面边界条件上均与正演波动方程不一致，增加了tti介质全波形反演一体化实现的难度。具体而言，由于算子矩阵为非对称形式，直接转置波动方程算子得到的伴随状态方程中包含弹性系数与伴随状态变量相乘后再进行空间微分的复杂项，且伴随状态变量所对应的自由表面边界条件形式较为繁琐，进一步增加了精确数值实现的难度。值得注意的是，tti介质qp波方程与弹性波方程的数学结构基本一致，因此上述困难同样存在于tti介质声波全波形反演中。

4、针对上述方程非对称带来的问题，已有研究者在vti介质中提出了相应的解决方案。通过对波动方程中的正应力进行变量替换，能够实现vti介质波动方程的对称化，使得正演与伴随计算能够在统一的数值框架下实现，采用相同的离散格式和程序流程即可完成两者的数值求解，从而显著降低vti介质全波形反演的实现复杂度。然而，该类方法均建立在介质对称轴竖直（倾角为零）的假设基础上，仅适用于对称轴方向保持不变的特殊情形，无法直接推广至对称轴具有空间变化倾角的tti介质。

5、因此，如何实现正演、伴随和梯度计算的统一数值实现，并降低tti介质全波形反演的工程实现复杂度，成为亟待解决的技术难题。

技术实现思路

1、本技术实施例的目的在于提供一种基于对称波动方程的tti介质全波形反演方法及装置，实现正演、伴随和梯度计算的统一数值实现，并降低tti介质全波形反演的工程实现复杂度。

2、为达到上述目的，一方面，本技术实施例提供了一种基于对称波动方程的tti介质全波形反演方法，包括以下步骤：

3、获取地震观测数据、震源子波以及待反演区域的初始地下模型；

4、基于局部坐标系与全局坐标系之间的应力张量旋转关系，以vti介质弹性系数和倾角参数重构tti介质波动方程，以简化波动方程中的参数；

5、对重构后的tti介质波动方程进行对称化处理，通过拆分正应力项并引入辅助变量进行替换来构造正演算子为对称矩阵的对称形式tti介质波动方程，使所述对称形式tti介质波动方程的伴随状态方程在数学结构及边界条件上与正演波动方程保持一致；

6、基于伴随状态法，根据所述对称形式tti介质波动方程在统一数值框架下求解正传波场和伴随波场，并计算模型参数的梯度；

7、基于所述梯度对所述初始地下模型进行迭代更新，直至满足预设收敛条件，输出反演结果。

8、本技术实施例的方法中，所述基于局部坐标系与全局坐标系之间的应力张量旋转关系，以vti介质弹性系数和倾角参数重构tti介质波动方程，包括：

9、根据全局坐标系与局部坐标系下voigt格式应力张量之间的旋转关系，将tti介质弹性波方程等价转换为以vti介质弹性系数和倾角参数表达的形式；

10、

11、其中，为震源项，，、分别为震源项在水平方向、垂直方向的分量，为声压源项，，为对应的标量值，为全局坐标系下voigt格式应力张量，，、、分别为全局坐标系下voigt格式的水平正应力、垂直正应力、剪切应力，为局部坐标系下voigt格式的应力张量，，、、分别为局部坐标系下voigt格式的水平正应力、垂直正应力、剪切应力，为与tti介质对称轴相对于竖直方向的倾角相关的bond变换矩阵，，为tti介质对称轴相对于竖直方向的倾角，为介质密度，为单位矩阵，为时间，为空间微分算子矩阵，，、分别为水平方向坐标、垂直方向坐标，为质点振动速度矢量，，、分别为水平方向、垂直方向的质点振动速度，为vti介质弹性系数矩阵，，、、、为vti介质弹性系数，为时间，上标表示矩阵转置。

12、本技术实施例的方法中，所述对重构后的tti介质波动方程进行对称化处理，并通过拆分正应力项并引入辅助变量构造正演算子为对称矩阵的对称形式tti介质波动方程，包括：

13、对重构后的tti介质弹性波方程中与正应力相关的方程项进行拆分：

14、

15、其中，为局部坐标系下的水平方向正应力，为拆分后的第一部分，为拆分后的第二部分，、分别为质点振动速度在水平和垂直方向的分量，为tti介质对称轴相对于竖直方向的倾角，、、为vti介质弹性系数，为声压震源项，为时间，、分别为水平方向坐标、垂直方向坐标；

16、构建满足的辅助变量和，并将所述辅助变量和对拆分后的方程进行变量替换，构造得到对称形式的tti介质弹性波方程。

17、本技术实施例的方法中，所述对称形式tti介质弹性波方程的矩阵形式为：

18、

19、其中，为介质密度，为2×2单位矩阵，为空间微分算子矩阵，为与tti介质对称轴相对于竖直方向的倾角相关的bond变换矩阵，，为选择矩阵，，为质点振动速度矢量，，为变量替换后的应力相关变量向量，，为震源项，，、分别为震源项在水平方向、垂直方向的分量，为变量替换后的声压震源项，，为由vti介质弹性系数组合的对角矩阵，，，和，上标表示矩阵转置。

20、本技术实施例的方法中，所述对称形式tti介质弹性波方程的伴随状态方程的矩阵形式为：

21、

22、其中，、分别为伴随波场的速度分量、应力相关变量向量，、分别为速度伴随源项、应力伴随源项；所述伴随状态方程与正演波动方程具有相同的算子矩阵和边界条件，使得正演与伴随在统一数值框架下采用相同的求解器实现。

23、本技术实施例的方法中，在声学近似条件下所述对称形式的tti介质弹性波方程退化为对称形式的tti介质qp波方程。

24、本技术实施例的方法中，所述基于伴随状态法，根据所述对称形式tti介质波动方程在统一数值框架下求解正传波场和伴随波场，并计算模型参数的梯度，包括：

25、根据所述对称形式的tti介质波动方程进行正演模拟，得到正传波场和预测数据；

26、根据目标函数和所述预测数据计算数据残差；

27、根据所述对称形式的tti介质波动方程反传所述数据残差，得到伴随波场；

28、根据梯度公式计算目标函数关于各模型参数的梯度；

29、其中，为目标函数，为模型参数向量的第个模型参数，为观测数据的最大记录时间，为时间，为介质密度，为由vti介质弹性系数组合的对角矩阵，为正传波场的速度分量，为伴随波场的速度分量，为正传波场的应力分量，为伴随波场的应力分量，、分别为、的转置矩阵；所述模型参数包括纵波速度、横波速度、介质密度、thomsen各向异性参数和。

30、本技术实施例的方法中，所述目标函数为最小二乘目标函数：

31、

32、其中，为以为变量的最小二乘目标函数，为空间区域，为预测数据，为观测数据，为检波点坐标，为单位脉冲函数，为空间坐标；所述模型参数向量包括纵波速度、横波速度、介质密度、thomsen各向异性参数和。

33、另一方面，本技术实施例还提供了一种基于对称波动方程的tti介质全波形反演装置，包括：

34、数据获取单元，用于获取地震观测数据、震源子波以及待反演区域的初始地下模型；

35、方程重构单元，用于基于局部坐标系与全局坐标系之间的应力张量旋转关系，以vti介质弹性系数和倾角参数重构tti介质波动方程，以简化波动方程中的参数；

36、方程对称化单元，用于对重构后的tti介质波动方程进行对称化处理，通过拆分正应力项并引入辅助变量构造正演算子为对称矩阵的对称形式tti介质波动方程，使所述对称形式tti介质波动方程的伴随状态方程在数学结构及边界条件上与正演波动方程保持一致；

37、梯度计算单元，用于基于伴随状态法，根据所述对称形式tti介质波动方程在统一数值框架下求解正传波场和伴随波场，并计算模型参数的梯度；

38、迭代更新单元，用于基于所述梯度对所述初始地下模型进行迭代更新，直至满足预设收敛条件，输出反演结果。

39、另一方面，本技术实施例还提供了一种计算机设备，包括存储器、处理器、以及存储在所述存储器上的计算机程序，所述计算机程序被所述处理器运行时，执行上述方法的指令。

40、由以上本技术实施例提供的技术方案可见，本技术实施例中，在重构tti介质波动方程后，可通过拆分正应力项并引入辅助变量构造正演算子为对称矩阵的对称形式tti介质波动方程，由于正演算子矩阵为对称矩阵，意味着该矩阵的转置等于自身；根据伴随状态法，伴随方程的算子矩阵在数学上即为正演算子矩阵的转置，二者自然完全相同；而且对称形式tti介质波动方程的伴随状态方程在数学结构及边界条件上与正演波动方程保持一致，使正演与伴随在数值离散格式和边界处理方式上也完全一致。因此，正演与伴随求解的区别在于震源项不同（正演为物理震源，伴随为数据残差），在计算机实现中完全可以采用同一套数值求解程序，仅需将输入的源项从物理震源替换为数据残差即可完成伴随波场的求解，无需为伴随方程单独编写求解代码，从而显著降低了实现复杂度和工程实现成本。