CN116088048A - 包含不确定性分析的各向异性介质多参数全波形反演方法 - Google Patents

Abstract

本发明为一种包含不确定性分析的各向异性介质多参数全波形反演方法，包括：建立具有倾斜对称轴的横向各向异性介质中各独立弹性模量的初始模型；采用弹性模量的初始模型进行正演数值模拟得到合成地震数据集；计算观测地震数据与模拟地震数据之间的数据残差；基于贝叶斯概率体系推导最大后验解；由低频到高频，利用迭代扩展卡尔曼滤波法进行下一频率模型参数以及后验协方差矩阵的迭代更新，前一个频率反演结果作为下一个频率反演的初始模型，直到计算频率值达到所设定的最大频率值，且最终的数据残差满足精度要求。本发明解决了现有弹性各向异性介质多参数反演与不确定性定量估计问题。

Description

包含不确定性分析的各向异性介质多参数全波形反演方法

技术领域

本发明属于勘探地球物理学领域，具体涉及到一种包含不确定性分析的各向异性介质多参数全波形反演方法。

背景技术

随着地球物理学科的进步以及勘探水平的不断提升，目前地震储层预测偏向于解决更复杂的地球物理难题。在常见的沉积地层以及复杂岩层构造中，地震波传播速度随传播方向变化而发生改变，因此如果充分考虑地下介质的各向异性特征，可以获得更为精确的反演结果。但是随着参数数量的增加，参数反演的难度增加且不同参数之间的串扰效应会对反演结果造成很大的影响。根据Tarantola,A.,2005,Inverse problem theory andmethods for model parameter estimation:SIAM,89.，从地震反射数据中确定地球物性参数的过程，可以被理解为一个逆散射问题，而散射理论实质上是一种针对扰动介质的分析。地震散射理论在处理地震波长较小的局部非均匀体，如含孔洞、裂缝等小尺度介质以及强散射体时具有很大优势。将地震散射理论用于全波形反演可以更为精确地刻画复杂各向异性介质的物性参数信息。

全波形反演方法充分利用了地震波的运动学和动力学特征来获取地下模型参数信息，具有复杂构造成像精度高、物性参数反演效果好等优点(Lailly,P.,1983,Theseismic inverse problem as a sequence of before stack migrations:Conferenceon Inverse Scattering:Theory and Application,Society for Industrial andApplied Mathematics,Expanded Abstracts,206–220.)。而实际工业生产中，初始模型的欠准、低频信息的缺失和大偏移距观测数据的限制、照明度不足以及多参数串扰效应等问题的存在，会给全波形反演结果带来很大的不确定性。根据Eikrem,K.S.,G.Naevdal,andM.Jakobsen,2019,Iterated Extended Kalman Filter Method for Time-lapse SeismicFull Waveform Inversion:Geophysical Prospecting,67,no.2,379-394,doi:10.1111/1365-2478.12730.，与传统反演方法相比，贝叶斯反演方法的优点在于可以量化反演结果的不确定性，使先验信息得到有效地利用。散射积分方法将正演结果以李普曼-施温革方程的形式给出，取代了传统有限差分正演方法，避免了其可能会引起的网格差分误差，更有利于实现地下介质精细构造成像(Jakobsen,M.,E.Ivan,I.Psencik,and B.Ursin,2020,Transition operator approach to seismic full-waveform inversion in arbitraryanisotropic elastic media:Communications in Computational Physics,27,no.1,1–31,doi:10.4208/cicp.OA-2018-0197.)。基于散射理论的贝叶斯全波形反演方法可以更可靠地识别地下介质参数并提供不确定性信息，适用于解决各向异性介质多参数反演问题，对于解释地质构造和更新储层模型都具有非常重要的意义。

发明内容

本发明所要解决的技术问题在于提供一种包含不确定性分析的各向异性介质多参数全波形反演方法，在地震数据缺失低频信息、信噪比较低的情况下反演得到分辨率较高的地下各向异性介质多参数信息，并定量估计全波形反演结果的不确定性。

本发明是这样实现的，

一种包含不确定性分析的各向异性介质多参数全波形反演方法，该方法包括：

步骤1：输入观测地震数据，提取地震数据的走时信息，进行层析成像，得到层析速度模型；

步骤2：利用傅里叶变换，将时间域地震数据转化到频率域，绘制其频谱，确定有效信息较为丰富的频带范围，得到分频地震数据d_obs；

步骤3：使用Gardner公式求得密度模型，结合层析速度模型建立各向异性VTI介质中各独立弹性模量模型，利用Bond变换得到各向异性TTI介质的各独立弹性模量，作为反演的初始模型；

步骤4：基于散射理论，求解各向异性TTI介质频率域弹性动力学波动方程进行正演数值模拟，利用三阶、四阶格林函数构建Lippmann-Swinger散射积分方程形式的位移-应变耦合积分方程；

步骤5：利用包含离散小波变换预处理算子的预处理广义最小残差法，在克雷洛夫子空间求解应变散射积分方程，进而求解粒子位移场u，得到模拟地震数据集d_cal；

步骤6：计算观测地震数据与模拟地震数据之间的数据残差，基于模型与数据之间的映射关系，根据数据残差得到敏感性核；

步骤7：根据分频地震数据d_obs，基于贝叶斯概率体系构建各向异性TTI介质全波形反演最大后验解以及后验协方差；

步骤8：利用迭代扩展卡尔曼滤波法迭代求取贝叶斯全波形反演模型参数最大后验解以及后验协方差更新量，若满足迭代停止条件，则输出这一频率的模型参数以及后验协方差，作为下一频率的输入。

步骤9：重复步骤4-步骤8进行下一频率模型参数以及后验协方差矩阵的迭代更新，直到计算频率值达到所设定的最大频率值，且最终的数据残差满足精度要求，此时所得到的模型参数即为各向异性TTI介质最终的反演结果。

进一步地，所述步骤3中，使用Gardner公式求得密度模型，结合层析速度模型建立各向异性VTI介质中各独立弹性模量模型，利用Bond变换得到各向异性TTI介质的各独立弹性模量，作为反演的初始模型，具体包括：

用初始速度模型中的P波速度v_P、和S波速度v_S，密度ρ以及汤姆森各向异性参数模型的Thomsen参数ε和δ表征各向异性VTI介质弹性模量，构成VTI弹性模量矩阵:

假设z轴与各向异性介质对称轴之间的夹角为θ，将各向异性VTI介质弹性模量矩阵相对于对称轴的方向进行Bond变换

得到一组新的刚度张量，即为公式(2)表示的各向异性TTI介质弹性模量矩阵，作为全波形反演各向异性TT I介质中各独立弹性模量的初始模型:

进一步地，所述步骤4中，求解各向异性TTI介质频率域弹性动力学波动方程进行正演数值模拟，利用三阶、四阶格林函数构建Lippmann-Swinger散射积分方程形式的位移-应变耦合积分方程，具体包括：

采用公式(3)的频率域弹性动力学波动方程进行正演模拟，求取粒子位移场：

式中，ρ和ω分别代表密度和角频率，u(x)为粒子位移量，S(x)为震源。

为点x处的应变场，I表示单位矩阵；将弹性常数C(x)分解为背景介质C⁽⁰⁾(x)和扰动介质δC(x)两个分量，引入背景格林函数

求解公式(3)，再利用弹性模量的局部积分对称性与弹性动力学粒子位移-应变关系，得到位移-应变场的积分方程为公式(4)和公式(5):

其中，三阶与四阶格林函数张量可以写为

进一步地，根据公式(4)和公式(5)将位移-应变方程写为一对积分耦合算子的形式:

其中，其中，

分别为三阶、四阶格林函数张量，u⁽⁰⁾、ε⁽⁰⁾为背景介质位移量和应变量，V(x₁,x₂)＝δC(x₁)(x₁-x₂)为散射势算子，ε为应变场，x₁和x₂为模型区域中相近但不重合的两点；δC(x₁)表示x₁点处迭代更新前后的弹性模量差值，通过求解公式(6)得到模拟地震数据集d_cal。

进一步地，

采用迭代求解方法求解公式(7)的应变场，公式(7)写为:

令

为

并将矩阵ε⁽⁰⁾分成n个列矩阵块，另

将公式8的求解转化为求解匹配ε⁽⁰⁾和

之间最小化的近似解x^(m)，其数学表达式为：

通过Arnoldi循环构造矩阵

的正交克雷洛夫子空间

对任意

设x＝x⁽⁰⁾+V_my，残差r的计算公式为:

其中，V_m＝{v₁,v₂,...,v_m}是

的标准正交基，满足AV_m＝V_m+1H_m+1,m，将公式(9)中的最小二乘问题简化为:

其中，e₁指单位向量，β指残差的二次范数。

进一步地，在达到迭代步数k时,令x⁽⁰⁾＝x^(k)，重新开始新的一轮迭代直至收敛，并引入

作为预处理算子，使得：

进一步地，所述步骤6，具体包括以下步骤：

步骤6.1：求解第i次迭代的数据残差

用于衡量目前模型的精度是否满足迭代停止需求，其中ΔC⁽ⁱ⁺¹⁾(x)为当前和前一次迭代之间的弹性模量之差:

步骤6.2：引入弹性模量扰动量的计算公式

将式(14)改写为：

其中，m^(p)＝(C(x)-C⁽⁰⁾(x))/C⁽⁰⁾(x)为模型参数，p为弹性模量个数；B^(p)(x)为与模型参数相关的结构矩阵，p表示21个弹性张量之一；

步骤6.3：基于模型与数据之间的关系δu⁽ⁱ⁾＝J⁽ⁱ⁾·δm⁽ⁱ⁾，计算敏感性核为:

步骤6.4：将敏感性核和扰动更新量写成分块矩阵的形式:

J⁽ⁱ⁾＝[J^(i,1),J^(i,2),...,J^(i,21)]^T   (17)

δm⁽ⁱ⁾＝[δm^(i,1),δm^(i,2),...,δm^(i,21)]^T   (18)。

进一步地，所述步骤7，具体包括：

通过公式(19)和公式(20)计算基于贝叶斯概率体系构建模型参数最大后验解以及后验协方差：

其中，C_D、C_M分别为数据协方差和模型协方差，m_pri为先验模型，d_k表示分频地震数据集d_obs的第k组观测数据集，d_cal为计算数据，m_pri为先验模型参数。

进一步地，所述步骤8，具体包括：

利用迭代扩展卡尔曼滤波法求解贝叶斯概率体系构建模型参数最大后验解以及后验协方差的迭代解为：

其中λ为动态因子，I为单位矩阵。

本发明与现有技术相比，有益效果在于：

本发明用于处理常规频率域地震数据，引入贝叶斯概率框架，建立基于概率体系的最大后验解，代替传统以观测数据和模拟数据最小二乘作为目标函数的全波形反演方法，同时反演各向异性TTI介质的各独立弹性参数，在得到反演结果的同时可以获得不确定性信息。模型与协方差信息更新策略采用“低频到高频，后验即先验”的多尺度反演策略。基于散射理论，充分考虑地下介质的多散射性，实现地下结构的精细成像。采用散射积分方法对模型进行正演数值模拟，通过使用格林函数计算多散射敏感性核，得到近似海森信息，进而构建协方差矩阵。在背景模型复杂程度较低的情况下，散射积分方程自动满足索末菲辐射条件，无需设置吸收边界，减小了网格离散误差给反演结果带来的不确定性。散射积分方法使得反演只对异常体区域进行离散。在此基础上，结合线性方程迭代求解算子，一定程度上降低了大尺度模型反演结果的计算成本。提高反演结果质量的同时，提高了计算效率，加大了其工业应用的有效性。本发明即使在缺失精确的初始模型和低频地震记录的情况下，其反演结果也具有良好的鲁棒性，适合于陆地和海洋等多种领域的研究和工业生产应用。

附图说明

图1为传统全波形反演算法的流程图；

图2为本发明提出的一种包含不确定性分析的各向异性介质多参数全波形反演方法的流程图；

图3为本发明所提供的模拟测试中的真实参数模型，其中(a)-(f)分别显示了P波、S波速度、密度、Thomsen参数以及模型各向异性倾斜角度；

图4为本发明所提供的模拟测试中的真实弹性模量模型，其中(a)-(e)分别显示了不同的弹性模量C₁₁，C₃₃，C₅₅，C₆₆，C₁₃；

图5为本发明所提供的模拟测试中的初始弹性模量模型，其中(a)-(e)分别显示了不同的弹性模量C₁₁，C₃₃，C₅₅，C₆₆，C₁₃；

图6为本发明所提供的模拟测试中的弹性模量模型的最终反演结果，其中(a)-(e)分别显示了不同的弹性模量C₁₁，C₃₃，C₅₅，C₆₆，C₁₃；

图7为本发明所提供的模拟测试中的弹性模量模型反演结果与真实模型之间的误差，其中(a)-(e)分别显示了不同的弹性模量C₁₁，C₃₃，C₅₅，C₆₆，C₁₃；

图8为本发明所提供的模拟测试中的弹性模量模型反演结果的后验标准差，其中(a)-(e)分别显示了不同的弹性模量C₁₁，C₃₃，C₅₅，C₆₆，C₁₃；

图9为本发明所提供的模拟测试中的任一弹性模量的先验协方差矩阵；

图10为本发明所提供的模拟测试中的所有弹性模量的先验协方差矩阵，从左上到右下对角线上的块矩阵分别为C₁₁，C₃₃，C₅₅，C₆₆，C₁₃；

图11为本发明所提供的模拟测试中的所有弹性模量的后验协方差矩阵，从左上到右下对角线上的块矩阵分别为C₁₁，C₃₃，C₅₅，C₆₆，C₁₃。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明了，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

传统全波形反演方法参见图1所示，本发明提出的一种包含不确定性分析的各向异性介质多参数全波形反演方法的流程如图2所示，带入理论测试实例，实施过程包括如下步骤：

步骤1：输入观测地震数据，提取地震数据的走时信息，进行层析成像，得到层析速度模型；

步骤2：利用傅里叶变换，将时间域地震数据转化到频率域，绘制其频谱，确定有效信息较为丰富的频带范围，得到分频地震数据d_obs；

步骤3：使用Gardner公式求得密度模型，结合层析速度模型建立各向异性VTI介质中各独立弹性模量模型，利用Bond变换得到各向异性TTI介质的各独立弹性模量，作为反演的初始模型；

步骤4：基于散射理论，求解各向异性TTI介质频率域弹性动力学波动方程进行正演数值模拟，利用三阶、四阶格林函数构建Lippmann-Swinger散射积分方程形式的位移-应变耦合积分方程；

步骤5：利用包含离散小波变换预处理算子的预处理广义最小残差法，在克雷洛夫子空间求解应变散射积分方程，进而求解粒子位移场u，得到模拟地震数据集d_cal；

步骤6：计算观测地震数据与模拟地震数据之间的数据残差，基于模型与数据之间的映射关系，根据数据残差得到敏感性核；

步骤7：根据分频地震数据d_obs，基于贝叶斯概率体系构建各向异性TTI介质全波形反演最大后验解以及后验协方差；

步骤8：利用迭代扩展卡尔曼滤波法迭代求取贝叶斯全波形反演模型参数最大后验解以及后验协方差更新量，若满足迭代停止条件，则输出这一频率的模型参数以及后验协方差，作为下一频率的输入。

步骤9：重复步骤4-步骤8进行下一频率模型参数以及后验协方差矩阵的迭代更新，直到计算频率值达到所设定的最大频率值，且最终的数据残差满足精度要求，此时所得到的模型参数即为最终的各向异性TTI介质反演结果。

步骤3中，具体包括以下步骤：

步骤1.1：用P波速度v_P、S波速度v_S，密度ρ,以及Thomsen参数ε和δ表征各向异性VTI介质弹性模量:

步骤1.2：各向异性VTI介质弹性模量矩阵构成为:

步骤1.3：假设z轴与各向异性介质对称轴之间的夹角为θ，将各向异性VTI介质弹性模量矩阵相对于对称轴的方向进行Bond变换，得到一组新的刚度张量，即为各向异性TTI介质弹性模量矩阵:

倾斜角矩阵的计算公式为：

步骤1.5：根据各向异性TTI介质弹性模量矩阵得到各向异性TTI介质刚度矩阵中各独立弹性模量，作为全波形反演的各向异性模型参数:

步骤4中，基于散射理论对各向异性TTI介质初始模型进行正演数值模拟，构建Lippmann-Swinger散射积分方程形式的位移-应变耦合积分方程，并利用基于克雷洛夫子空间的广义最小残差法，加速求解各向异性TTI介质的应变量，进而求解粒子位移场，获得模拟地震数据集，包括以下步骤：

引入背景格林函数

求解频率域弹性动力学波动方程，得到粒子位移场表达式:

式中，x′为x附近某一点，u⁽⁰⁾(x)表示背景介质的波场；

利用弹性模量的局部积分对称性，将式(4)改写为:

其中

为三阶格林函数张量；

由弹性动力学粒子位移-应变关系，得到应变场的积分方程：

其中

为四阶格林函数张量。

根据公式(4)和公式(5)将位移-应变方程写为一对积分耦合算子的形式:

其中，V(x₁,x₂)＝δC(x₁)(x₁-x₂)为散射势算子；

引入迭代求解方法来解决线性系统求解问题求解应变场:

令

为

并将矩阵ε⁽⁰⁾分成n个列矩阵块，则应变量求解方程改写为:

求解匹配ε⁽⁰⁾和

之间最小化的近似解x^(m)为：

通过Arnoldi循环构造矩阵

的正交克雷洛夫子空间κ_m:

对任意x＝x⁽⁰⁾+κ_m，设x＝x⁽⁰⁾+V_my，残差r的计算公式为公式(11):

其中，V_m＝{v₁,v₂,...,v_m}是

的标准正交基，满足AV_m＝V_m+1H_m+1,m；

步骤2.11：最小二乘问题化简为:

其中，e₁指的是单位向量，β指的是残差的二次范数。

使得

作为预处理算子，在达到重启迭代步数k时,令x⁽⁰⁾＝x^(k)，重新开始新的一轮迭代直至收敛：

步骤6，具体包括以下步骤：

步骤6.1：求解第i次迭代的数据残差

用于衡量目前模型的精度是否满足迭代停止需求，其中ΔC⁽ⁱ⁺¹⁾(x)为当前和前一次迭代之间的弹性模量之差:

步骤6.2：引入弹性模量扰动量的计算公式

将式(14)改写为：

其中，m^(p)＝C⁽⁰⁾(x)(I+ΔC(x))为模型参数，B^(p)(x)为与模型参数相关的结构矩阵；

步骤6.3：基于模型与数据之间的关系，得到敏感性核为:

步骤6.4：将敏感性核和扰动更新量写成分块矩阵的形式:

数据残差用模型扰动量以及敏感性核的乘积表示为：

δu⁽ⁱ⁾＝J⁽ⁱ⁾·δm⁽ⁱ⁾

步骤7具体包括以下步骤：

步骤7.1：基于贝叶斯概率体系，求得全波形反演模型参数最大后验解以及后验协方差的表达式为:

其中，C_D、C_M分别为数据协方差和模型协方差，m_pri为先验模型，d_k表示分频地震数据集d_obs的第k组观测数据集，d_cal为计算数据，m_pri为先验模型参数。

利用迭代扩展卡尔曼滤波法迭代求解模型参数，得到最大后验解的以及后验协方差迭代解为:

实施例

建立各向异性TTI介质弹性模量反演初始参数模型。

根据图3所示的真实参数信息，设置倾斜角θ，构建各向异性TTI介质真实弹性模量模型，如图4所示，计算公式如下：

其中

根据已有的地震剖面结果和测井资料、地质岩性物理资料等信息，建立所测试的各向异性TTI介质的弹性模量初始参数模型，此例中初始模型为类似光滑梯度模型，如图5所示。

根据真实参数模型大小设置适当的地震观测系统，并加入噪声干扰，使得实验更贴向实际情况。

设置模型网格数量为40×20，格间距25米，模型大小为1.0×0.5km。采用主频为7.5Hz的雷克子波为震源子波，采用40个震源和40个检波器，所有震源与检波器都均匀置于模型的顶部。背景介质模型为高斯平滑后的均匀各向同性介质模型。

设置信噪比为4，信噪比的计算公式为：

其中d为频率数据，w为角频率，E为期望值。

设置模型反演频率范围为4Hz到20Hz。

反演过程由低频到高频依次进行，前一个频率得到的模型将作为下一个频率更新时所用的初始模型，按照设置频率范围由低到高，反演结果如图6所示。最大后验解和真实模型之间的误差如图7所示，标准差如图8所示。模型中任一弹性模量和所有弹性模量的先验协方差矩阵分别如图9和图10所示，反演结果的后验协方差矩阵如图11所示。

从反演结果可以得到，即使初始速度模型与真实模型相差较大，且信噪比较低的情况下，本发明仍然可以很好地反演出各向异性TTI介质的五个独立的弹性模量。反演结果虽然在模型深部以及边界处仍然存在一定误差，但相比较初始模型有了大幅提升，取得了较为理想的效果。并且由后验协方差以及反演结果的分辨率可以观察到各个参数之间的串扰关系。参数C₃₃和C₅₅的反演结果分辨率最高，其次是C₁₁和C₆₆，最后是C₁₃。因为C₃₃和C₅₅可以被认为只与P波和S波的垂向速度有关，因此这两个系数在更新过程中比其他系数独立性更强，这主要是因为地表数据是由垂直分量地震记录组成的。在5个弹性常数中，只有C₆₆与SH波有关，其他四个弹性参数C₁₁,C₅₅,C₃₃,C₁₃与P波和SV波的相速度和极化有关。然而，C₁₁和C₃₃也存在一定的串扰关系，一定程度上两者的反演结果会互相产生影响。这是因为C₁₁和C₃₃接近所谓的P波各向异性，这与P波在水平方向和垂直方向上的速度差有关。此外，由于SH波的各向异性强度影响相较于其他4个弹性模量偏弱。这是因为C₁₃不仅会受到C₅₅的影响，还会受到C₃₃的影响。这三个系数与垂直入射时P波相速度的二阶导数密切相关。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种包含不确定性分析的各向异性介质多参数全波形反演方法，其特征在于，该方法包括：

步骤1：输入观测地震数据，提取地震数据的走时信息，进行层析成像，得到层析速度模型；

步骤2：利用傅里叶变换，将时间域地震数据转化到频率域，绘制其频谱，确定有效信息较为丰富的频带范围，得到分频地震数据d_obs；

步骤3：使用Gardner公式求得密度模型，结合层析速度模型建立各向异性VTI介质中各独立弹性模量模型，利用Bond变换得到各向异性TTI介质的各独立弹性模量，作为反演的初始模型；

步骤4：基于散射理论，求解各向异性TTI介质频率域弹性动力学波动方程进行正演数值模拟，利用三阶、四阶格林函数构建Lippmann-Swinger散射积分方程形式的位移-应变耦合积分方程；

步骤5：利用包含离散小波变换预处理算子的预处理广义最小残差法，在克雷洛夫子空间求解应变散射积分方程，进而求解粒子位移场u，得到模拟地震数据集d_cal；

步骤6：计算观测地震数据与模拟地震数据之间的数据残差，基于模型与数据之间的映射关系，根据数据残差得到敏感性核；

步骤7：根据分频地震数据d_obs，基于贝叶斯概率体系构建各向异性TTI介质全波形反演最大后验解以及后验协方差；

步骤8：利用迭代扩展卡尔曼滤波法迭代求取贝叶斯全波形反演模型参数最大后验解以及后验协方差更新量，若满足迭代停止条件，则输出这一频率的模型参数以及后验协方差，作为下一频率的输入；

步骤9：重复步骤4-步骤8进行下一频率模型参数以及后验协方差矩阵信息的迭代更新，直到计算频率值达到所设定的最大频率值，且最终的数据残差满足精度要求，此时所得到的模型参数即为最终的各向异性TTI介质反演结果。

2.根据权利要求1所述的一种包含不确定性分析的各向异性介质多参数全波形反演方法，其特征在于，所述步骤3中，使用Gardner公式求得密度模型，结合层析速度模型建立各向异性VTI介质中各独立弹性模量模型，利用Bond变换得到各向异性TTI介质的各独立弹性模量，作为反演的初始模型，具体包括：

用初始速度模型中的P波速度v_P、和S波速度v_S，密度ρ以及汤姆森各向异性参数模型的Thomsen参数ε和δ表征各向异性VTI介质弹性模量，构成VTI弹性模量矩阵:

假设z轴与各向异性介质对称轴之间的夹角为θ，将各向异性VTI介质弹性模量矩阵相对于对称轴的方向进行Bond变换

得到一组新的刚度张量，即为公式(2)表示的各向异性TTI介质弹性模量矩阵，作为全波形反演各向异性TTI介质中各独立弹性模量的初始模型:

3.根据权利要求1所述的一种包含不确定性分析的各向异性介质多参数全波形反演方法，其特征在于，所述步骤4中，求解各向异性TTI介质频率域弹性动力学波动方程进行正演数值模拟，利用三阶、四阶格林函数构建Lippmann-Swinger散射积分方程形式的位移-应变耦合积分方程，具体包括：

采用公式(3)的频率域弹性动力学波动方程进行正演模拟，求取粒子位移场：

式中，ρ和ω分别代表密度和角频率，u(x)为粒子位移量，S(x)为震源，

为点x处的应变场，I表示单位矩阵；将弹性常数C(x)分解为背景介质C⁽⁰⁾(x)和扰动介质δC(x)两个分量，引入背景格林函数

求解公式(3)，再利用弹性模量的局部积分对称性与弹性动力学粒子位移-应变关系，得到位移-应变场的积分方程为公式(4)和公式(5):

其中，三阶与四阶格林函数张量写为

4.根据权利要求3所述的一种包含不确定性分析的各向异性介质多参数全波形反演方法，其特征在于，根据公式(4)和公式(5)将位移-应变方程写为一对积分耦合算子的形式:

其中，其中，

分别为三阶、四阶格林函数张量，u⁽⁰⁾、ε⁽⁰⁾为背景介质位移量和应变量，V(x₁,x₂)＝δC(x₁)(x₁-x₂)为散射势算子，ε为应变场，x₁和x₂为模型区域中相近但不重合的两点；δC(x₁)表示x₁点处迭代更新前后的弹性模量差值，通过求解公式(6)得到模拟地震数据集d_cal。

5.根据权利要求4所述的一种包含不确定性分析的各向异性介质多参数全波形反演方法，其特征在于，

采用迭代求解方法求解公式(7)的应变场，公式(7)写为:

令

为

并将矩阵ε⁽⁰⁾分成n个列矩阵块，另

将公式8的求解转化为求解匹配ε⁽⁰⁾和

之间最小化的近似解x^(m)，其数学表达式为：

通过Arnoldi循环构造矩阵

的正交克雷洛夫子空间κ_m:

对任意x＝x⁽⁰⁾+κ_m，设x＝x⁽⁰⁾+V_my，残差r的计算公式为:

其中，V_m＝{v₁,v₂,...,v_m}是

的标准正交基，满足AV_m＝V_m+1H_m+1,m，将公式(9)中的最小二乘问题简化为:

其中，e₁指单位向量，β指残差的二次范数。

6.根据权利要求5所述的一种包含不确定性分析的各向异性介质多参数全波形反演方法，其特征在于，在达到迭代步数k时,令x⁽⁰⁾＝x^(k)，重新开始新的一轮迭代直至收敛，并引入

作为预处理算子，使得：

7.根据权利要求1所述的一种包含不确定性分析的各向异性介质多参数全波形反演方法，其特征在于，所述步骤6，具体包括以下步骤：

步骤6.1：求解第i次迭代的数据残差

用于衡量目前模型的精度是否满足迭代停止需求，其中ΔC⁽ⁱ⁺¹⁾(x)为当前和前一次迭代之间的弹性模量之差:

步骤6.2：引入弹性模量扰动量的计算公式

将式(14)改写为：

其中，m^(p)＝(C(x)-C⁽⁰⁾(x))/C⁽⁰⁾(x)为模型参数，p为弹性模量个数；B^(p)(x)为与模型参数相关的结构矩阵，p表示21个弹性张量之一，B^(p)(x)为与模型参数相关的结构矩阵；

步骤6.3：基于模型与数据之间的关系δu⁽ⁱ⁾＝J⁽ⁱ⁾·δm⁽ⁱ⁾，计算敏感性核为:

步骤6.4：将敏感性核和扰动更新量写成分块矩阵的形式:

J⁽ⁱ⁾＝[J^(i,1),J^(i,2),...,J^(i,21)]^T          (17)

δm⁽ⁱ⁾＝[δm^(i,1),δm^(i,2),...,δm^(i,21)]^T         (18)。

8.根据权利要求7所述的一种包含不确定性分析的各向异性介质多参数全波形反演方法，其特征在于，所述步骤7，具体包括：

通过公式(19)和公式(20)计算基于贝叶斯概率体系构建模型参数最大后验解以及后验协方差：

其中，C_D、C_M分别为数据协方差和模型协方差，m_pri为先验模型，d_k表示分频地震数据集d_obs的第k组观测数据集，d_cal为计算数据，m_pri为先验模型参数。

9.根据权利要求8所述的一种包含不确定性分析的各向异性介质多参数全波形反演方法，其特征在于，所述步骤8，具体包括：

利用迭代扩展卡尔曼滤波法求解贝叶斯概率体系构建模型参数最大后验解以及后验协方差的迭代解为：

其中λ为动态因子，I为单位矩阵。

Images