稀疏反褶积正则化策略优选：数学模型构建与性能评估-凯发娱乐官网

1. 前言

1.1. 研究背景及意义

地震勘探是探测地下资源(如石油、天然气和矿产)的关键技术，通过分析地震波的传播和反射特性，可以绘制地下结构的详细图像。然而，地震数据在采集和处理过程中受到地震子波复杂性和随机噪声的影响，导致反射波难以直观解释，限制了勘探的分辨率和准确性[1] [2]。

随着油气勘探的深入，小尺度岩性油气藏的识别变得更加重要。在褶积模型下，需从地震信号中准确反褶积出反射系数信息，包括其位置和振幅。有限频地震信号的反褶积常常面临多解问题，实际应用中地质模型的成层性要求解具有稀疏性[3] [4]。

传统反褶积方法，如基于wiener的最小均方误差准则，虽然有效但在处理异常值和保持数据稀疏性方面存在不足。为此，研究人员提出了基于l1正则化的稀疏反褶积方法，通过l1范数增强反射系数序列的稀疏性。但l1正则化在噪声和异常值面前仍显敏感，进而出现了l1/2正则化方法，但其仍有局限。最近几年来，cauchy正则化逐渐受到关注，它通过对异常值进行温和惩罚，有效提高了反褶积结果的稳定性和鲁棒性[2] [5]-[8]。

本文重点比较了cauchy正则化与l1和l2联合正则化下的稀疏反褶积，旨在分析不同方法在稀疏反褶积中的效果，并提供理论支持和实用指导。

1.2. 本文结构安排

本文结构如下：首先介绍地震数据反褶积的稀疏性概念及常用正则化技术的理论基础，如l1范数、l2范数和cauchy正则化。接着探讨cauchy分布正则化及l1和l2联合正则化的方法，包括其数学模型和求解过程。然后，通过模拟和实际地震数据实验，比较这些方法在稀疏反褶积中的特点与优势。最后，讨论各正则化方法的适用性、实际应用中的问题、参数选择对结果的影响以及计算效率。

2. 国内外研究现状

反褶积理论最早由enders robinson提出，主要用于提高地震资料的垂向分辨率。然而，随着地震勘探需求的不断提高，传统的线性维纳滤波方法逐渐难以满足薄层勘探的精度要求。为此，学者们在过去的四十多年里一直致力于研究稀疏约束的反褶积方法[6] [9]-[12]。

wiggins首次提出了最小熵反褶积(med)，这是应用稀疏准则进行反褶积的重要开端。随后，gray在稀疏性假设的基础上，提出了变量范数约束的反褶积方法。youzwishen则进一步改进了cauchy准则，并将其成功应用于贝叶斯反演框架中。canadas提出了盲反褶积理论，使得反演反射系数序列和地震子波成为可能。之后，meng和zhang分别利用cauchy准则和改进的cauchy分布，进行了盲反褶积的求解研究。danilo指出，稀疏脉冲反褶积的目的是从噪声记录中提取有效脉冲。liu通过cauchy分布约束，改进了稀疏反射序列的求解方法[13]-[16]。

wang提出了一种基于geman范数约束的频率域反射系数反演方法，并且基于toeplitz稀疏矩阵分解进一步发展了盲反褶积技术。kang等人引入l1/2正则化理论，推进了地震稀疏反褶积的发展。ni等通过引入基于cnn的深度学习子波整形反褶积方法，大大提升了反褶积处理的效果[16]-[18]。

3. 基于正则化约束的反褶积方法

3.1. 最小平方反褶积

由震源激发的纵(横)波经地下传播并被人们在地面或井中接收到的地震波，通常是一个有一定长度的脉冲振动，它描述了介质质点的振动规律，应用信号分析领域中的广义术语，可称为振动信号，在地球物理领域称为地震子波[2]。

在自激自收条件下，robison提出利用褶积模型描述地震响应，地震记录可近似看作是地震子波与地层反射系数的褶积结果，又因为野外实际地震数据采集时不可避免地会含有噪声，故地震记录可表示为：

s( t )=w( t )∗r( t ) n( t ) ，

式中s( t ) 表示地震记录，w( t ) 表示子波，r( t ) 表示反射系数序列，n( t ) 表示随机噪声，*表示卷积操作。

地震子波与反射系数褶积生成无噪地震数据的过程亦可用图1表示。

figure 1. convolution model of earthquake records: (a) earthquake wavelet, (b) reflection coefficient, (c) earthquake record

图1. 地震记录褶积模型：(a) 地震子波，(b) 反射系数，(c) 地震记录

我们考虑其中一道反射系数序列，变为如下的矩阵运算形式：

d=wr n ，

式中d 为地震道列向量，w 为子波褶积矩阵，r 为反射系数列向量，n 为随机噪声列向量。

由于地震子波的影响，来自相邻界面的反射波重叠在一起，难以区分开来，因此需要压缩子波来提高地震勘探资料的纵向分辨率。传统的反褶积方法是设计一个与子波性质相反的滤波因子，然后与地震记录褶积，得到期望的窄脉冲，这个过程称为反褶积。若假设地震道长度为k ，则：

[ d 1 d 2 ⋮ d k−1 d k ]=[ ω 0 ω −1 ⋯ ⋯ ω 1−k ω 1 ω 0 ω −1 ⋱ ⋮ ⋮ ω 1 ω 0 ω − 1 ⋮ ⋮ ⋱ ω 1 ω 0 ω −1 ω k ⋯ ⋯ ω 1 ω 0 ][ r 1 r 2 ⋮ r k−1 r k ] [ n 1 n 2 ⋮ n k−1 n k ] 。

从反演的角度考虑，在子波褶积矩阵w 已知的条件下(可以通过测井数据或地震记录中提取)我们希望找到一个模型r ，使得正演的记录与实际地震记录最小平方误差最小，即寻找使得j 最小的r 。

j= min r 1 2 d−w r 2 2 ，

对上式中的r 进行求求偏导，并令其导数为0：

∂j ∂ r i =0 ，

上式可变形化简为：

w t ( w t r−d )=0 ，

即：

r= ( w t w ) −1 w t d ，

式中w t w 称为协方差矩阵，( w t w ) −1 w t 称为w 的最小平方逆，也称为广义线性逆。在地球物理问题中，虽然协方差矩阵通常是可以求出来的，但它可能并不唯一存在，因此上式为一个病态问题，无法得到适定的、唯一的解。故需要对模型r 进行正则化约束。

3.2. 基于不同范数正则化约束的反褶积方法

对于最小平方反褶积模型，反演目标函数可修改为：

j= min r 1 2 ‖ d−wr ‖ 2 2 h( r ) ，

式中的h( r ) 为正则化项，通过其对模型r 的解一些约束。

3.2.1. 基于l1范数正则化约束的反褶积方法

正则化项的选择通常取决于反射系数所满足的先验条件。一般认为反射系数序列由较大的稀疏脉冲反射系数和服从高斯分布的微小反射系数背景组成。如果只是想重点突出较大的反射系数，忽略微小的高斯背景，就可以认为反射系数序列具有稀疏性。

l1正则化是经典的稀疏求解方法，其核心思想是通过引入l1范数惩罚项，促使解中的大部分系数趋向于零，从而实现稀疏性。具体来说，l1正则化倾向于对解的非零元素进行惩罚，而较小的系数则倾向于被压缩为零。这种特性使得l1正则化在稀疏信号重构领域表现优异，我们用l1范数衡量稀疏性，所以反演目标函数中的正则化项可设置为：

h( r )= α 1 ‖ r ‖ 1 ，

其中α 1 为正则化因子参数，用来调节惩罚项的比重。因此，目标函数修改为：

j= min r 1 2 ‖ d−wr ‖ 2 2 α 1 ‖ r ‖ 1 。

3.2.2. 基于l2范数正则化约束的反褶积方法

除了l1范数正则化外，l2范数正则化也被广泛应用于反褶积问题中。l2正则化通过加入l2范数惩罚项，促进解的平滑性。这种正则化倾向于对解中所有的系数施加均匀的惩罚，从而避免某些系数过大。与l1正则化相比，l2正则化不会压缩系数到零，而是趋向于减小系数的幅度，导致解的“稀疏性”不明显，但其平滑性更好。因此，l2正则化经常用于那些不需要特别强调稀疏性的场景，如图像处理中的平滑处理或者物理模型的稳健解。

在这种情况下，正则化项被设置为l2范数形式，具体表示为：

h( r )= α 2 ‖ r ‖ 2 2 ，

其中，α 2 是正则化因子，用于控制l2正则化项对目标函数的影响程度。将此正则化项代入反演目标函数后，修改后的目标函数为：

j= min r 1 2 ‖ d−wr ‖ 2 2 α 2 ‖ r ‖ 2 2 。

3.2.3. 基于cauchy分布正则化约束的反褶积方法

cauchy分布正则化是一种介于l1和l2正则化之间的方法，适用于反射系数序列既含有较大脉冲又有较小高斯噪声背景的情况。cauchy正则化源于cauchy分布的概率模型。cauchy分布是一种具有重尾特性的概率分布，这意味着它能更好地捕捉那些具有大幅度变化的信号分量。与l1和l2正则化不同，cauchy正则化的目标函数对大幅度的解施加较少的惩罚，从而保留信号中的重要信息，同时对较小的噪声分量有较强的压制作用，从而实现更精确的稀疏解。cauchy正则化项通常以以下形式表示：

h( r )= ‖ r ‖ cauchy �= ∑ i=1 k ln( 1 r i 2 β 2 ) ，

将cauchy正则化项引入到目标函数中，得到以下优化问题：

j= min r 1 2 ‖ d−wr ‖ 2 2 α 3 ‖ r ‖ cauchy 。

其中，α 3 是正则化因子，β 是一个尺度参数，控制cauchy分布的宽度。与l1和l2正则化不同，cauchy正则化项通过对反射系数平方的对数函数进行建模，使得大的反射系数受到较小的惩罚，而小的反射系数受到较大的抑制，从而保持了稀疏性。

3.2.4. 基于l1和l2联合范数正则化约束的反褶积方法

单独使用l1范数或l2范数正则化各有其优势和局限性。l1范数正则化适合突出稀疏的反射系数，但可能对噪声不够敏感，而l2范数正则化能够有效抑制高频噪声和保持解的平滑性，但会导致解不够稀疏。因此，结合3.2.1和3.2.2中讨论的两种范数的特点，同时利用l1正则化的稀疏性和l2正则化的平滑性，平衡解的稀疏性与其抗噪能力，采用l1和l2联合范数的正则化方法(也称为弹性网正则化，elastic net)可以结合两者的优势，提高反褶积结果的质量。联合正则化方法的目标是通过引入一个同时包含l1和l2范数的正则化项，既保持反射系数的稀疏性，又能够有效抑制噪声，获得更为平滑且精确的解。其正则化项可以表示为：

h( r )= α 1 ‖ r ‖ 1 α 2 ‖ r ‖ 2 2 ，

其中α 1 和α 2 为正则化因子参数，用来调节惩罚项的比重。因此，目标函数修改为：

j= min r 1 2 ‖ d−wr ‖ 2 2 α 1 ‖ r ‖ 1 α 2 ‖ r ‖ 2 2 。

3.3. 对于不同正则化反褶积的求解

本节深入探讨基于l1/2正则化和cauchy正则化的反褶积求解方法，旨在增强反射系数序列的稀疏性恢复，同时有效抑制噪声影响。

3.3.1. 基于交替方向乘子法的cauchy分布正则化约束的反褶积方法求解

cauchy正则化方法因其在处理异常值时的温和惩罚特性，逐渐在地震反褶积中获得应用。cauchy正则化项的引入使得目标函数变得非凸，因此求解过程通常采用交替方向乘子法(alternating direction method of multipliers, admm)。admm是一种适用于大规模优化问题的有效算法，通过将复杂的优化问题分解为多个简单的子问题，从而简化求解过程。具体的求解步骤如下：

1) 目标函数分解：将目标函数分解为与反射系数序列和正则化项相关的独立子问题；

2) 引入拉格朗日乘子：通过增广拉格朗日函数分离线性项和非线性项，引入辅助变量，转化为多个拉格朗日子问题的求解；

3) 逐步更新和迭代：交替更新反射系数序列和拉格朗日乘子，直至收敛。

admm具有良好的分布式计算特性和稳定性，特别适用于大规模数据集的非凸问题。其交替更新策略有助于在保持稀疏性的同时抑制噪声，提升反褶积的准确性。

3.3.2. 基于迭代阈值算法的l1和l2联合范数正则化约束的反褶积方法解法

l1/2正则化结合了l1和l2正则化的优势，适合稀疏信号恢复问题。通常采用迭代阈值算法(iterative threshold algorithm, ita)来求解。具体步骤如下：

1) 初始化：设定初始反射系数序列，选择适当的阈值参数和正则化因子；

2) 计算残差和阈值操作：每次迭代中，计算残差并应用迭代软阈值公式更新反射系数，以保证稀疏性和精确性；

3) 更新反射系数序列：根据残差和阈值操作的结果，更新反射系数，调整正则化因子以平衡稀疏性和解的平滑性；

4) 迭代至收敛：继续迭代直至满足收敛条件。

ita方法计算简单、收敛快速，特别适用于大规模地震数据处理。通过l1/2正则化，可以有效突出稀疏特征，提升反褶积的分辨率和稳定性。

4. 模型试验与结果分析

4.1. 实验数据和参数设置

本文假设子波已知，实验中模型的子波和频率为5~10~60~80 hz的ormsby小波(图2)。在实际处理中，子波可以从地震记录中通过盲提取获得，如利用高阶统计量的方法[19]。

figure 2. ormsby wavelet in time domain (red) and frequency domain (blue) with a frequency range of 5~10~60~80 hz

图2. 频率为5~10~60~80 hz的ormsby小波时域图(红色)和频域图(蓝色)

4.2. 单道地震信号实验结果

我们取4.1中提到的ormsby子波，生成一个模拟的地震信号与其进行褶积操作，并对其加入了标准差为0.15的高斯噪声，得到的褶积结果作为一个模拟的原始地震记录(如图3)，并对其进行3.2.3和3.2.4中提到的方法进行反褶积操作，求出对应的反射系数序列，最后将反射系数序列和子波褶积操作(重构地震信号)，分别对应图4和图5。

figure 3. wavelet (top) and noisy seismic record (bottom)

图3. 子波(上)与含噪地震记录(下)

figure 4. comparison between the reflection coefficient obtained using the cauchy method and the true reflection coefficient (top); comparison between the reconstructed seismic signal and the original seismic signal (bottom)

图4. 使用cauchy方法求得的反射系数结果与真实反射系数的对比(上)；重构地震信号与原始地震信号的对比(下)

figure 5. comparison between the reflection coefficient obtained using the l1/2 method and the true reflection coefficient (top); comparison between the reconstructed seismic signal and the original seismic signal (bottom)

图5. 使用l1/2方法求得的反射系数结果与真实反射系数的对比(上)；重构地震信号与原始地震信号的对比(下)

我们可以直观地看到两种方法分别求得的反射系数都比较贴合原始的反射系数，都在恢复稀疏序列和抗噪性方面有较好的表现，但在细节方面l1/2的方法明显好于cauchy的方法。对于一维简单模型我们暂时只做一些直观的观察，对与下面章节中的多道模型，加入定量分析。

4.3. 多道地震记录实验结果

在本节中，我们对多道地震记录进行了反褶积实验，评估提出的反褶积方法在不同地震道上的效果。为了确保实验结果的可比性和一致性，我们对每一列地震道分别进行反褶积处理，并记录反射系数序列和重建的地震记录。

我们基于marmousi2模型中的p波速度模型的局部区域(图6中的红色方框)，构建了一个地震信号模型，如图7(左)所示，图7(右)则展示了加入标准差为0.05的高斯噪声的地震记录，共301道，时间采样间隔为1 ms。

图8分别展示了使用cauchy方法对于含噪声的地震记录的反褶积结果(左)和重构的地震记录(右)，图9则展示了l1/2方法的结果。从红色框中可以直观地看出，后者方法对于较小反射系数的保护优于前者，且对于低频噪声的压制优于前者；从整体来看，后者的处理结果纹理更加清晰，而前者的结果有一些过度锐化的痕迹。

figure 6. marmousi2 p wave velocity model

图6. marmousi2 p波速度模型

4.4. 多道地震记录实验结果分析

为了进一步量化对比不同方法的效果，我们使用均方误差和相关系数作为评价标准，对含噪地震记录的反演结果(反射系数和重构地震记录)，分别对两种方法的结果进行计算。

figure 7. reflection coefficient model (left) and model with noise (right)

图7. 反射系数模型(左)及含噪声的模型(右)

figure 8. reflection coefficient results using cauchy regularization constraints (left) and reconstructed seismic records (right)

图8. 使用cauchy正则化约束的反射系数结果(左)和重构的地震记录(右)

均方误差(mean squared error)：计算重建地震记录与原始地震记录之间的均方误差。均方误差越小说明拟合效果越好。

相关系数(correlation coefficient)：评估计算得出的反射系数序列与真实反射系数序列的相关性。高相关系数表明恢复的反射系数与真实反射系数之间的一致性更高。

信噪比(signal-to-noise ratio)：信噪比是用于评估信号质量的关键指标，它定义为信号强度与噪声强度的比值。具体来说，信噪比越大，说明信号在噪声中的占比越高，恢复的反射系数和重构地震记录的质量越好。

figure 9. reflection coefficient results (left) and reconstructed seismic records (right) using l1/2 joint regularization constraints

图9. 使用l1/2联合正则化约束的反射系数结果(左)和重构的地震记录(右)

对于不同方法的以上评价指标计算结果如表1所示。

table 1. comparison of mean square error and structural similarity of deconvolution results of the same model by different methods

表1. 不同方法对同一模型做反褶积操作结果的均方误差和结构相似性对比

从表中可以看出，针对上述的marmousi2模型，基于l1和l2范数联合约束的方法的均方误差比基于cauchy约束的方法要低，并且相关系数和信噪比也更高。

实验结果表明，基于l1和l2范数联合约束的反褶积方法相对于使用cauchy正则化的方法，在处理地震数据时，能够更有效地压缩子波，从而提高地震数据的分辨率，并且在存在噪声的情况下相对于后者仍然能够保持较好的稳定性和鲁棒性。上述结果说明了基于l1和l2范数联合约束的反褶积方法相对于使用cauchy正则化的方法表现出更明显的优越性。

5. 结论

5.1. 研究总结

本研究探索了地震反演中不同正则化方法的效果的对比，通过比较l1和l2范数正则化联合约束与cauchy正则化约束在提升地震反演模型的精度方面的表现得出了一些结论。研究结果表明，l1和l2范数正则化联合约束在提高反演结果的分辨率和稳定性方面优于单独使用cauchy约束。这一结论在对多道数据和复杂模型(如marmousi2模型)的测试中得到了验证。通过引入l1和l2联合正则化，反演模型能够更有效地抑制噪声，同时保持地下结构的真实特征。

l1和l2联合正则化方法通过将l1正则化的稀疏性与l2正则化的平滑性相结合，成功地在稀疏求解问题中找到了平衡。这种结合方法不仅能够提高解的稀疏性，减少噪声影响，还能通过平滑化提升解的稳定性和抗噪能力。因此，l1和l2联合正则化在地震反褶积、图像处理、机器学习等多个领域中都显示出了广泛的应用前景，尤其适用于需要同时解决稀疏性和抗噪性的复杂问题。

5.2. 不足与局限性

尽管本研究表明了l1和l2联合正则化的优势，但仍存在一些不足与局限性。首先，本研究的结果主要基于对marmousi2模型和多道数据的测试，对于其他类型的地质模型和单道数据，效果有待进一步验证。其次，在处理极高噪声数据时，正则化参数的选择变得尤为敏感，这可能会影响反演结果的稳定性。此外，模型计算复杂度的增加也可能在实际应用中带来额外的计算开销。在模型适用性方面，本文方法假设反射系数序列具有稀疏性，但在一些地质条件复杂的实际应用中，该假设可能不完全成立，影响方法的适用范围和效果。

5.3. 未来展望

未来的研究可以在以下几个方面进行改进和拓展。首先，可以研究更为自动化的正则化参数选取方法，例如使用交叉验证或贝叶斯优化技术来动态调整l1和l2联合正则化参数，从而提高反演模型的适应性和稳定性。其次，将本研究的l1和l2联合正则化方法应用于其他复杂地质模型(如盐丘模型)和不同类型的数据集(如实际野外勘探数据)，以检验其广泛适用性和鲁棒性。此外，可以考虑将正则化方法与机器学习技术相结合，利用深度神经网络进行特征提取，从而进一步提高反演精度和计算效率。

通过本文的研究，进一步丰富了反褶积方法的理论和实践，为提高地震勘探资料的分辨率和精度提供了新的思路和方法。希望未来能够在实际应用中发挥更大的作用，促进地球物理勘探技术的发展。

参考文献