CN108037531B - 一种基于广义全变分正则化的地震反演方法及系统 - Google Patents

Abstract

本发明公开了一种基于广义全变分正则化的地震反演方法及系统，旨在解决现有地震反演技术采用全变分正则化会在地层内部产生阶梯效应的问题；本发明中方法利用广义全变分除了使用待反演参数的一阶偏导信息外，还使用了二阶甚至更高阶信息这一特点，能够在清晰勾勒地层分界面的同时减弱地层内部的阶梯效应，同时本发明发现了可以利用子波矩阵、差分矩阵和待反演参数的自然对数三者的乘积变换为子波矩阵对应的卷积核、差分矩阵对应的卷积核以及待反演参数的自然对数三者的卷积的现象，将时域的卷积运算转化为频域的点乘运算，大大提高了反演速度，且本发明无需先反演反射系数，在频域直接反演得到待反演参数；本申请适用于地震勘探技术领域。

Description

一种基于广义全变分正则化的地震反演方法及系统

技术领域

本发明涉及属于地震勘探技术领域，涉及一种基于广义全变分正则化的地震反演方法及系统。

背景技术

波阻抗、纵横波速度、密度、孔隙度等地球物理参数与岩石性质密切相关。因此，地球物理参数反演一直是勘探地震学领域的研究热点。近年来，研究人员发现利用全变分(Total Variation，TV)正则化方法可以得到块状的反演结果，因此该方法在地震反演方面受到了重视。此外，研究人员还发现，与传统的逐道反演方法相比，多道同时反演方法考虑了地震数据多道相邻数据之间的联系，能充分利用地震数据的横向连续性信息，提高反演结果的横向连续性和抗噪声性能，从而揭示更多地质结构特征。由于叠前、叠后参数反演方法具有类似的形式，为了简便起见，下面以叠后地震波阻抗反演为例，介绍近年来基于全变分正则化的地震反演方法的特点。

Gholami发表在GEOPHYSICS期刊2015年第5期第R217-R224页，题为“Nonlinearmultichannel impedance inversion by total-variation regularization”的文献，公开了一种基于全变分(Total Variation，TV)正则化的非线性多道地震波阻抗反演方法。该方法结合了多道同时反演在利用地震数据时空信息方面的特点和全变分正则化能获得块状波阻抗的优势。在求解方面，该方法采用DCT变换使得计算速度得以提升。但是，该方法需要先反演反射系数，再由反射系数与波阻抗之间的积分关系，利用积分运算求得波阻抗。反射系数反演的结果是含有偏差的，而要由反射系数得到波阻抗还需要知道t＝0时刻的波阻抗值，这个值是很难精确知道的。从反射系数得到波阻抗的过程是一个积分过程，由于反射系数反演结果的不精确、t＝0时刻波阻抗的值也不能精确知道，所以，这个积分过程易导致累积误差。

另一篇文献是申请人于2017年3月发表在Journal of Geophysics andEngineering第3期520-522页的，题为“Seismic acoustic impedance inversion withmulti-parameter regularization”的文献。该文献中，公布了一种基于多参数正则化的地震波阻抗反演方法。这种方法综合了对反射系数的L1范数稀疏约束、对波阻抗的全变分约束、初始模型约束和多道同时反演等手段。实验结果显示，该方法与基于L1范数正则化多道同时反演方法、基于全变分正则化多道同时反演方法相比，具有更优的抗噪声性能、更小的反演误差、更能清晰刻画地层的分界面。同时，作者的研究也指出，与基于L1范数正则化多道同时反演方法相比，基于全变分正则化多道同时反演的反演效果有明显提升。与基于全变分正则化多道同时反演方法相比，同时具有L1范数正则化和全变分正则化的多道同时反演对反演效果只是略有提升，但是却带来了更复杂的正则化参数选择问题。作者还指出，这类多道同时反演方法的计算效率还有待提高。

近年来，图像与信号处理领域的研究表明全变分正则化虽然能较好的保护图像的边缘信息，但是在内部的光滑区域却会产生明显的阶梯效应。对应到地震反演领域，会发现全变分正则化在清晰勾勒出地层分界面的同时会在地层内部光滑区域产生阶梯效应。学者们研究发现，通过改进全变分模型，可以减弱阶梯效应。在对全变分模型的改进方面，最有名的当属广义全变分(Total Generalized Variation，TGV)模型。广义全变分除了使用待反演参数的一阶偏导信息外，还使用了二阶甚至更高阶信息。Bredies等人(2010)在名为SIAM Journal on Imaging Sciences的期刊，公开了题为“Total GeneralizedVariation”的文献。他们在该文献中，通过实验给出了广义全变分正则化比全变分正则化更能保持图像内部区域光滑的结论。但是，基于广义全变分(TGV)正则化的反问题的求解比较困难。因为广义全变分除了含有一阶项外，还包含二阶的复杂形式，这导致目标函数的形式也比较复杂。这一点阻止了广义全变分正则化约束在地震反演领域的推广使用。

综合上述，现有的地震反演方法存在如下问题：(1)现在的全变分正则化会在地层内部产生阶梯效应；(2)现有的多道同时反演方法计算复杂度和空间复杂度高。

本发明巧妙地利用交替方向乘子法(Alternating Direction Method ofMultipliers，ADMM)解决了广义全变分正则化约束的地震反演的求解问题,同时本发明在频域直接反演，能降低计算复杂度和空间复杂度，提高反演速度。

发明内容

本发明的目的在于：针对现有地震反演技术采用全变分正则化会在地层内部产生阶梯效应的问题，本申请提供了一种基于广义全变分正则化的地震反演方法及系统。

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

本申请提供了一种基于广义全变分正则化的地震反演方法，建立关于待反演参数的对数的广义全变分正则化约束的目标函数，再求解目标函数得到反演参数；

上述方案中，具体地，利用交替方向乘子法求解目标函数得到反演参数；

上述方案中，具体步骤为：

步骤1：构建目标函数J：目标函数J包括子波矩阵W、待反演参数经对数运算后的矩阵L、差分矩阵D、原始地震记录S所组成的数据保真项、初始模型约束项以及待反演参数的对数的广义全变分TGV(L)，差分矩阵D是待反演参数的自然对数到其对应的反射系数之间的转换矩阵且为托普利兹矩阵；

目标函数J为：

其中W为子波矩阵；L为待反演参数的对数；S为观测到的地震记录矩阵；L′为待反演参数初始模型的对数；μ为广义全变分正则化系数；η为初始模型约束正则化系数；v为辅助变量；TGV(L)表示L的广义全变分，D为差分矩阵，α₀、α₁分别为L的一、二阶偏导项的权重系数；

D为差分矩阵，其形式为：

ξ(v)的具体形式为：

ν_x是v在x方向的偏导，ν_y是v在y方向的偏导，是求水平方向的偏导的算子，是求垂直方向的偏导的算子；

步骤2：变换目标函数J：将目标函数J中的广义全变分TGV(L)变换为矩阵形式，并将子波矩阵W、差分矩阵D和待反演参数的自然对数L三者的乘积WDL变换为子波矩阵对应的卷积核w、差分矩阵对应的卷积核d以及待反演参数的自然对数L三者的卷积对变换后的目标函数J再进行傅里叶变换得到目标函数的频域表达式；

步骤2中变换后的目标函数J为：

其中K_h＝[1,-1],K_v＝[1,-1]^T；w为子波对应的列向量，也是子波矩阵W对应的卷积核；d＝[0.5， -0.5]^T是D矩阵对应的卷积核；符号表示二维矩阵卷积运算；

目标函数的频域表达式为：

式中，表示L的傅里叶变换结果，为d的傅里叶变换结果，为w的傅里叶变换结果，为S的傅里叶变换结果，为V_y的傅里叶变换结果，为V_x的傅里叶变换结果，为K_v的傅里叶变换结果，为K_h的傅里叶变换结果，符号“ο”表示点乘算子；

步骤3：求解反演结果：利用交替方向乘子法对步骤2中的目标函数的频域表达式进行求解并迭代更新待反演参数的对数，当待反演参数的自然对数L的L₂范数的相对变化小于预先设定的阈值tol时输出反演结果；

所述步骤3的具体步骤为：

步骤3.1：引入辅助变量并对辅助变量及正则化参数数值进行初始化：

引入辅助变量Z_i(i＝1,2,…,5)，C_i(i＝1,2,…,5)；

并初始化辅助变量及正则化参数数值：

Z_i＝0(i＝1,2,…,5)；C_i＝0(i＝1,2,…,5)；L＝L',μ₀,β₀,μ₁,β₁,η,tol；

步骤3.2：更新L,V_x,V_y；

运用ADMM算法可以得到L，V_x，V_y的求解公式为：

其中，符号“/”表示按元素点除，表示2维傅里叶逆变换，如下所示：

其中，式中各符号的意义为：

其中“*”表示共轭运算；

步骤3.3：更新辅助变量Z_i(i＝1,2,…,5)：

步骤3.4：更新辅助变量C_i(i＝1,2,…,5)：

步骤3.4：输出最终的反演结果：当待反演参数的自然对数L的L₂范数的相对变化小于预先设定的阈值tol时，迭代过程结束，反演得到的波阻抗为：exp(L)，“exp()”表示求指数运算。

本申请还提供了一种基于广义全变分正则化的地震反演系统，包括：

目标函数生成模块：用于生成具有待反演参数的对数的广义全变分正则化约束的目标函数J；

所述目标函数J包括子波矩阵W、待反演参数经对数运算后的矩阵L、差分矩阵D、原始地震记录S所组成的数据保真项、初始模型约束项以及待反演参数的对数的广义全变分TGV(L)，差分矩阵D是待反演参数的自然对数到其对应的反射系数之间的转换矩阵且为托普利兹矩阵；

所述目标函数变换模块：将广义全变分TGV(L)变换为矩阵形式，并将子波矩阵W、差分矩阵D和待反演参数的自然对数L三者的乘积WDL变换为子波矩阵对应的卷积核w、差分矩阵对应的卷积核d以及待反演参数的自然对数L三者的卷积并对目标函数J进行变换，再对变换后的目标函数J行傅里叶变换得到目标函数的频域表达式；

目标函数J：

其中W为子波矩阵；L为待反演参数的对数；S为观测到的地震记录矩阵；L′为待反演参数初始模型的对数；μ为广义全变分正则化系数；η为初始模型约束正则化系数；v为辅助变量；TGV(L)表示L的广义全变分，D为差分矩阵，α₀、α₁分别为L的一、二阶偏导项的权重系数；

D为差分矩阵，其形式为：

ξ(v)的具体形式为：

^ν _x是v在x方向的偏导，ν_y是v在y方向的偏导，是求水平方向的偏导的算子，是求垂直方向的偏导的算子；

目标函数变换模块：对目标函数J进行傅里叶变换得到目标函数的频域表达式；

将目标函数J变换后的目标函数J为：

其中K_h＝[1,-1],K_v＝[1,-1]^T；w为子波对应的列向量，也是子波矩阵W对应的卷积核；d＝[0.5，-0.5]^T是D矩阵对应的卷积核；符号表示二维矩阵卷积运算；

将目标函数J做傅里叶变换后的目标函数的频域表达式为：

式中，表示L的傅里叶变换结果，为d的傅里叶变换结果，为w的傅里叶变换结果，为s的傅里叶变换结果，为V_y的傅里叶变换结果，为V_x的傅里叶变换结果，为K_v的傅里叶变换结果，为K_h的傅里叶变换结果，符号“ο”表示点乘算子；

反演参数求解模块：利用交替方向乘子法对目标函数的频域表达式进行求解并迭代不断更新待反演参数的对数，当待反演参数的自然对数L的L₂范数的相对变化小于预先设定的阈值tol时输出反演结果；

其中反演参数求解模块的具体步骤：

步骤3.1：引入辅助变量并对辅助变量及正则化参数数值进行初始化：

引入辅助变量Z_i(i＝1,2,…,5)，C_i(i＝1,2,…,5)；

并初始化辅助变量及正则化参数数值：

Z_i＝0(i＝1,2,…,5)；C_i＝0(i＝1,2,…,5)；L＝L',μ₀,β₀,μ₁,β₁,η,tol；

步骤3.2：更新L,V_x,V_y；

运用ADMM算法可以得到L，V_x，V_y的求解公式为：

其中，符号“/”表示按元素点除，表示2维傅里叶逆变换，如下所示：

其中，式中各符号的意义为：

其中“*”表示共轭运算；

步骤3.3：更新辅助变量Z_i(i＝1,2,…,5)：

步骤3.4：更新辅助变量C_i(i＝1,2,…,5)：

步骤3.4：输出最终的反演结果：当L的L₂范数的相对变化小于预先设定的阈值tol时，迭代过程结束，反演得到的波阻抗为：exp(L)，“exp()”表示求指数运算。

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

1.本发明中方法利用广义全变分除了使用待反演参数的一阶偏导信息外，还使用了二阶甚至更高阶信息这一特点，因此能够在清晰勾勒地层分界面的同时，减弱地层内部的阶梯效应，同时本发明可以由目标函数直接求得待反演参数值，而无需先反演反射系数，在频域直接反演得到待反演参数，应用广义全变分正则化可以在保持地层分界面的同时，减弱阶梯效应的影响，使得地层内部区域较光滑；

2.本发明巧妙地利用FFT变换，在频域进行反演可以降低算法的空间复杂度和时间复杂度，提高反演效率；

3.本发明中从待反演参数的自然对数L矩阵到其对应的反射系数之间的转换矩阵也就是D矩阵是Toeplitz矩阵，而该类型矩阵与其他矩阵相乘等价于该类型矩阵对应的卷积核与其他矩阵的卷积，从而使得本申请巧妙地利用频域的点乘运算来代替时域的矩阵相乘；通过利用子波矩阵及各个差分矩阵所具有的卷积性质，将时域的卷积运算转化为频域的点乘运算，可以大大提高反演速度；

4.本发明具有普遍适用性，适用于叠后参数反演或叠前弹性参数同步反演，同时适用于二阶甚至更高阶广义全变分正则化的地震反演。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。附图并未刻意按实际尺寸等比例缩放绘制附图，重点在于示出本发明的主旨。

图1本发明的流程图；

图2本发明用于测试的波阻抗真实值图；

图3本发明用于测试的波阻抗初始模型；

图4本发明合成地震数据添加20％的高斯白噪声剖面图；

图5本发明用于测试的子波；

图6本发明的波阻抗反演结果图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例一

本实施例提供了一种基于广义全变分正则化的地震反演方法，其中利用二阶广义全变分为例进行说明，结合图1对本发明方法的流程进行说明：

步骤1：输入和创建反演所需的数据

步骤1.1：输入地震数据、测井资料、层位解释信息；

步骤1.2：根据步骤1.1输入的数据提取地震子波W，并创建待反演参数的初始模型；

步骤2：建立和求解目标函数

目标函数J为：

其中W为子波矩阵；L为待反演参数的对数；S为观测到的地震记录矩阵；L′为待反演参数初始模型的对数；μ为广义全变分正则化系数；η为初始模型约束正则化系数；v为辅助变量；TGV(L)表示L的广义全变分，D为差分矩阵，α₀、α₁分别为L的一、二阶偏导项的权重系数；

D为差分矩阵，其形式为：

ξ(v)的具体形式为：

ν_x是v在x方向的偏导，ν_y是v在y方向的偏导，是求水平方向的偏导的算子，是求垂直方向的偏导的算子；

为了使用傅里叶变换提高算法速度，我们将TGV正则项写成矩阵形式，并将差分矩阵和子波矩阵表示为卷积形式，即：

其中K_h＝[1,-1],K_v＝[1,-1]^T；w为子波对应的列向量，也是子波矩阵W对应的卷积核；d＝[0.5， -0.5]^T是D矩阵对应的卷积核；符号表示二维矩阵卷积运算；

步骤2.2：对变换后的目标函数J进行傅里叶变换得到目标函数的频域表达式：

目标函数的频域表达式为：

式中，表示L的傅里叶变换结果，为d的傅里叶变换结果，为w的傅里叶变换结果，为S的傅里叶变换结果，为V_y的傅里叶变换结果，为V_x的傅里叶变换结果，为K_v的傅里叶变换结果，为K_h的傅里叶变换结果，符号“ο”表示点乘算子；

步骤2.3：引入辅助变量并对辅助变量及正则化参数数值进行初始化：

引入辅助变量Z_i(i＝1,2,…,5)，C_i(i＝1,2,…,5)；

并初始化辅助变量及正则化参数数值：

Z_i＝0(i＝1,2,…,5)；C_i＝0(i＝1,2,…,5)；L＝L',μ₀,β₀,μ₁,β₁,η,tol；

步骤2.4：更新L,V_x,V_y；

运用ADMM算法可以得到L，V_x，V_y的求解公式为：

其中，符号“/”表示按元素点除，表示2维傅里叶逆变换，如下所示：

其中，式中各符号的意义为：

其中“*”表示共轭运算；

步骤2.4.1：更新辅助变量Z_i(i＝1,2,…,5)：

步骤2.4.2：更新辅助变量C_i(i＝1,2,…,5)：

步骤3：输出最终的反演结果：当待反演参数的自然对数L的L₂范数的相对变化小于预先设定的阈值tol时，迭代过程结束，反演得到的波阻抗为：exp(L)，“exp()”表示求指数运算。

实施例中的波阻抗真实模型(如图2)取自Marmousi2模型的一部分：深度方向共有393个采样点，采样间隔为2.5m，距离方向共有453道，相邻道之间的距离是3.75m；

波阻抗初始模型(如图3)是对真实模型的高斯低通滤波；用于反演的地震数据(如图4)是用波阻抗真实模型求得反射系数，然后用反射系数与主频为40Hz的雷克子波褶积，再添加20％的高斯白噪声而得到。实施例中的子波(如图5)为主频为40Hz的雷克子波；波阻抗反演结果如图6所示；

参考图2与图6可以看出，采用本发明方法得到的波阻抗反演结果(如图6)与真实的波阻抗真实模型(图2)十分相近，且从图6中可以看到地层的分界面刻画得十分清晰。

实施例二

本实施例提供了一种基于上述实施例的基础的一种基于广义全变分正则化的地震反演系统，主要包括：

目标函数生成模块：用于生成具有待反演参数的对数的广义全变分正则化约束的目标函数J；

目标函数变换模块：对目标函数J进行傅里叶变换得到目标函数的频域表达式；

反演参数求解模块：利用交替方向乘子法对目标函数的频域表达式进行求解并迭代不断更新参数直到待反演参数的自然对数L的L₂范数的相对变化小于预先设定的阈值tol时输出反演结果。

实施例三

在实施例二的基础上，进一步地详细说明本申请的系统：

一、目标函数生成模块的工作流程如下：

需要输入地震数据、测井资料、层位解释信息，同时还要根据数据提取地震子波W，并创建待反演参数的初始模型，再根据上述参数创建目标函数J：

其中W为子波矩阵；L为待反演参数的对数；S为观测到的地震记录矩阵；L′为待反演参数初始模型的对数；μ为广义全变分正则化系数；η为初始模型约束正则化系数；v为辅助变量；TGV(L)表示L的广义全变分，D为差分矩阵，α₀、α₁分别为L的一、二阶偏导项的权重系数；

D为差分矩阵，其形式为：

ξ(v)的具体形式为：

ν_x是v在x方向的偏导，ν_y是v在y方向的偏导，是求水平方向的偏导的算子，是求垂直方向的偏导的算子；

二、目标函数变换模块的工作流程如下：

1：将TGV正则项写成矩阵形式，并将差分矩阵和子波矩阵表示为卷积形式，也即：

其中K_h＝[1,-1],K_v＝[1,-1]^T；w为子波对应的列向量，也是子波矩阵W对应的卷积核；d＝[0.5， -0.5]^T是D矩阵对应的卷积核；符号表示二维矩阵卷积运算；

2：对变换后的目标函数J进行傅里叶变换得到目标函数的频域表达式：

式中，表示L的傅里叶变换结果，为d的傅里叶变换结果，为w的傅里叶变换结果，为S的傅里叶变换结果，为V_y的傅里叶变换结果，为V_x的傅里叶变换结果，为K_v的傅里叶变换结果，为K_h的傅里叶变换结果，符号“ο”表示点乘算子；

3：引入辅助变量Z_i(i＝1,2,…,5)，C_i(i＝1,2,…,5)，并对辅助变量、正则化参数等数值进行初始化：

Z_i＝0(i＝1,2,…,5)；C_i＝0(i＝1,2,…,5)；L＝L',μ₀,β₀,μ₁,β₁,η,tol；

三、反演参数求解模块的工作流程如下：

1：求解目标函数

更新L,V_x,V_y；

运用ADMM算法可以得到L，V_x，V_y的求解公式为：

其中，符号“/”表示按元素点除，表示2维傅里叶逆变换，如下所示：

其中，式中各符号的意义为：

其中“*”表示共轭运算；

1.2：更新辅助变量Z_i(i＝1,2,…,5)：

1.3：更新辅助变量C_i(i＝1,2,…,5)：

步骤1.4：输出最终的反演结果，当待反演参数的自然对数L的L₂范数的相对变化小于预先设定的阈值tol时，迭代过程结束，反演得到的参数结果为：exp(L)，“exp()”表示求指数运算。

本系统中方法利用广义全变分除了使用待反演参数的一阶偏导信息外，还使用了二阶甚至更高阶信息这一特点，因此能够在清晰勾勒地层分界面的同时，减弱地层内部的阶梯效应，同时本发明可以由目标函数直接求得待反演参数值，而无需先反演反射系数，在频域直接反演得到待反演参数，用广义全变分正则化可以在保持地层分界面的同时，减弱阶梯效应的影响，使得地层内部区域较光滑。

其中应当说明的是：

初始模型约束项在文中具体指的是：

Toeplitz矩阵是指文中提到的托普利兹矩阵；

FFT指的是傅里叶变换；

其中本文中“辅助变量、正则化参数数值进行初始化”中的“辅助变量”有两个，分别为：Z_i(i＝1,2,…,5)，C_i(i＝1,2,…,5)；正则化参数有：μ₀,β₀,μ₁,β₁,η；

其中μ₀,β₀,μ₁,β₁,η视具体的地震数据，经过与井数据的测试对比之后设定；

其中上述中的tol为预先设定的迭代阈值。

S为观测到的地震记录，根据本领域人员公知常识可以得出S的内容，同时本申请的实施例内容的S采用的是二维矩阵，但是S并不限于二维；

本实施例中的广义变分正则化采用的是二阶，但是本方法并不限于二阶；

当然系统在具体设计与执行的时候并不一定按照上述内容进行划分执行，但是在不脱离本系统的思想的技术方案均应在保护范围内。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何属于本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

Claims (8)

1.一种基于广义全变分正则化的地震反演方法，其特征在于，建立关于待反演参数的对数的广义全变分正则化约束的目标函数，再求解目标函数得到反演参数；

利用交替方向乘子法求解目标函数得到反演参数；

具体步骤为：

步骤1：构建目标函数J：目标函数J包括子波矩阵W、待反演参数经对数运算后的矩阵L、差分矩阵D、原始地震记录S所组成的数据保真项、初始模型约束项以及待反演参数的对数的广义全变分TGV(L)，差分矩阵D是待反演参数的自然对数到其对应的反射系数之间的转换矩阵且为托普利兹矩阵；

步骤2：变换目标函数J：将目标函数J中的广义全变分TGV(L)变换为矩阵形式，并将子波矩阵W、差分矩阵D和待反演参数的自然对数L三者的乘积WDL变换为子波矩阵对应的卷积核w、差分矩阵对应的卷积核d以及待反演参数的自然对数L三者的卷积对变换后的目标函数J再进行傅里叶变换得到目标函数的频域表达式；

步骤3：求解反演结果：利用交替方向乘子法对步骤2中的目标函数的频域表达式进行求解并迭代更新待反演参数的对数，当待反演参数的自然对数L的L₂范数的相对变化小于预先设定的阈值tol时输出反演结果。

2.如权利要求1所述的一种基于广义全变分正则化的地震反演方法，其特征在于，所述步骤1中的目标函数J为：

其中W为子波矩阵；L为待反演参数的对数；S为观测到的地震记录矩阵；L'为待反演参数初始模型的对数；μ为广义全变分正则化系数；η为初始模型约束正则化系数；v为辅助变量；TGV(L)表示L的广义全变分，D为差分矩阵，α₀、α₁分别为L的一、二阶偏导项的权重系数；

D为差分矩阵，其形式为：

ξ(v)的具体形式为：

ν_x是v在x方向的偏导，ν_y是v在y方向的偏导，是求水平方向的偏导的算子，是求垂直方向的偏导的算子。

3.如权利要求1所述的一种基于广义全变分正则化的地震反演方法，其特征在于，所述步骤2中变换后的目标函数J为：

其中K_h＝[1,-1],K_v＝[1,-1]^T；w为子波对应的列向量，也是子波矩阵W对应的卷积核；d＝[0.5，-0.5]^T是D矩阵对应的卷积核；符号表示二维矩阵卷积运算。

4.如权利要求1所述的一种基于广义全变分正则化的地震反演方法，其特征在于，目标函数的频域表达式为：

式中，表示L的傅里叶变换结果，为d的傅里叶变换结果，为w的傅里叶变换结果，为S的傅里叶变换结果，为V_y的傅里叶变换结果，为V_x的傅里叶变换结果，为K_v的傅里叶变换结果，为K_h的傅里叶变换结果，符号“o”表示点乘算子。

5.如权利要求4所述的一种基于广义全变分正则化的地震反演方法，其特征在于，所述步骤3的具体步骤为：

步骤3.1：引入辅助变量并对辅助变量及正则化参数数值进行初始化：

引入辅助变量Z_i(i＝1,2,…,5)，C_i(i＝1,2,…,5)；

并初始化辅助变量及正则化参数数值：

Z_i＝0(i＝1,2,…,5)；C_i＝0(i＝1,2,…,5)；L＝L',μ₀,β₀,μ₁,β₁,η,tol；

步骤3.2：更新L,V_x,V_y；

运用ADMM算法可以得到L，V_x，V_y的求解公式为：

其中，符号“/”表示按元素点除，表示2维傅里叶逆变换，如下所示：

其中，式中各符号的意义为：

其中“*”表示共轭运算；

步骤3.3：更新辅助变量Z_i(i＝1,2,…,5)：

步骤3.4：更新辅助变量C_i(i＝1,2,…,5)：

步骤3.4：输出最终的反演结果：

迭代过程结束，反演得到的波阻抗为：exp(L)，“exp()”表示求指数运算。

6.一种基于广义全变分正则化的地震反演系统，其特征在于，包括：

目标函数生成模块：用于生成具有待反演参数的对数的广义全变分正则化约束的目标函数J；

目标函数变换模块：对目标函数J进行傅里叶变换得到目标函数的频域表达式；

反演参数求解模块：利用交替方向乘子法对目标函数的频域表达式进行求解并迭代不断更新待反演参数的对数，当待反演参数的自然对数L的L₂范数的相对变化小于预先设定的阈值tol时输出反演结果。

7.如权利要求6所述的一种基于广义全变分正则化的地震反演系统，其特征在于，所述目标函数J包括子波矩阵W、待反演参数经对数运算后的矩阵L、差分矩阵D、原始地震记录S所组成的数据保真项、初始模型约束项以及待反演参数的对数的广义全变分TGV(L)，差分矩阵D是待反演参数的自然对数到其对应的反射系数之间的转换矩阵且为托普利兹矩阵；

所述目标函数变换模块：将广义全变分TGV(L)变换为矩阵形式，并将子波矩阵W、差分矩阵D和待反演参数的自然对数L三者的乘积WDL变换为子波矩阵对应的卷积核w、差分矩阵对应的卷积核d以及待反演参数的自然对数L三者的卷积并对目标函数J进行变换，再对变换后的目标函数J行傅里叶变换得到目标函数的频域表达式。

8.如权利要求7所述的一种基于广义全变分正则化的地震反演系统，其特征在于，所述目标函数J：

其中W为子波矩阵；L为待反演参数的对数；S为观测到的地震记录矩阵；L'为待反演参数初始模型的对数；μ为广义全变分正则化系数；η为初始模型约束正则化系数；v为辅助变量；TGV(L)表示L的广义全变分；D为差分矩阵；α₀、α₁分别为L的一、二阶偏导项的权重系数；

D为差分矩阵，其形式为：

ξ(v)的具体形式为：

ν_x是v在x方向的偏导，ν_y是v在y方向的偏导，是求水平方向的偏导的算子，是求垂直方向的偏导的算子；

将目标函数J变换后的目标函数J为：

其中K_h＝[1,-1],K_v＝[1,-1]^T；w为子波对应的列向量，也是子波矩阵W对应的卷积核；d＝[0.5，-0.5]^T是D矩阵对应的卷积核；符号表示二维矩阵卷积运算；将目标函数J做傅里叶变换后的目标函数的频域表达式为：

式中，表示L的傅里叶变换结果，为d的傅里叶变换结果，为w的傅里叶变换结果，为s的傅里叶变换结果，为V_y的傅里叶变换结果，为V_x的傅里叶变换结果，为K_v的傅里叶变换结果，为K_h的傅里叶变换结果，符号“o”表示点乘算子。