基于交叠组稀疏广义全变分的地震信号随机噪声衰减

陈颖频 彭真明 李美惠 喻 飞

(①电子科技大学信息与通信工程学院，四川成都 610054；②闽南师范大学物理与信息工程学院，福建漳州 363000；③电子科技大学信息地学研究中心，四川成都 610054)

0 引言

实际地震信号中存在大量随机噪声，这给地震资料的处理和解释带来较大的干扰。目前，去除地震信号随机噪声的方法主要有基于全变分模型的正则化去噪方法[1-4]、基于稀疏表示和字典学习的随机噪声去除方法[5-7]和基于变换域的滤波去噪方法[8-9]等。其中，基于稀疏表示和字典学习的去噪方法需观察数据进行字典训练，在获得较好去噪效果的同时，要耗费大量时间在字典学习上，只要数据量稍大就会引起较大计算量。基于变换域的滤波去噪方法通常需将地震信号从空域映射到变换域中进行截断滤波，易产生吉布斯效应。

全变分(Total variation，TV)模型被证明是一种可有效去除随机噪声的模型，自从它被Rudin等[1]提出之后，就引起学者们的广泛关注[10-16]。TV模型充分挖掘了二维图像的横向纵向梯度信息，较好地契合了自然图像的局部光滑和梯度稀疏等先验知识。该模型被广泛应用于地震图像去噪[10-11]、地震波阻抗反演[12-13]、地震波全波形反演[14]、弱小目标检测[15]和超分辨率分析[16]等众多领域。

由于TV模型假设图像是分片光滑常数，导致TV去噪模型存在较为严重的阶梯效应。为抑制阶梯效应和提高TV模型的去噪效果，Bredies等[17]提出广义全变分(Total generalized variation，TGV)模型。该模型是全变分模型的推广，具有凸性、下半连续性、旋转不变性等众多优秀的数学性质，并能逼近任意多项式[18]。TGV模型同时约束了图像的一阶梯度与二阶梯度，因此有效缓解了全变分模型的阶梯效应。Knoll等[18]将广义全变分模型用于核磁共振成像，取得较好的应用效果。随后，Guo等[19]提出一种基于TGV和剪切波变换的细节保留方法，用于MRI图像重构中。Qin等[20]将这种多约束正则项用于图像解卷积。Kong等[11]则将Shearlet和TGV正则项用于地震信号去噪。

近年来，Liu等[21]和Chen等[22-23]提出了交叠组稀疏正则项。该正则项是一种非分离正则项[24]。它不仅考虑了图像差分域的稀疏性，还挖掘了每个点邻域差分信息，因此可利用图像梯度的结构化稀疏特性。通过交叠组合梯度可提高平滑区域与边界区域的差异，从而抑制TV模型的阶梯效应。Liu等[25]借鉴Chen等[22-23]的成果，将一维交叠组稀疏正则项推广到二维领域，并将其引入各向异性全变分模型，用于椒盐噪声的去噪和解卷积问题。Liu等[26]将交叠组稀疏正则项用于Speckle噪声的去除。

上述工作全部是基于各向异性全变分模型做的推广，然而各向异性全变分模型没有二阶差分的约束，交叠组稀疏梯度的引入对阶梯效应的抑制能力有限。需指出的是，TGV和交叠组稀疏收敛算子对阶梯效应的抑制机理并不一样，前者利用图像一阶、二阶差分约束缓解阶梯效应，后者则通过图像梯度的结构特性抑制阶梯效应。目前，TGV模型与交叠组稀疏收敛算子的交叉研究仍处于起步阶段。

由于经典的TGV模型未考虑图像梯度的邻域结构特性，受Liu等[25]的启发，本文提出一种基于交叠组稀疏收敛算子的改进广义全变分模型(简称为TGV-OGS)，并将其应用于地震信号去噪中。鉴于该模型的复杂性，采用交替乘子迭代法(Alternating direction method of multipliers，ADMM)[27-28]对其求解。为提高算法效率，假定处理图像满足周期性边界条件，并采用快速解卷积方法[29-31]，巧妙地避免了大型差分矩阵的相乘运算，将差分操作理解为卷积，再利用卷积定理，从频域恢复图像。

在后续实验中，将对比常规各向异性全变分方法(Anisotropic total variation，ATV)[32]、各向异性组稀疏方法(Anisotropic total variation based on overlapping group sparsity，ATV-OGS)[25]、TGV去噪方法及本文去噪方法，并从峰值信噪比(Peak signal to noise ratio，PSNR)、结构相似性(Structural similarity，SSIM)[33]、信噪比(Signal to noise ratio，SNR)[7]及运算时间等指标客观地对比这几种算法。从实验结果可见，本文提出的算法及所构建的模型能进一步改进TGV的去噪性能。

1 预备知识

经典的二阶TGV模型[34]定义为

‖Dv*S-Vv‖1)+α1(‖Dh*Vh‖1+

‖Dv*Vv‖1+‖Dv*Vh+Dh*Vv‖1)]}

(1)

1.2 二维矩阵交叠组稀疏邻近算子与MM算法

假定V0为待收敛矩阵，则其交叠组稀疏邻近算子记为

P(V) =proxγφ(V0)

(2)

(3)

图1 二维交叠组稀疏示意图

利用优化—最小化方法(Majorization-minimization，MM)[36]可有效计算式(2)。根据MM算法，最小化P(V)需首先找到一个函数Q(V,U)(U为辅助变量)，该函数对所有V,U，都有Q(V,U)≥P(V)，并且当且仅当U=V时等号成立。据此，每次计算的Q(V,U)最小值为P(V)的优化值，式(2)的计算可转化为如下迭代运算

(4)

考虑到组稀疏正则项的特殊性，并注意到下列不等式的存在

(5)

式中，当且仅当U=V时，等号成立。

(6)

将S(V,U)重写为

(7)

式中：v是矩阵V的向量形式；C(U)与V无关，可视为关于V的常数项；D(U)∈RN2×N2，是一个对角矩阵，其对角元素定义为

[D(U)]m,m

(8)

结合式(4)、式(7)和式(8)，可将式(2)转化为如下迭代优化问题

(9)

该式为二次规划问题，其迭代最优解为

V(k+1)=mat{[I+γD2(V(k))]-1v0}

(10)

式中：I∈RN2×N2，表示单位矩阵；v0是V0的向量形式； mat表示向量矩阵化运算。

2 模型构建及其求解方法

2.1 基于交叠组稀疏的改进TGV模型

将TGV中的L1约束项改进为交叠组稀疏约束项，从而更好地挖掘图像一阶梯度与二阶梯度的差分信息，建模如下

φ(Dv*S-Vv)]+α1[φ(Dh*Vh)+

φ(Dv*Vv)+φ(Dv*Vh+Dh*Vv)]}}

(11)

对比经典的TGV模型，本文构建的模型将一阶梯度与二阶梯度矩阵的每个像素点邻域K2个信息点交叠组合，从而更充分挖掘一阶梯度矩阵和二阶梯度矩阵结构特性。显然，每个像素点的邻域梯度与二阶邻域梯度存在一定相似性，该结构稀疏先验知识被交叠组稀疏模型较好地刻画。

2.2 构建模型的ADMM求解

为求解式(11)定义的改进TGV模型，利用ADMM框架求解模型。令

则根据ADMM原理，可将式(11)的增广拉格朗日函数表达为

μ0[φ(Z1)+φ(Z2)]+μ1[φ(Z3)+

φ(Z4)+φ(Z5)]-

〈β0Λ1,Z1-(Dh*S-Vh)〉-

〈β0Λ2,Z2-(Dv*S-Vv)〉-

〈β1Λ3,Z3-(Dh*Vh〉-〈β1Λ4,Z4-Dv*Vv〉-

〈β1Λ5,Z5-(Dv*Vh+Dh*Vv)〉+

(12)

式中：μ0=μα0；μ1=μα1；Λi(i=1,2,…,5)是拉格朗日乘子； 〈X,Y〉表示两个矩阵X与Y的内积。

在ADMM框架下，分离变量及其对偶变量之间与三元组S、Vh、Vv是去耦合的。而三元组中三个变量是相互耦合的，因此需要建立关于三个变量的三元一次方程组。对于S子问题，其子目标函数为

(13)

本文假定处理数据满足周期性边界条件，因此可以利用快速傅里叶变换进行卷积计算。

根据卷积定律，两个矩阵在空域卷积，对应到频域，卷积结果的频谱为两个矩阵频谱的点乘。将式(13)做傅里叶变换

(14)

(15)

整理可得

(16)

其中

式中：1是元素全为1的矩阵；上角“*”表示共轭运算。

同理，对Vh子目标函数关于Vh求导为零，可得Vh子问题的约束方程为

(17)

其中

详细推导过程不再赘述。

同理，对Vv的子目标函数关于Vv求导为零，可得Vv子问题的约束方程为

(18)

其中

(19)

该式可用克莱姆法则和快速反傅里叶变换求解

(20)

式中：除法都为按元素点除运算；F-1表示二维逆快速傅里叶变换算子；定义矩阵行列式算子为

K12。K23。K31+K13。K21。K32-K13。K22。K31-

K12。K21。K33-K11。K32。K23

对于Z1子问题，其目标子函数为

(21)

根据式(10)，可得Z1的更新公式为

(22)

同理，Z2、Z3、Z4和Z5的更新公式为

(23)

且有

式中vec(·)表示将矩阵列向量化。

Λ1的目标子函数为

(24)

利用梯度上升法可得其更新公式

(25)

其中γ为学习率。

类似地，拉格朗日乘子变量Λ2、Λ3、Λ4和Λ5的更新公式对应为

(26)

至此，本文构建模型的所有子问题都得以解决。将该算法总结为表1。

表1 TGV-OGS算法伪代码

其中： tol=10-4表示算法停止迭代的阈值；Max=30

3 数值实验

将本文算法应用于地震信号，并通过峰值信噪比、结构相似性指示系数、信噪比及运行时间客观评价该算法的去噪性能，并与ATV、ATV-OGS方法以及TGV去噪方法做全面对比。

3.1 评价指标

去噪领域最常用评价指标有：峰值信噪比、结构相似性指标、信噪比和运行时间。其中峰值信噪比、结构相似性指标[33]为

(27)

SSIM(X,Y)

(28)

为更好地表征去噪效果，还选用了信噪比参数，定义[7]如下

(29)

3.2 合成地震信号去噪性能对比测试

本节给出合成地震记录，该记录由20 Hz雷克子波与人工合成反射系数卷积获得。图2为截取的前50道合成记录，可清晰观察到去噪效果。

图2 合成地震信号

3.2.1 交叠组稀疏正则项效果测试

为验证交叠组稀疏正则项的有效性，在合成记录中加入0均值、标准差为σ的高斯白噪声，且设定K=3。ATV、ATV-OGS、TGV以及本文方法的测试结果如表2所示，最优指标用粗体标出。

表2 合成地震信号算法性能测试表

从表2中可看到，交叠组稀疏正则项不仅对经典的各向异性TV模型算法性能有所改进，而且对TGV去噪方法有较好的改进。为直观反映其效果，以伪彩色形式(图3)展示σ=30时上述四种算法的去噪效果。

观察图3中由矩形红色框圈定区域(图3d)，显然是平滑区域，在放大区域中可以看到若干个被重噪声污染的像素；ATV去噪方法在该区域有较明显的块效应。从放大图(图3f)中可见，仍然有两个像素被重噪声干扰。从图3h可见，利用交叠组稀疏正则项则非常好地抑制了块效应，且两个重噪声点被去除得较为干净。值得指出的是，各向异性交叠组稀疏方法并未完全抑制ATV的块效应。再观察TGV的去噪结果(图3i、图3j)，可见TGV去噪方法也对ATV的块效应有较好的抑制作用，在平滑区域获得更加光滑的重构信号。但从图3j可见，TGV去噪方法对于重噪声污染点也是无能为力的。而本文提出的方法则将交叠组稀疏和TGV的优点结合起来，获得比ATV-OGS和TGV更好的去噪效果。观察图3l可以看到，本文算法获得的去噪结果在该区域更接近原图。值得注意的是，设定K=3时，算法性能尚未达到最优。

为进一步观察去噪效果，将图3中六个子图的第40道抽出进行观察(图4)。从图4可见ATV存在明显的阶梯效应，而ATV-OGS去噪方法(图4d)较好地缓解了ATV的阶梯效应，但是去噪结果仍然不够光滑。传统TGV方法也存在类似问题。本文提出的方法在TGV理论框架上结合了组稀疏收敛方法，从而挖掘了数据邻域梯度的结构信息，将两者的优点充分结合，大幅提高了信号的重构质量。在图4b中50～100ms区间，地震信号被一强幅度噪声严重干扰。该强噪声污染在图4c和图4e中依然存在，这反映传统ATV和TGV方法的重大缺陷，即两种算法都是以单一像素点为处理对象，孤立地迭代，这样就无法充分考虑到信号邻域的梯度相似信息。而ATV-OGS方法和本文提出方法则不存在这种问题，观察图4d和图4f可知，基于交叠组稀疏收敛方法的ATV-OGS和TGV-OGS都能有效去除50～100ms存在的尖峰干扰。显然，本文提出方法充分利用了地震信号的邻域相似性，这对去除大幅度的噪声污染尤其有效。通过上述实验，得出如下结论，各向异性交叠组稀疏去噪方法能有效利用图像一阶梯度的结构特性，而广义全变分模型中存在一阶和二阶图像梯度，这些梯度同样含有邻似性，将交叠组稀疏正则项与TGV模型充分结合，进一步提升了TGV模型的去噪能力。

图3 四种算法去噪效果对比

(a)原图；(b)图a的局部放大；(c)带噪声图；(d)图c的局部放大；(e)ATV去噪结果；(f)图e的局部放大；(g)ATV-OGS去噪

结果；(h)图g的局部放大图；(i)TGV去噪结果；(j)图i的局部放大图；(k)TGV-OGV(K=3)去噪结果；(l)图k的局部放大

3.2.2 参数敏感性分析

对本文算法的交叠组合数K进行测试对比，以评价其对算法的整体效果。采用PSNR和SSIM两个指标对算法进行客观评价，针对不同高斯噪声水平，将K从1到13连续变化，并记录PSNR和SSIM值，详见图5。从图5可以看到，随着K的增大，在不同噪声下，K对PSNR和SSIM起到的作用是不一样，显然，在低噪声的时候(例如σ=10)，K取11时，PSNR和SSIM都达到峰值，若K继续加大，PSNR则下降。当σ=10时，本文算法的PSNR达到最高值37.7846 dB，高出TGV算法1.905 dB。从图5可见，当噪声较低的时候，图像的邻域信息对算法性能起到了正面作用。然而，K也不宜取得过大，否则可能会取到邻域中图像变化剧烈的区域，这种情况下，PSNR和SSIM就要下降。而当噪声较大时，邻域梯度的结构特性被破坏得比较严重，这种情况下，K稍大就会导致算法性能下降，例如当σ=40时，K=7时获得了PSNR和SSIM的最高值。从图5a中可见，当标准差为30时，设定K=9可以令PSNR值达到最高。图6展示了在标准差为30时，

图4 四种算法单道去噪效果对比

图5 不同σ值的PSNR(a)和SSIM(b)随K值变化曲线的对比

图6 TGV-OGV(K=9)去噪效果图

算法参数K=9的去噪效果，其PSNR为30.4824dB，SSIM为0.9789，SNR为15.8382dB，均远超过经典TGV去噪模型。对比图6b与图4l可见，当K设置得足够合理，恢复出来的地震信号将更加符合梯度的稀疏先验知识。再对比图6c与图4f可知，当K取值合适时，本文算法的保护边缘能力及对重噪声的抗噪能力更强。

3.3 实际地震信号去噪性能测试

以二维地震信号作为测试信号(图7a)。该地震信号为叠后地震信号，已经去除了噪声。为客观评价算法性能，在地震信号中加入σ=10，20，30，40的高斯白噪声，对比ATV、ATV-OGS、TGV以及TGV-OGS等算法的去噪性能指标，各项指标被记录在表3中。表3显示，对于实际地震信号，本文提出算法的去噪效果最好。

图7展示了σ=30时各种算法的二维地震信号去噪效果，可见本文算法有效去除了人为增加的高斯噪声。对比原地震信号可以发现，本文方法将地震信号中的高频干扰去除得较彻底，不存在ATV方法的块效应问题，且同相轴具有较好的横向连续性。

图8展示了图7中各个子图中道号为70的单道地震信号。从图8b可以观察到，在1850～1900ms区间内，地震信号被一个大幅度噪声污染，显然，基于逐点迭代的ATV和TGV方法都无法有效去除

该位置的噪声(图8c和图8e)。而基于交叠组稀疏收敛方法的ATV-OGS和TGV-OGS方法则非常好地去除了该位置的噪声(图8d和图8f)。对比图8d与图8f可见，本文提出方法恢复的地震信号更接近原地震信号。从图8还可知，交叠组稀疏正则项对于平滑区域的高噪声污染有较为显著的效果。

表3 实际地震信号算法性能测试表

图7 四种算法去噪效果对比(a)二维地震信号；(b)带噪数据；(c)ATV去噪结果；(d)ATV-OGS去噪结果；(e)TGV去噪结果；(f)TGV-OGV去噪结果

图9给出图8各个子图的频谱，可见加噪地震信号存在大量高频干扰，而经过各种去噪方法去噪后的地震信号频谱则一定程度压制了高频干扰。从图9f可见，在大于100Hz范围，本文方法压制噪声的效果明显好于ATV、ATV-OGS、TGV方法。

图10显示图7b中加入的噪声(σ=30)和各种方法去除的噪声剖面的对比。观察图10a的色标可发现，实际加入噪声的幅度范围是[-120,120]；而从图10b可知，ATV方法估计的噪声范围是[-60,60]，显然低估了噪声幅度；从图10c可见，ATV-OGS方法估计的噪声范围是[-70,70]，结果更准确些；观察图10d可知，TGV方法估计的噪声范围是[-50,50]；而本文提出方法(图10e)估计的噪声范围是[-100,100]，最接近实际加入噪声的幅度。因此，从噪声分布来看，本文方法估计的噪声范围是最接近实际添加噪声的。

图8 四种算法单道去噪效果对比

图9 图8数据对应的频谱

图10 四种算法去除噪声的对比

4 结论

本文从交叠组稀疏正则项出发，结合广义全变分定义，在ADMM框架下提出一种改进广义全变分去噪算法，并将其应用于地震信号去噪。该算法充分利用图像一阶、二阶梯度的邻域相似性，提高平滑区域与边界区域的差异性，从而提高去噪的鲁棒性，获得相比于经典TGV更好的去噪性能。

将本文方法与ATV、ATV-OGS、TGV算法进行比较，从实验结果得出如下结论和认识：

(1)将差分算子视为卷积算子，并结合ADMM能将复杂的模型转换为一系列简单的数学问题，并估计出较为准确的地震信号。

(2)本文算法去噪能力高于其他各类全变分去噪算法，尤其对平滑区域的大幅度噪声污染有显著的去噪效果，因此，提出算法对噪声更加鲁棒。

(3)取适当的K值可以有效提升整体去噪性能，在实际应用中需要调节参数取值，若K值取得过小，则邻域信息利用不够充分；反之，若取值过大，则会引入过多不相似的图像块，导致去噪性能下降。

值得注意的是，本文方法提出的通用正则项同样适用于其他各类图像重构问题，如地震波阻抗反演、自然图像去噪、椒盐噪声去除、图像解卷积，核磁共振图像重构等反问题的求解中，针对此于问题有待于深入研究。另外，由于引入MM算法，导致涉及两重循环，运算效率较低，因此算法的效率优化无疑是未来研究热点。