CN105549076B - 一种基于交替方向法和全变分理论的地震数据处理方法 - Google Patents

Abstract

本发明提出一种基于交替方向法和全变分理论的地震数据处理方法，该方法建立基于全变分的地震数据模型，利用经典的二次罚函数方法将此模型转换为满足能够迭代与交替最小化的求解问题，通过交替方向法来迭代获得恢复后的地震数据。该技术可以同时对地震数据进行去噪和去模糊(反褶积)处理，能够有效地提高地震数据的信噪比和分辨率，即便是在信号被严重噪声污染的情况下，也能够得到较好的处理效果。

Description

一种基于交替方向法和全变分理论的地震数据处理方法

技术领域

本发明属于地球物理勘探领域，涉及一种地震数据处理的方法，具体涉及一种基于交替方向法和全变分理论的地震数据处理方法，其适用于对地球物理勘探数据提高信噪比和分辨率的处理。

背景技术

随着精细的油气地震勘探技术对地震数据信噪比和分辨率要求越来越高，在处理地震数据的同时保持高信噪比和高分辨率已成为本行业的一项重大课题，它直接影响到能源勘探的技术进步。目前，传统提高地震数据信噪比和分辨率技术已经难以适应当前高精度勘探的需求，并且二者一直相互制约。因此，寻找并开发出适应于当前复杂地质条件下的高信噪比和高分辨率处理技术，并使之有助于提高油气勘探的精度具有更现实的意义。然而，由于地震数据处理方面长期没有取得明显的进展，寻找并开发出适应于当前复杂地质条件下的高信噪比和高分辨率处理技术还仍然是本领域一项等待填补空白的前沿技术。此类方法的研究与技术开发不但是本学科重大课题，在应用领域也具有影响能源勘探技术进步的重大意义。

目前，国内外研究学者和机构提出了多种不同的去噪和反褶积方法，以提高地震资料的信噪比和分辨率，如维纳(Wiener)滤波，中值滤波，f-k滤波、f-x反褶积、K-L变换、奇异值分解(Singular Value Decomposition,SVD)、基于现代信号分析的小波(Wavelet)变换、曲波(Curvelet)变换、非二次幂Curvelet变换、离散余弦变换(Discrete CosineTransform, DCT)、基于独立成分分析的信号分离、基于非局部均值的去噪算法以及常规的L¹范数稀疏脉冲反褶积、参数化稀疏脉冲反褶积方法等[①Treitel S.The complexWiener filter.Geophysics,1974,39(2):169-173.②李月,马海涛,林红波等.基于核函数主成分的维纳滤波方法研究.地球物理学报,2010,53(5):1226-1233.③刘洋,王典,刘财等.局部相关加权中值滤波技术及其在叠后随机噪声衰减中的应用.地球物理学报,2011,54(2):358-367.④Stewart R R,Schieck D G.3-D f-k filtering.Journal of SeismicExploration,1993,2:41-54.⑤Canales L L.Random noise reduction.SEG TechnicalProgram Expanded Abstracts,1984,3.⑥Shen H Y,Li Q C.Seismic wave fieldssparation and noise attenuation in frequeccy domain var singular valuedecomposition.Near-Surface Geophysics and Human Activity,2008:178-181.⑦张恒磊,张云翠,宋双等.基于Curvelet域的叠前地震资料去噪方法.石油地球物理勘探,2008,43(5):508-513.⑧Bonar D,Sacchi M.Denoising seismic data using the nonlocalmeans algorithm.Geophysics,2012,77(1):A5-A8.⑨Danilo R V.Parametric sparse-spike deconvolution and the recovery of the acoustic impedance.The 76thAnnual International Meeting of Society of Exploration Geophysicists,NewOrleans,Tulsa:Society of Exploration Geophysicists,2006:2141-2144.]。

在提高信噪比方面：1)基于傅立叶变换的去噪技术：一些学者研究了关于f-x预测技术，能够有效地压制随机噪声，衰减面波和线性干扰，保护和加强有效波高频成分，为后续高分辨率处理打下良好基础。一些学者提出了三维f-k滤波技术，对分离出面波和线性干扰有较好的效果，但该类技术存在丢失有效信号的缺陷。2)基于小波分解与重建的去噪技术：一些学者利用小波变换具有分频和局部分析能力的特点，对地震资 料中的噪声进行压制，很好地提高了地震资料的信噪比。一些学者利用二维小波变换易于分离噪声和有效信号的特点，对地震资料中面波进行了消除，取得了满意的效果。3)基于曲波Curvelet变换的去噪技术：一些学者利用Curvelet变换具有良好的方向特性，将地震数据变换到曲波域，并运用窗口萎缩算法，较好地压制了随机噪声，提高了地震资料的质量。一些学者提出了非二次幂Curvelet变换技术，对地震资料中面波进行了处理测试，获得了较好的测试效果。4)基于奇异值分解SVD的去噪技术：一些学者研究了奇异值分解SVD算法在地震资料处理中的应用。针对感兴趣的地震信号成分进行SVD波场分离与去噪，有效地提高了地震资料的信噪比。一些学者研究了局部频率域SVD噪声压制方法，能够较好地对倾斜或弯曲同向轴进行处理，获得了满意的去噪效果。

在提高分辨率方面：1)一些学者研究了关于L1和L2范数的稀疏脉冲反褶积技术，使用该技术可以使分辨率得到较好地改善。然而，该类方法在提高分辨率的同时往往会降低信噪比。2)基于曲波Curvelet变换的稀疏反褶积技术：一些学者提出了基于曲波Curvelet变换的稀疏反褶积技术。曲波变换对多维信号的稀疏表示，能获得最优的非线性逼近。在提高地震分辨率的同时，保持了同相轴的连续性。3)基于维纳滤波理论的反褶积技术：该技术是是各商业处理系统中经常使用的处理模块，其根据期望输出的不同，会得到不同类型的反褶积，如最小平方反褶积、预测反褶积、脉冲反褶积、子波整形反褶积。该类反褶积的应用需要以下几点假设：①地震子波一致；②反射系数呈高斯白噪；③地震子波为最小相位；④无噪声；在满足以上假设的条件下，该类反褶积均会得到满意的效果。然而，由于地下介质中传播的地震子波通常是混合相位，地层反射系数也不呈高斯白噪分布，因此，会限制该类技术的应用。此外，该类技术往往在提高分辨率的同时，会降低信噪比，即信噪比和分 辨率存在相互制约的影响。4)地震盲反褶积技术：近些年，为了解决传统反褶积技术其假设条件与实际情况不符的问题，一些学者研究了地震盲反褶积问题。该类技术不需要对地震子波和反射系数进行先验假设，仅从观测数据中分离子波，提取反射系数。一些学者研究了基于互信息率的盲反褶积方法，适应于非最小相位子波，能够有效拓展地震记录的频带，提高地震资料的分辨率。一些学者提出了基于粒子群算法的盲反褶积，适应于最小相位子波及混合相位子波的反褶积，能够更好地从地震数据中估计反射系数，有效拓宽地震资料的频谱，得到高分辨率的地震资料。因此，研究地震盲反褶积技术将具有重要的现实意义和发展前景。

通过对以上现有技术研究与分析，发现目前该领域主要存在一个共性或关键技术瓶颈问题：常规提高地震资料信噪比和分辨率技术存在二者相互制约的问题，即在提高分辨率的同时，往往会降低信噪比。两者相互制约影响的弊端已成为地震资料处理领域长期困扰的技术难题。

因此，为了更好地提高地震数据的品质质量，急需一种能够解决上述问题的地震数据处理技术。

发明内容

为了解决上述问题，我们提出一种基于交替方向法和全变分理论的地震数据处理方法，其在提高地震数据的信噪比的同时，可以保持高的数据分辨率。

依据本发明的技术方案，提供一种基于交替方向法和全变分理论的地震数据处理方法，其基于建立基于全变分的地震数据模型，利用经典的二次罚函数方法将此模型转换为满足能够迭代与交替最小化的求解问题，通过交替方向法来迭代获得恢复后的地震数据。

其主要包括以下步骤：首先在野外勘探目标区中在地表以人工方法激发地震波，利用地震采集设备，获得地震数据；将野外采集到的地震数据进行格式转换，满足数据处理的格式需求；

然后采用下述步骤进行数据处理：

步骤1：建立基于全变分的地震数据模型；将全变分理论融入到带有噪声的原始地震反褶积模型中，建立成基于全变分的地震数据模型，能免于受噪声的影响，其目标是：在给定合适的褶积算子情况下，可从含有噪声的观测地震数据中恢复原始地震数据；

步骤2：对建立的地震数据恢复模型(目标函数)进行重写和转化，即利用二次罚函数方法将“有约束优化”转化为“无约束优化”问题，使其能够迭代与交替方向最小化求解；

步骤3：利用交替方向法(ADM)对目标函数(8)进行迭代与交替方向最小化求解；即在步骤2的基础上，利用交替方向法(ADM)对转化后的模型进行迭代与交替方向最小化求解，这种求解方法可以得到每一步迭代过程中的近似最优解；

步骤4：利用第三步获得数据，通过交替方向法来迭代获得恢复后的地震数据，能够在合理的时间内得到全局最优解。

更近一步地，所述建立基于全变分的地震数据模型具体为：设有n个地震道记录的地震剖面，每个地震道记录有m个采样点，则用矩阵r来表示此n·m个数据，其元素为x_ij(i＝1,2,…,n；j＝1,2,…,m)，则地震数据r被表示为：

令x^*＝vec(r)为矩阵r的拉伸向量，即x^*∈R^nm；K∈R^nm×nm为褶积算 子，w∈R^nm为添加的随机噪声，f∈R^nm为观测的地震数据。满足如下关系：

f＝Kx^*+w. (2)

对于给定的褶积算子K，其目标是从含有噪声的观测地震数据f中恢复原始地震数据x^*，即视为去噪和去模糊(反褶积)处理。若K为单位算子，则此操作则视为去噪处理。由于||Kx-f||₂≈||w||₂，即恢复原始地震数据x^*为最小化公式(3)：

其中，||·||₂表示L₂范数，Φ_fid(x,f)称为随机噪声的数据保真项。众所周知，对于公式(2)，直接从f中恢复x^*为不适定问题，因为它极其敏感于噪声w。因此，为了使恢复x^*具有稳定性，引入TV正则项，即全变分(TV)的离散形式：

其中，D_i为局部有限差分算子，D_ix∈R²为x在第i个采样点处包含水平和垂直两个方向的一阶有限差分。TV正则化也有利于保留数据锋利的边缘和边界。因此，地震数据的恢复为解决如下最小化问题：

即地震数据的重建模型为：

其中，μ>0为平衡二者的参数。注意：x^*为原始正演模型或实际地震资料理论上的原始数据；x为被恢复的实际地震数据。

本发明所提出的基于交替方向法和全变分理论的地震数据处理方法，通过建立基于全变分的地震数据模型，利用经典的二次罚函数方法将此模型转换为满足能够迭代与交替最小化的求解问题，通过交替方向法来迭代获得恢复后的地震数据。本发明方法可以同时对地震数据进行去噪和去模糊(反褶积)处理，能够有效地提高地震数据的信噪比和分辨率，即便是在信号被严重噪声污染的情况下，也能够得到较好的处理效果。

附图描述

图1为人工合成的原始地震正演模型。

图2为利用Gaussian核函数与原始正演模型(图1)进行褶积后的模型，Gaussian核函数的标准差σ＝2。

图3(a)为本发明方法对图2模型的处理结果(μ＝5)；图3(b)为本发明方法对图2模型进行处理，其迭代次数与信噪比SNR及目标函数值的关系曲线。

图4为Wiener滤波对图2模型的处理结果。

图5为Laplacian滤波对图2模型的处理结果。

图6是利用Gaussian核函数与原始模型(图1)进行褶积后，再由MATLAB中“randn”函数添加随机噪声后的模型。

图7(a)为本发明方法对图6模型的处理结果(μ＝4)；

图7(b)为本发明方法对图6模型进行处理，其迭代次数与信噪比SNR及目标函数值的关系曲线。

图8为Wiener滤波对图6模型的处理结果。

图9为Laplacian滤波对图6模型的处理结果。

图10(a)、图10(b)、图10(c)和图10(d)分别是图6、图7(a)、图8和 图9模型的第130道频谱图。

图11是在原始正演模型(图1)上加入了很强的随机噪声得到的模型。

图12(a)为本发明方法对图11模型的处理结果(μ＝4)；图12(b)为本发明方法对图11模型进行处理，其迭代次数与信噪比SNR及目标函数值的关系曲线。

图13(a)为一实际地震资料初始叠加剖面；图13(b)为本发明方法对图13(a)处理的结果；图13(c)和图13(d)分别为应用某商业软件的径向预测滤波(RadiPredicFilt)和随机噪声衰减(RNA)两个模块对图13(a)处理的结果。

图14(a)是一实际地震资料CDP道集；图14(b)、图14(c)、图14(d)分别为本发明方法、径向预测滤波(RadiPredicFilt)和随机噪声衰减(RNA)模块对图14(a)处理的结果。

图15(a)是一实际地震资料CMP道集；图15(b)是本发明方法对图15(a)处理的结果；

图16(a)是一实际地震资料叠前时间偏移剖面；图16(b)是本发明方法对图16(a)处理的结果。

具体实施方式

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

地震资料数字处理，是指用数字计算机对野外采集的地震信息进行处理和分析，以获得高质量的、可靠的地震剖面，为下一步资料解释提供直观的、可靠的依据和有关的地质信息。本发明是提出一种基于交替方向法和全变分理论的地震数据处理方法。建立基于全变分的地震数据模型，利用经典的二次罚函数方法将此模型转换为满足能够迭代与交替最小化的求解问题，通过交替方向法来迭代获得恢复后的地震数据。

其中，在下面所使用的术语“全变分”为本领域的通用词语，指代在地震勘探中的原始反褶积求解模型中融入数学领域的全变分理论，使其免于受噪声的影响，其优势是它可以从含有噪声的观测地震数据中有效地恢复原始地震数据。术语“交替方向法”指在地震资料处理过程中，求解基于全变分理论的地震数据恢复模型的一种方法，即首先固定模型中三个变量的其中两个变量值，求解另一个变量值；其次，利用刚刚求解得到的变量值和其中一个变量值，反过来求解另一个变量值；最后，利用刚刚求解得到的两个变量值，反过来求解第三个变量值；以此类推，不断迭代得到最优的地震数据恢复结果。

以下各个数学符号的物理意义如下：

m：地震道记录的采样点个数；

n：地震道记录个数；

r：各地震道与每个地震道采样点组成的n×m矩阵；

x_ij：r中的元素，指第i个地震道、第j个采样点值，其中，i＝1,2,…,n；j＝1,2,…,m；

vec(r)：矩阵r的拉伸向量；

x：被恢复的实际地震数据；

K：褶积算子；

f：观测的地震数据；

R_nm：n·m维度空间；

w：添加的随机噪声；

||·||₂：指L₂范数；

Φ_fid(x,f)：随机噪声的数据保真项；

D：一阶全局有限差分算子；

D¹,D²：分别为水平和垂直方向的一阶全局有限差分算子；

D_i：局部有限差分算子，以D¹的第i行为第一行、以D²的第i行为第二行的两行矩阵。

TV(x)：全变分的离散形式；

μ：平衡数据保真项Φ_fid(x,f)和全变分形式TV(x)的参数；

y_i＝D_ix＝((y¹)_i；(y²)_i)：x在第i个采样点处包含水平和垂直两个方向的一阶有限差分，是一2维向量；

y＝(y¹；y²)：2nm维向量；

y¹,y²：长度为n·m维的向量，分别指水平和垂直两个方向的一阶有限差分；

β：二次罚函数方法中的罚参数；

λ：Lagrangian函数参数；

λ_i：R²空间中的向量；

x^k,y^k,λ^k：分别是x,y,λ的第k次迭代更新值；

L_A(x,y,λ)：增广Lagrangian函数；

Ν(·)：表示矩阵的零空间；

F(x^k+1)：指x^k+1的二维离散傅里叶变换；

γ：更新λ过程中被赋予的步长；

ε：算法中用于控制结束条件的公差。

本发明的基于交替方向法和全变分理论的地震数据处理方法，其主题旨在建立基于全变分的地震数据模型，优势是使其免受噪声的影响，更好地实现去噪和反褶积处理，以同时提高地震数据的信噪比和分辨率，模型如下：

其中，n为地震剖面地震道数，m为每个地震道记录采样点数；x为被恢复的实际地震数据；K为褶积算子；f为观测的地震数据；D_i为局部有限差分算子，D_ix∈R²为x在第i个采样点处包含水平和垂直两个方向的一阶有限差分；μ>0为平衡二者的参数。

为了解决基于全变分的地震数据模型，考虑利用交替方向法(AlternatingDirection Method,ADM)进行求解。首先，给出交替方向法(ADM)的基本思想[Gabay D,Mercier B.A dual algorithm for the solution of nonlinear variational problemsvia finite-element approximations.Computing Mathematics Application,1976,2:17-40.]，即令θ₁(·)和θ₂(·)为凸函数，A为连续的线性算子，则最小化能量函数为：

通过引入辅助变量v，上式可等价转换为：

因此，交替方向最小化迭代策略如下：

其中，Θ(u,v,λ)为式(2)的增广Lagrangian函数：

式(3)的交替最小化迭代思想：首先由初始化(v^k,λ^k)获取u^k+1，再由求得的u^k+1和λ^k来获取v^k+1，最后由求得的(u^k+1,v^k+1)来获取λ^k+1。最终通过多次迭代交替方向最小化策略来获取恢复后的地震数据。

通过交替方向法对目标函数进行迭代求解。将基于全变分的地震数据模型重写为式(2)形式，并利用二次罚函数方法将其转化为无约束目标函数：

其中，β>>0为罚参数。利用下式的交替方向法(ADM)对上式目标函数进行求解(初始x＝x^k，λ＝λ^k)：

更具体地，依据本发明的技术方案，一种基于交替方向法和全变分理论的地震数据处理方法，其主要包括以下步骤：首先在野外勘探目标区中在地表以人工方法激发地震波，利用地震采集设备，获得地震数据； 将野外采集到的地震数据进行格式转换，满足数据处理的格式需求；

然后采用下述步骤进行数据处理：

步骤1、将全变分理论融入到带有噪声的原始地震反褶积模型中，建立成基于全变分的地震数据模型，能免于受噪声的影响。其目标是：在给定合适的褶积算子情况下，可从含有噪声的观测地震数据中恢复原始地震数据。

步骤2、对步骤1所建立的地震数据恢复模型(目标函数)进行重写和转化，即利用二次罚函数方法将“有约束优化”转化为“无约束优化”问题，使其能够迭代与交替方向最小化求解。

步骤3、在步骤2的基础上，利用交替方向法(ADM)对转化后的模型进行迭代与交替方向最小化求解，这种求解方法可以得到每一步迭代过程中的近似最优解；

步骤4、设计一种有效的实现地震数据恢复的算法，能够在合理的时间内得到全局最优解。

结合发明内容以及上述步骤概述，给出详细说明和解释，如下：

第一步：建立基于全变分的地震数据模型：设有n个地震道记录的地震剖面，每个地震道记录有m个采样点，则用矩阵r来表示此n·m个数据，其元素为x_ij(i＝1,2,…,n；j＝1,2,…,m)，则地震数据r被表示为：

令x＝vec(r)为矩阵r的拉伸向量，即x∈R^nm；K∈R^nm×nm为褶积算子，w∈R^nm为添加的随机噪声，f∈R^nm为观测的地震数据。满足如下关系：

f＝Kx+w. (6)

对于给定的褶积算子K，其目标是从含有噪声的观测地震数据f中恢复原始地震数据x，即视为去噪和去模糊(反褶积)处理。若K为单位算子，则此操作则视为去噪处理。由于||Kx-f||₂≈||w||₂[Rudin L I,Osher S,Fatemi E.Nonlinear variation based imageremoval algorithms.Physica D,1992,60:259-268.]，即恢复原始地震数据x^*为最小化公式(7)：

其中，||·||₂表示L₂范数，Φ_fid(x,f)称为随机噪声的数据保真项。众所周知，对于公式(6)，直接从f中恢复x为不适定问题，因为它极其敏感于噪声w。因此，为了使恢复x具有稳定性，引入TV正则项[Rudin L I,Osher S,Fatemi E.Nonlinear variation basedimage removal algorithms.Physica D,1992,60:259-268.]，即全变分(TV)的离散形式：

其中，D_i为局部有限差分算子，D_ix∈R²为x在第i个采样点处包含水平和垂直两个方向的一阶有限差分。TV正则化也有利于保留数据锋利的边缘和边界。因此，地震数据的恢复为解决如下最小化问题：

即地震数据的重建模型为：

其中，μ>0为平衡二者的参数。

第二步：对建立的地震数据恢复模型(目标函数)进行重写和转化。具体地，将目标函数(10)重写为式(11)形式，则有：

其中，y＝(y¹；y²)∈R^2nm，y¹和y²是长度为n·m维向量，满足((y¹)_i；(y²)_i)＝y_i∈R²。利用二次罚函数方法[Vogel C R,Oman M E.Fast,robust total variation basedreconstruction of noisy,blurred images.IEEE Transactions on Image Processing,1998,7:813-824.]将(11)转化为无约束目标函数：

其中，β>>0为罚参数。

第三步：利用交替方向法(ADM)对目标函数(12)进行迭代与交替方向最小化求解(初始x＝x^k，λ＝λ^k)：

其中，L_A(x,y,λ)为式(12)的增广Lagrangian函数：

λ_i为R²空间中的向量。式(13)中arg min_yL_A(x^k,y,λ^k)等价于如下形式：

式(15)通过二维收缩来解决：

在固定y^k+1和λ^k基础上，求解式(13)中arg min_xL_A(x,y^k+1,λ^k)为一最小二乘问题，其对应的标准方程为：

其中，D＝(D¹ D²)∈R^2nm×nm

ERROR:rangecheck

OFFENDING COMMAND:string

STACK:

66038

33018

32512

33019。

Claims (1)

1.一种基于交替方向法和全变分理论的地震数据处理方法，其基于建立全变分的地震数据模型，利用经典的二次罚函数方法将此模型转换为满足能够迭代与交替方向最小化的求解问题，通过交替方向法来迭代获得恢复后的地震数据；

所述基于交替方向法和全变分理论的地震数据处理方法主要包括以下步骤：

首先在野外勘探目标区中在地表以人工方法激发地震波，利用地震采集设备，获得地震数据；将野外采集到的地震数据进行格式转换，满足数据处理的格式需求；

然后采用下述步骤进行数据处理：

步骤1：建立基于全变分的地震数据模型；将全变分理论融入到带有噪声的原始地震反褶积模型中，建立成基于全变分的地震数据模型，能免于受噪声的影响，其目标是：在给定合适的褶积算子情况下，可从含有噪声的观测地震数据中恢复原始地震数据；

步骤2：对建立的地震数据恢复目标函数进行重写和转化，即利用二次罚函数方法将“有约束优化”转化为“无约束优化”问题，使其能够迭代与交替方向最小化求解；

步骤3：利用交替方向法(ADM)对目标函数公式(8)进行迭代与交替方向最小化求解；即在步骤2的基础上，利用交替方向法(ADM)对转化后的模型进行迭代与交替方向最小化求解，这种求解方法可以得到每一步迭代过程中的近似最优解；

步骤4：利用第三步获得数据，通过交替方向法来迭代获得恢复后的地震数据，能够在合理的时间内得到全局最优解；

其中，所述建立基于全变分的地震数据模型具体为：设有n个地震道记录的地震剖面，每个地震道记录有m个采样点，则用矩阵r来表示此n·m个数据，其元素为x_ij，i＝1,2,…,n；j＝1,2,…,m，则地震数据r被表示为：

令x^*＝vec(r)为矩阵r的拉伸向量，即x^*∈R^nm；K∈R^nm×nm为褶积算子，w∈R^nm为添加的随机噪声，f∈R^nm为观测的地震数据，满足如下关系：

f＝Kx^*+w (2)

对于给定的褶积算子K，其目标是从含有噪声的观测地震数据f中恢复原始地震数据x^*，即视为去噪和反褶积处理；由于||Kx-f||₂≈||w||₂，恢复原始地震数据x^*为最小化公式(3)：

${\mathrm{\Φ}}_{fid}(x,f)=\left|\right|Kx-f|{|}_2^2---\left(3\right)$

其中，||·||₂表示L₂范数，Φ_fid(x,f)称为随机噪声的数据保真项；对于公式(2)，直接从f中恢复x^*为不适定问题；为了使恢复x^*具有稳定性，引入TV正则项，即全变分(TV)的离散形式：

$TV\left(x\right)=\underset{i=1}{\overset{n\·m}{\Σ}}\left|\right|D_{i}x|{|}_2---\left(4\right)$

其中，D_i为局部有限差分算子，D_ix∈R²为x在第i个采样点处包含水平和垂直两个方向的一阶有限差分；TV正则化也有利于保留数据锋利的边缘和边界；地震数据的恢复为解决如下最小化问题：

$\underset{x}{min}TV\left(x\right)+{\mathrm{\μ\Φ}}_{fid}(x,f)---\left(5\right)$

即地震数据的重建模型为：

$\underset{x}{\min }\underset{i=1}{\overset{n\·m}{\Σ}}\left|\right|D_{i}x|{|}_2+\frac{\mathrm{\μ}}{2}||Kx-f|{|}_2^2---\left(6\right)$

其中，μ＞0为平衡二者的参数；x^*为原始地震数据；x为被恢复的地震数据；

对建立的地震数据恢复目标函数进行重写和转化具体为，将目标函数公式(6)重写为公式(7)形式，则有：

$\left\{\begin{array}{c}\underset{x,y}{\min }{\mathrm{\Σ}}_{i=1}^{n\·m}\left|\right|y_i|{|}_2+\frac{\mathrm{\μ}}{2}||Kx-f|{|}_2^2\\ \begin{array}{cc}s.t.& y_i=D_{i}x,i=1,...,n\·m\end{array}\end{array}\right.---\left(7\right)$

其中，y＝(y¹；y²)∈R^2nm，y¹和y²是长度为n·m维向量，满足((y¹)_i；(y²)_i)＝y_i∈R²；利用二次罚函数方法将公式(7)转化为无约束目标函数：

$\underset{x,y}{\min }\underset{i=1}{\overset{n\·m}{\Σ}}\left(\right|\left|y_i\right|{|}_2+\frac{\mathrm{\β}}{2}\left|\right|y_i-D_{i}x\left|{|}_2^2\right)+\frac{\mathrm{\μ}}{2}\left|\right|Kx-f|{|}_2^2---\left(8\right)$

其中，β＞＞0为罚参数；

利用交替方向法(ADM)对目标函数公式(8)进行迭代与交替方向最小化求解，初始x＝x^k，λ＝λ^k：

$\left\{\begin{array}{c}{y}^{k+1}\←\arg {\min }_y{L}_A(x^k,y,{\mathrm{\λ}}^k)\\ x^{k+1}\←\arg {\min }_x{L}_A(x,y^{k+1},{\mathrm{\λ}}^k)\\ {\mathrm{\λ}}^{k+1}\←{\mathrm{\λ}}^k-\mathrm{\β}(y^{k+1}-{\mathrm{Dx}}^{k+1})\end{array}\right.---\left(9\right)$

其中，L_A(x,y,λ)为公式(8)的增广Lagrangian函数：

$\begin{array}{c}{L}_A(x,y,\mathrm{\λ})=\underset{i=1}{\overset{n\·m}{\Σ}}\left(\left|\right|y_i|{|}_2-{\mathrm{\λ}}_i^T\left(y_i-D_{i}x\right)+\frac{\mathrm{\β}}{2}||y_i-D_{i}x|{|}_2^2\right)\\ +\frac{\mathrm{\μ}}{2}\left|\right|Kx-f|{|}_2^2\end{array}---\left(10\right)$

λ_i为R²空间中的向量；公式(9)中argmin_yL_A(x^k,y,λ^k)等价于如下形式：

$\underset{y_i\∈R^2}{\min }\left|\right|y_i|{|}_2+\frac{\mathrm{\β}}{2}||y_i-(D_i{x}^k+\frac{1}{\mathrm{\β}}{\left({\mathrm{\λ}}^k\right)}_i)|{|}_2^2,i=1,...,n\·m---\left(11\right)$

公式(11)通过二维收缩来解决：

$y_i^{k+1}=\max \left\{\left|\right|D_i{x}^k+\frac{1}{\mathrm{\β}}{\left({\mathrm{\λ}}^k\right)}_i|{|}_2-\frac{1}{\mathrm{\β}},0\right\}\frac{D_i{x}^k+\frac{1}{\mathrm{\β}}{\left({\mathrm{\λ}}^k\right)}_i}{\left|\right|D_i{x}^k+\frac{1}{\mathrm{\β}}{\left({\mathrm{\λ}}^k\right)}_i|{|}_2}---\left(12\right)$

在固定y^k+1和λ^k基础上，求解公式(9)中arg min_xL_A(x,y^k+1,λ^k)为一最小二乘问题，其对应的标准方程为：

$(D^{T}D+\frac{\mathrm{\μ}}{\mathrm{\β}}{K}^{T}K)x=D^T(y^{k+1}-\frac{1}{\mathrm{\β}}{\mathrm{\λ}}^k)+\frac{\mathrm{\μ}}{\mathrm{\β}}{K}^{T}f---\left(13\right)$

其中，D＝(D¹；D²)∈R^2nm×nm；

以上各个数学符号的物理意义

m：地震道记录的采样点个数；

n：地震道记录个数；

r：各地震道与每个地震道采样点组成的n×m矩阵；

x_ij：r中的元素，指第i个地震道、第j个采样点值，其中，i＝1,2,…,n；j＝1,2,…,m；

vec(r)：矩阵r的拉伸向量；

x：被恢复的地震数据；

K：褶积算子；

f：观测的地震数据；

R^nm：n·m维度空间；

w：添加的随机噪声；

||·||₂：指L₂范数；

Φ_fid(x,f)：随机噪声的数据保真项；

D：一阶全局有限差分算子；

D¹,D²：分别为水平和垂直方向的一阶全局有限差分算子；

D_i：局部有限差分算子，以D¹的第i行为第一行、以D²的第i行为第二行的两行矩阵；

TV(x)：全变分的离散形式；

μ：平衡数据保真项Φ_fid(x,f)和全变分形式TV(x)的参数；

y_i＝D_ix＝((y¹)_i；(y²)_i)：x在第i个采样点处包含水平和垂直两个方向的一阶有限差分，是一2维向量；

y＝(y¹；y²)：2nm维向量；

y¹,y²：长度为n·m维的向量，分别指水平和垂直两个方向的一阶有限差分；

β：二次罚函数方法中的罚参数；

λ：Lagrangian函数参数；

λ_i：R²空间中的向量；

x^k,y^k,λ^k：分别是x,y,λ的第k次迭代更新值；

L_A(x,y,λ)：增广Lagrangian函数。