0 引言

地震数据在采集过程中通常无法避免地受到随机噪声干扰，去除随机噪声对后期的地震数据处理和解释至关重要.其中，基于稀疏变换的去噪方法由于较易实现、速度较快、效果较好而得到广泛的应用，如采用傅里叶变换的去噪方法(Liu and Fomel, 2013).但是，傅里叶变换是一种全局变换，缺乏对局部细节描述的能力，信噪分离时经常会损失大量的有效信号，因此，发展了拥有局部描述能力的短时傅里叶变换去噪方法(Vassiliou and Garossino, 1998).短时傅里叶变换虽然可以通过时窗来达到局部描述的目的，但是时窗大小难以自适应，小波变换去噪很好地解决了这个问题，信噪分离比短时傅里叶更加彻底(Xing et al., 2007; Gaci, 2014；汪金菊等，2016).当小波基扩展到二维的情况下时，小波变换对二维信号进行稀疏表示不能有效地描述二维信号的复杂纹理，去噪时难以达到令人满意的效果.为了弥补这一缺陷，在小波变换的基础上又发展了超小波变换的去噪方法，如采用curvelet变换(Tang and Ma, 2011; Gomes et al., 2017)、contourlet变换(Shan et al., 2009; Li and Gao, 2013; Zhao et al., 2016)、shearlet变换(刘成明等, 2014; Kong and Peng, 2015)、seislet变换(刘洋等，2009；Liu et al., 2015)等方法进行去噪.

各类超小波基可以较好的描述二维复杂信号，但是信号的纹理是复杂多变的，采用单一的固定基仍然难以对信号进行自适应的稀疏表示，在追求高精度的情况下寻找一种能够自适应稀疏表示信号的方法是必要的.字典学习能够根据不同的地震信号训练出相应的稀疏表示字典，如：K-奇异值分解(K-Singular Value Decomposition，K-SVD)的字典学习方法(Aharon et al., 2006).基于KSVD的地震去噪方法可以自适应地对地震信号进行稀疏表示，达到更好的去噪效果(Tian et al., 2017; Hou et al., 2018)，但是KSVD具有耗时过长的缺点，在地震数据较大时付出的时间代价太大.为了解决KSVD耗时过长的问题，在确保高精度的条件下，又发展了多种字典学习方法(Cai et al., 2014; Tong et al., 2016; Chen, 2017).

字典学习方法可以自适应地对地震信号进行稀疏表示，但是在去噪成像过程中，通常会因为缺少信号结构的先验约束信息而造成瑕疵，结合字典学习和固定基变换优点的双稀疏字典可以在自适应稀疏表示地震信号的同时消除瑕疵，获得更高的信噪比(Chen et al., 2016)，如：采用KSVD和小波变换的双稀疏字典(Zhu et al., 2015).

本文提出了一种由DDTF和contourlet变换相结合的双稀疏字典，并采用FISTA恢复系数的去噪方法.该方法首先对含噪数据进行contourlet变换获得contourlet系数，然后在稀疏编码阶段采用FISTA更新字典系数，再于字典更新阶段采用DDTF得到DDTF字典，接着再次使用FISTA可以得到与DDTF字典相对应的字典系数，并通过DDTF字典和与之对应的字典系数可以得到新的contourlet系数，最后对新的contourlet系数进行阈值处理和反contourlet变换得到去噪后的地震数据.通过模拟数据和实际数据的实验证明：与contourlet变换、KSVD字典、DDTF字典去噪方法相比，本文方法不仅拥有最高信噪比，而且在运算速度上远快于KSVD字典，与DDTF字典相差不大.

1 基于DDTF的双稀疏字典FISTA去噪

DDTF是一种高效的字典学习方法，与KSVD相比，DDTF在紧标架的约束条件下，只做一次SVD分解就可达到更新字典的目的，在速度上优于KSVD，并且精度与KSVD相差不大.目前，DDTF已广泛应用于地震数据处理中(Liang et al., 2014; Yu et al., 2016; Nazari Siahsar et al., 2017)，基于DDTF的双稀疏字典去噪可以表示为

(1)

其中，y_i为含噪数据Y的第i列，Φ和Φ^－1表示contourlet正、反变换算子，D和D^T表示DDTF字典和对应的转置，为原始数据对应稀疏系数的第i列，t为稀疏约束参数，控制着的个数，D^TD=I为紧标架约束条件，I为单位矩阵.

可将(1)式改写为

(2)

(2)式中，λ为正则化参数，一般情况下，L0约束可以等价为L1约束(Donoho, 2006).对(2)式中的变量D和分别进行求解，求解的过程称为稀疏编码阶段，首先设定一个初始字典D，固定D对进行求解：

(3)

对于(3)式，本文采用FISTA(Daubechies et al., 2004; 曲中党等，2015)进行求解，并按照Beck和Teboulle(2009)给出迭代结果:

(4)

(4) 式中，上角标k表示迭代次数，Ω_i^k－1表示前两次系数的叠加组合，β^k－1是初始值为1的参数，α为控制步长的参数，a_i^k为通过FISTA更新的系数，soft为软阈值函数，表示为

(5)

(5) 中x为待阈值函数，T为软阈值，然后固定a_i^k对D进行求解，求解D的过程称为字典更新阶段：

(6)

对于(6)式，本文按照Zhan和Dong(2016)进行求解，在得到更新的字典D^k+1后，再次使用FISTA，并进行硬阈值和contourlet反变换可得到去噪后的数据：

(7)

(7) 式中，表示对进行SVD分解，表示y_i去噪后的结果，T_h为硬阈值，hard为硬阈值函数，表示为

(8)

2 Contourlet变换

本文中的Φ选取contourlet变换，contourlet变换是一种超小波变换，拉普拉斯金字塔(Laplacian pyramid, LP)和方向滤波器(directional filter bank, DFB)是实现contourlet变换的两个主要步骤，图 1显示了contourlet变换的实现流程：

图 1中，↓(2, 2)表示下采样矩阵f，信号Y通过下采样分解成低通部分和差值部分，低通部分可称为信号的低频信息Y_l，差值部分是低通部分与上一层信号的差值，可称为信号的高频信息Y_e.LP通过不断地对Y_l进行下采样完成地震信号的多尺度分解，捕获信号中的线性特征，LP中的第j层低频信息Y_j^j和高频信息Y_e^j可以表示为

(9)

(9)式中，m, n, o, p表示位置参数，接着采用DFB捕获Y_e^j的方向信息，把各尺度上的方向进行分解，并将同一方向上的线性特征组合成相应的稀疏系数s，DFB矩阵可以表示为

(10)

(10)式中，d表示方向参数，DFB可以看成为具有方向性的基函数族：

(11)

h_d^j为DFB的分解合成函数，DFB的作用可以看成是不同方向上的基函数与Y_e^j进行内积，引入φ_{d, j}表示DFB在各方向上的基函数，可由内积的形式表示s：

(12)

3 算法框架

输入：含噪数据Y、FISTA参数(λ、α)、初始字典D、硬阈值T_h.

① 对Y进行contourlet变换，得到contourlet系数：s=ΦY；

② 根据s各尺度、各方向下的维度各自选择大小合适的窗口(如4×4、8×8等)，对s进行样本选取，并将选取的样本重排成新样本矩阵中的列，记新样本矩阵为s_r；

③ 根据式子(3)进行稀疏编码，得到更新后的字典系数a_i^k；

④ 根据式子(6)进行字典更新，得到更新后的字典D^k+1；

⑤ 根据式子(7)，在得到更新的字典D^k+1后，再次使用FISTA，并进行硬阈值，可得到去噪后的新样本系数；

⑥ 接着对进行重排复原，得到去噪后的contourlet系数，最后进行contourlet反变换得到去噪后的数据.

输出：更新后的字典D^k+1、去噪后的数据.

4 模拟数据实验

本文采用信噪比检验地震数据的去噪效果：

(13)

式子(13)中表示原始数据，本文所有的数据预处理为：对含噪数据进行8×8的时间窗口取样，DDTF、KSVD和本文方法均取总样本数的10%作为字典的训练样本，为了减少偶然性对实验造成的干扰，对每个实验均重复10次，信噪比取平均值.实验环境均为：Windows 7 SP1 64位系统、4G内存、英特尔第二代酷睿i5处理器、Matlab2013a.

为了研究本文方法的去噪性能，采用的模拟数据如图 2a所示，共有90道，每道采样点为256个，道间距为10 m，时间采样间隔为0.001 s.模拟数据还包含了4个同相轴，其中1个能量较强，3个较弱，红框内为相互交叉的两个能量较弱的同相轴，并用红色箭头标注出能量较强同相轴的顶点和模拟数据中能量最强部分.对模拟数据加入随机噪声，如图 2b所示，信噪比为4 dB，可以从含噪数据中分辨出能量较强的同相轴，识别出红色箭头所指部分，但是能量较弱的同相轴已经模糊不清，红框内无法识别出相互交叉的同相轴.

对含噪数据分别采用contourlet变换、DDTF、KSVD和本文方法去噪，去噪结果分别如图 3a、b、c、d所示，各自的信噪比和耗时如表 1所示，图 4b、c、d别为contourlet变换、DDTF、KSVD和本文方法去除的噪声部分，图 5a、b分别为DDTF字典和KSVD字典，由于篇幅有限，本文方法仅给出第一层contourlet系数的字典(下同)，如图 5c所示.DDTF字典和KSVD字典都是直接对含噪数据进行训练的，字典容易受到噪声的干扰，而本文方法利用了contourlet变换多方向多尺度的特性，在不同的方向和尺度下设置不同的阈值参数，有效地对噪声进行了滤除，提高了字典对噪声的鲁棒性，从各自的字典图中也可以看出，DDTF字典和KSVD字典中的原子构造与模拟数据中的纹理更为贴近，而本文方法的字典显示的是contourlet域中稀疏系数的纹理构造，在contourlet域中模拟数据和随机噪声得到了有效的区分，本文方法的字典更能揭示模拟数据的本质.

从图 3、图 4和表 1可以看出，contourlet变换虽然运算速度最快，但是它在去噪的同时也引入了较为严重的条带状噪声，去噪效果较差，能量较弱同相轴模糊不清，只能识别出能量较强的同相轴，红框内无法识别出相互交叉的同相轴，可以识别出红色箭头所指的能量最强部分，但是对于顶点只能勉强识别出来，从去除的噪声部分来看，contourlet变换对同相轴削弱较多，对红色箭头所指部分破坏较大.DDTF和KSVD去噪结果相当，存在着少量噪声，不仅能识别出能量较强的同相轴，能量较弱的同相轴也可以勉强识别出来，不过红框内的去噪结果表现为“缺失”，很容易让人误以为左边的同相轴和右边的是同一个，对于红色箭头部分，DDTF和KSVD都可以较为清晰地识别出来，从去除的噪声部分来看，DDTF和KSVD都对红色箭头所指部分造成了能量削弱，虽然去噪效果差不多，但是DDTF对比KSVD有着速度上的优势.本文方法去除了大部分的噪声，虽然引入了微量的条带状噪声，但是综合来看去噪效果最好，剩余噪声最少，能量较强和较弱的同相轴都可以识别出来，还能分辨出同相轴个数，红框内可以勉强识别出相互交叉的同相轴，对于红色箭头部分也可以较为清晰地识别出来，从去除的噪声部分来看，本文方法对能量最强部分造成了削弱，但是对顶点的破坏程度要小于DDTF和KSVD，本文方法耗时虽然比contourlet变换和DDTF多，但是远少于KSVD的.

对信噪比为1~7 dB的含噪模拟数据使用contourlet变换、DDTF、KSVD和本文方法去噪，信噪比变化如图 6a所示，图 6b显示的是随着FISTA迭代次数的增加，使用本文方法处理4 dB含噪模拟数据的信噪比变化.

由图 6可以看出，相比于contourlet变化、DDTF和KSVD，当处理信噪比为1~7 dB的含噪模拟数据时，本文方法在去噪精度上占有优势，并且，随着FISTA迭代次数的增加，本文方法的去噪效果也随之提高，当迭代次数超过13次时，去噪结果基本不再变化.

5 实际数据实验 5.1 实际陆地数据

为了研究本文方法对实际数据的去噪效果，采用实际陆地数据如图 7a所示，共有106道，每道采样点为256个，道间距为25 m，时间采样间隔为0.001 s.与模拟数据对比，陆地数据的纹理结构更为复杂，红框内为陆地数据中较为复杂的间断点纹理，用红色箭头标注出陆地数据的顶点以及左右两个底部.对陆地数据加入随机噪声，如图 7b所示，含噪数据的信噪比为3 dB，随机噪声干扰了红框和红色箭头处的纹理构造识别，特别是左边底部，已经模糊不清.

对含噪陆地数据分别采用contourlet变换、DDTF、KSVD和本文方法去噪，去噪结果分别如图 8a、b、c、d所示，各自的信噪比和耗时如表 2所示，图 9a、b、c、d分别为contourlet变换、DDTF、KSVD和本文方法去除的噪声部分，图 10a、b、c分别为DDTF、KSVD和本文方法的字典，可看出DDTF、KSVD字典与实际陆地数据的纹理十分贴切，而本文方法的字典显示的是实际陆地数据在contourlet域中的构造.

从图 8、图 9和表 2可以看出，contourlet变换有速度上的优势，但是去噪精度较低，在去除随机噪声的同时还引入了相当强度的条带状噪声，从去除的噪声部分来看，很多纹理细节被当成噪声部分去除，对红框和红色箭头标注的纹理构造都有较为严重的削弱.DDTF和KSVD去噪效果差不多，精度均比contourlet变换高，去除了大量的随机噪声，从去除的噪声部分来看，对红框和红色箭头标注的纹理构造破坏程度较contourlet变换轻，同时耗时也比contourlet变换多，但是DDTF耗时远少于KSVD.本文方法拥有最好的去噪效果，去除了绝大部分的随机噪声，虽然也引入了微量的条带状噪声，但是综合来看去噪精度最高，损害的纹理细节部分最少，从去除的噪声部分来看，虽然对红框内的间断点纹理有微量的能量削弱，但是几乎没有削弱红色箭头标注的顶点和左右两个底部的能量，在耗时上，本文方法虽然比DDTF多，但是仍然远少于KSVD.

对信噪比为1~7 dB的含噪陆地数据使用contourlet变换、DDTF、KSVD和本文方法去噪，信噪比变化如图 11a所示，图 11b显示的是随着FISTA迭代次数的增加，使用本文方法处理3 dB含噪陆地数据的信噪比变化.

由图 11可以看出，本文方法对信噪比为1~7 dB的含噪陆地数据有效，当数据纹理较为复杂时，相比contourlet变换、DDTF和KSVD，本文方法拥有更好的去噪效果，并且，随着FISTA迭代次数的增加，本文方法的去噪效果也随之提高，当迭代次数超过15次时，去噪结果基本不再变化.

5.2 实际海洋数据

为了进一步研究本文方法的去噪效果，采用实际海洋数据如图 12a所示，共有120道，每道采样点为256个，道间距为20 m，时间采样间隔为0.001 s.与陆地数据相比，海洋数据的纹理结构比较简单，多为直线纹理构成，用红色箭头标注出海洋数据中能量较强的顶点和能量较弱的底部.对海洋数据加入随机噪声，如图 12b所示，含噪数据的信噪比为1 dB，随机噪声对红色箭头处的纹理识别造成了干扰，一些能量较弱的同相轴无法识别.

对含噪海洋数据分别采用contourlet变换、DDTF、KSVD和本文方法去噪，去噪结果分别如图 13a、b、c、d所示，各自的信噪比和耗时如表 3所示，图 14a、b、c、d分别为contourlet变换、DDTF、KSVD和本文方法去除的噪声部分，图 15a、b、c分别为DDTF、KSVD和本文方法的字典，可看出DDTF、KSVD字典与实际海洋数据的纹理十分相似，而本文方法的字典显示的是实际海洋数据在contourlet域中的情况.

从图 13、图 14和表 3可以看出，contourlet变换的去噪速度最快，但是去噪效果最差，虽然它几乎未损害纹理细节，但是去除的随机噪声最少，而且还引入了相当强度的条带状噪声.DDTF的去噪速度远快于KSVD，去噪效果两者差不多，都去除了相当强度的随机噪声，同时也对纹理细节造成了损害，削弱了顶点和底部的能量.本文方法的去噪速度慢于DDTF，却远快于KSVD，在去噪精度上仍然是本文方法最高，损害的纹理细节较少，稍微削弱了顶点处的能量，虽然引入微量的条带状噪声，但是却去除了绝大部分的随机噪声，综合来看本文方法拥有最好的去噪效果.

对信噪比为1~7 dB的含噪海洋数据使用contourlet变换、DDTF、KSVD和本文方法去噪，信噪比变化如图 16a所示，图 16b显示的是随着FISTA迭代次数的增加，使用本文方法处理1 dB含噪海洋数据的信噪比变化.

由图 16可以看出，本文方法对信噪比为1~7 dB的含噪海洋数据依然奏效，与contourlet变换、DDTF和KSVD相比，本文方法拥有最高的去噪精度，并且，随着FISTA迭代次数的增加，本文方法的去噪效果也随之提高，当迭代次数超过18次时，去噪结果基本不再变化.

6 结论

在模拟数据与实际数据的实验中，通过与contourlet变换、DDTF和KSVD进行对比，我们可以得到以下结论：在精度上，无论是模拟数据还是实际陆地、海洋数据，本文方法都拥有最高的信噪比，对噪声的鲁棒性最好，去噪时损害有效信号较少；在速度上，contourlet变换的耗时远少于字典学习方法，在字典学习方法中，DDTF和本文方法都要优于传统的KSVD，本文方法耗时虽然多于DDTF，但是付出的时间代价在可接受范围之内.虽然本文方法有着高精度并且兼顾速度的优点，但是在超小波变换后的多尺度下所需调整的参数过多，并且严格来说，本文方法在速度上还有提升的空间，如何简化参数和进一步地提升速度将会是下一步的研究方向.