0 引言

海底地震仪(Ocean Bottom Seismometer, OBS) (郝天珧和游庆瑜，2011；游庆瑜等，2003)采集数据的预处理是利用OBS研究海域地球内部构造的重要基础，尤其是在黄渤海等浅海地区受底质、观测环境以及仪器耦合等影响，数据噪声较大，如何有效去除噪声并提取有效地震信号是OBS数据处理的关键环节之一.OBS预处理一般需经以下几步：数据解编处理、数据裁截处理、频谱分析和环境噪声分析，然后选择合适的带通滤波器滤波，对各OBS的SEGY格式数据做速度折合、自动增益、反褶积、时间校正等常规处理消除噪声影响，最后进行震相拾取工作(刘丽华等，2012).近年来，随着压缩感知的飞速发展，各研究领域均从理论上分析了各自的稀疏恢复能力，也大力推动了稀疏表达方法的发展(Patel and Chellappa, 2013)，同时为OBS去噪提供了新的思路.在OBS数据处理中，采集到的信号包含了许多冗余信息，因此需要找到一个完备或者过完备字典(变换)对其进行稀疏表达，进而获得更多准确的震相信息，方便后续的反演工作.

压缩感知(Compressed Sensing)理论指出：当信号本身具有稀疏性或者可以被某字典稀疏表示时，可以通过设计一个与稀疏字典基不相关的观测矩阵来观测信号，从而降低观测维度，节省观测成本.通过适当的促稀疏算法求解反问题，就可以很高概率地重构出理想数据(Candès，2006b; Donoho, 2006).

变换(字典)的选择直接影响数据的稀疏表示，地震领域常用的变换方法有离散余弦变换(Discrete Cosine Transform，DCT)、傅里叶变换(FT)、小波变换(WT)(Daubechies，1992)、曲波变换(Curvelet Transform)(Candès and Donoho, 2000; Candès et al., 2006a)以及近来深度学习中广泛采用的学习型超完备字典等.信号经过离散余弦变换后，能量主要集中在低频部分，从而达到稀疏表达的目的，但离散余弦变换是全局变换，因而不能很好地刻画信号的局部特征.2006年压缩感知框架建立以来，在医学核磁共振成像(Lustig et al., 2008)以及雷达(Oka and Lampe, 2009)领域都得到广泛应用，该领域大都选择使用傅里叶变换作为稀疏基，因为傅里叶变换可以把时频信息有机地联系起来，并且有成熟的快速算法，但傅里叶变换是整个时间域内的积分，不能描述局部特征，因此在做OBS数据的稀疏表达时傅里叶变换并不理想.Gabor(1946)提出了短时傅里叶变换，通过加窗的方式将信号划分为许多小的时间间隔，然后在每个时间窗内做傅里叶变换以提取信号的局部信息，但由于窗口的大小和形状都是固定的，不能随时间频率的变化而变化，因而不能从本质上识别信号的局部信息.小波变换继承并发展了Gabor的局部化思想(Daubechies，1992)，窗口大小即尺度固定，使形状随频率的变化而变化，根据频率的不同来调整时间分辨率，弥补了短时傅里叶变换的缺点，同时小波变换能够有效识别点异常信号，因而在地球物理领域中被广泛应用(徐亚等，2006；姚胜利，2007；吴招才和刘天佑，2008；张华等，2011).OBS信号中有用信息(地震波前)表现为曲线状地震震相，傅里叶变换和小波变换都不能稀疏地对其进行表示.Candès于2000年基于脊波变换(小波变换中一种)提出了具有多尺度、多方向和线性识别能力的曲波变换(Candès and Donoho, 2000)，之后又提出了全新的第二代曲波变换(Candès and Donoho, 2004)，被认为是地震信号稀疏表达最优方法之一，并被广泛应用于噪声去除(Herrmann and Verschuur，2004；Herrmann et al., 2007; Herrmann and Hennenfent, 2008)、数据重建(Shahidi et al., 2013)和波场模拟(Sun et al., 2009；Ma and Plonka, 2010)等方面.Shan等(2009)对各种变换在地震去噪过程中的应用做了详细对比.

稀疏约束可划分为L₀范数约束和L₁范数约束两大类.Mallat和Zhang(1993)在IEEE会议上提出了基于L₀约束的促稀疏算法：匹配追踪算法(Matching Pursuit，MP)，又称贪婪算法.该方法通过迭代所有可能，在字典域找出与原始信号最为匹配(内积最大)的基函数.该方法绝对收敛，同时比较稳定，但是由于L₀是NP-hard问题(非确定性难问题)，需要尝试所有次数，导致该方法计算量巨大.同年，Pati等(1993)提出了正交匹配追踪算法，顾名思义，该方法中每次迭代得到近似解后，其余部分与该解正交，因此同一个基函数在不同迭代过程中就不会被重复选择，相应计算量增大了不少，但是解的精确性大大提高.随后又有出现正则化匹配追踪(Needell and Vershynin, 2009)、逐步匹配追踪(Donoho et al., 2012)等等.Chen等(2001)提出基追踪算法(Basis Pursuit，BP)，BP方法基于L₁范数约束信号稀疏度，将非凸优化问题转化为凸优化问题，利用线性规划求解，是一种全局优化原则，复杂度高于匹配追踪类算法.随后，Daubechies等(2004)提出迭代阈值法，该算法因其本身算法简单且计算复杂度低而被广泛研究，基于该算法衍生出的算法有CRSI (Curvelet Recovery by Sparsity-promoting Inversion)(Herrmann and Hennenfent, 2008)、Bregman快速迭代(Yin，2010；Ma，2011)以及联合迭代法(白兰淑等，2014).

本文结合压缩感知(CS)和离散曲波变换的基本理论，提出利用迭代阈值结合改进后的冷却阈值法(CRSI)进行OBS数据去噪分析.通过对模拟数据进行去噪处理，对比了小波变换和曲波变换去噪结果，最后通过实际OBS数据去噪应用，检验本文方法的有效性.

1 基于稀疏表达的OBS数据去噪方法 1.1 压缩感知

压缩感知思想跳出了香农采样定理的要求，它基于信号的稀疏性，直接获取信号中有意义的信息，忽略无用信息.其核心思想是要求信号或其变换是稀疏的，即L₀范数很小，通过一个与变换基互不相关的测量矩阵投影到一个低维空间，然后求解L₀范数或者L₁范数最优化问题重构出原始信号.该投影过程即为采集过程，投影之后的空间维度远远低于原始信号，从而实现了数据的压缩.图 1为压缩感知理论信号采集与恢复的示意图.

地震数据中同相轴通常表现为线性异常，可以通过Curvelet变换来稀疏表示，因此本文选择Curvelet变换基作为稀疏字典，利用OBS数据在Curvelet变换域中的稀疏性，对OBS数据进行稀疏表示，并对Curvelet系数进行稀疏约束进而反演求其稀疏解，对优化后的Curvelet系数进行反变换以获取更准确丰富的震相信息，便于后续的速度模型反演.

1.2 Curvelet变换

第一代曲波变换(Candèses and Donoho, 2000)的数字实现比较复杂，而且需要分块平滑和Ridgelet分析等一系列步骤，同时会造成很大程度的数据冗余.因此，Candès和Donoho(2004)提出了新的曲波变换，即第二代曲波变换.两代变换原理上完全不同，第一代曲波变换的基本思想是将图像分为足够小的区块，将曲线近似成每块中的直线来看待，利用局部Ridgelet变换分析其特性.而第二代曲波变换则无需利用Ridgelet变换，实现起来要更为方便.

本文研究的是第二代曲波变换(二维)，设在二维空间中x为空间位置坐标，ω为频域参量，r和θ为频域极坐标.

定义曲波变换为：

(1)

根据帕萨瓦尔定律，可得到频率域曲波变换形式：

(2)

其中，Mother Curvelet为ϕ_j(x)，定义尺度为2^-j、方位角为θ_l、位置为x_k^{(j, l)}=R^-1_θl(k₁·2^-j, k₂·2^-j/2)的曲波为：

(3)

其傅里叶变换为，则在尺度2^-j上所有的Curvelet都可以由ϕ_j旋转和平移得到.引入等间隔旋转角序列：θ_l=2π·2^-⌊j/2⌋·l，l=0, 1, …, 0≤θ_l＜2π，空间位移参量序列：k=(k₁, k₂)∈Z².其中R_θl表示以θ_l为弧度的旋转.

设“半径窗”W和“角度窗”V均为平滑、非负和实值，并且满足允许性条件：

(4)

(5)

对所有尺度j≥j₀，定义傅里叶域窗为：

(6)

式中⌊j/2⌋表示j/2的整数部分，U_j在极坐标下是楔形窗.图 2为连续曲波空间频率域分块图(Candès and Donoho, 2004).

将式(2)离散化得到：

(7)

其中S_θl为旋转矩阵，.离散曲波变换频率域区域分块见图 3 (Candès et al., 2006a).

最终得到曲波变换重构公式为：.曲波变换具有多尺度、抛物尺度及多方向性等特征，能够稀疏地表示OBS信号中的曲线异常(震相信息).

1.3 去噪方法

设原始数据为y，有效数据及震相信息为f，噪声为n，则有：

(8)

根据OBS数据在Curvelet域的稀疏性分析可知：

(9)

D^H表示曲波反变换，x为能够稀疏表示OBS数据地震震相信息的Curvelet系数，那么，利用L₁范数约束条件，OBS数据去噪问题可转化为求解以下反问题：

(10)

该问题求解很大程度上取决于拉格朗日乘子λ，因为它调节着两项之间的权重关系，可以不断调整λ，使其满足精度要求：

(11)

采用基于Landweber下降法的迭代阈值技术(Daubechies et al., 2004)可以求解上述非线性问题：首先强调L₁项，即选取较大的λ，然后逐渐减小λ，即相对增加L₂项权重，以不断逼近真实解.

每一次迭代更新减小阈值，xⁿ←S_λ(x^n-1+D^T(y-Dx^n-1))，都会将二次方项最小化，继而可以通过阈值法投影到L₁球上(Pope，2009).这样虽然收敛，但要求较多的迭代次数，需要很长时间才能得到最终解.

在此基础上发展的冷却阈值法(Herrmann and Hennenfent, 2008)，在迭代过程中逐渐减小λ的值，以此作为阈值迭代法外循环的热启动(唐刚，2010；Tang and Ma, 2010a)，从而加快了整个过程的收敛性.本文中采用的迭代冷却阈值法可描述为：

输入：不完整数据y、内外最大迭代次数N和M；

初始化：m=1, ‖Dy‖＞λ₁＞…＞λ_M, x⁰=0

迭代：1:while‖y-D^Hx‖₂＜ε, and, m＜M

    2:for n=1, 2… do

    3:xⁿ←S_λ(x^n-1+ D^H(y-Dx^n-1)),

      S_λ(x)=sgn(x)·max(0, |x|-λ)

    4:endfor

    5:m=m+1

    6:endwhile

    7:

输出：，且

冷却阈值法首先设置一个与变换域最大系数相近的初始阈值λ₁来加强稀疏性，在之后的迭代中拉格朗日乘子逐渐减小，本文采用依次减小初始阈值的十分之一进行计算，逐渐减小误差.

接下来本文将Curvelet冷却阈值迭代的方法应用到模型数据和实际OBS数据中，验证方法的可行性.

2 模拟数据去噪分析 2.1 模拟数据

首先，利用模拟数据(图 4)完成去噪实验，图 4a所示为该模拟数据的速度模型，网格大小为400×500，网格间距为15 m.图 4b所示为炮点位于第126个网格点位置时的正演地震记录(已去直达波)，震源子波为雷克子波，共500道，道间距15 m，每道3000个采样点，采样间隔为1 ms.图 4b中的线性异常表示地下各界面对应的纵波反射震相(OBS中还包含横波反射震相以及更为复杂的折射波震相)，OBS预处理的最终目的就是将尽量多且信息准确的地震震相提取出来，以方便后续速度建模和结构反演(刘丽华等，2012).

2.2 去噪效果对比

为了对比小波变换和曲波变换的抗噪性能，本文在模拟数据中添加高斯噪声(频率域为白噪特征，存在于整个频率域，不会影响信号在Curvelet域的稀疏性)，使地震数据信噪比很小(小于0.1).这里定义信噪比为：

(12)

式中，y表示模拟信号(无噪声)，表示去噪后的信号.为方便实验，本文截取时间维度上2048个采样点进行测试(图 5a)，图 6所示为局部放大的结果，图像上显示，与小波变换结果相比，曲波变换结果更好地保持了地震同相轴.在反问题计算过程中，内外循环各4次，总计迭代16次，同时记录每次迭代结果的信噪比(图 7)，对比Curvelet和Wavelet的收敛速度及去噪结果信噪比.

采用小波变换和本文方法对模拟数据进行了去噪分析，结果见图 5c，5d.图 5b所示为添加高斯噪声后的正演数据，Curvelet冷却阈值迭代去噪的结果为图 5c，信噪比为71.7039，大于经过相同迭代次数的Wavelet(65.9019，图 5d).单纯分析收敛速度的话，小波变换要比曲波变换收敛快，但最终滤波结果信噪比比曲波变换差，从图 5c，5d对比也可以明显看出能够提取线性异常的曲波变换抗噪能力要优于小波变换，图 6所示为局部(偏移距100~165(×15m)，时间1.05~1.45 s)放大信号对比.由于OBS整体的数据量较小，收敛速度不是制约预处理的主要因素，而高的信噪比才是OBS预处理的真正目的，因此基于曲波变换的冷却阈值滤波方法更有优势.图 8所示为模拟数据的曲波谱，展示了地震数据所有尺度和所有方向的Curvelet系数，突出了地震数据的方向性.

3 实际OBS数据去噪

为了进一步验证本文去噪方法的效果，本节选取了在我国渤海海域采集的OBS数据进行去噪分析. 图 9、图 10所示为渤海海域某OBS数据，炮线长120 km，横轴为炮点号，纵轴为经过6 km·s^-1速度折合后的时间.图 9a为经过原始解编、带通滤波、增益和反褶积等处理后的原始数据，图 9b为图 9a经过本文方法滤波后的结果，震相的连续性得到了增强，同时抑制了能量较弱的多次波，从而丰富了震相的信息.图 10所示为Pn震相(浅蓝)、PmP震相(绿)、Pg1震相(红)以及Pg2震相(深蓝)在常规处理后和曲波冷却阈值滤波后的拾取情况，比较图 10a和图 10b发现，经过本文方法处理后，Pn震相的连续性得到增强，同时另外三种震相的信息更为丰富，能够为后续速度结构反演提供更多有用信息.图 11所示为该数据的Curvelet谱，从图中能够看出，本文将该OBS数据分为5个尺度，且每个尺度水平方向上能量最强，这一点与OBS数据震相大多近似水平(折合后)的性质相吻合，经过多次迭代，谱中能量较弱的高频噪声得以压制，OBS数据达到去噪效果.

4 结论

本文基于CS和Curvelet变换建立了一种基于稀疏表达的OBS数据去噪方法，由于曲波变换具有抛物尺度(各向异性)，且对线性异常(震相信息)敏感，因而能更好地稀疏表示OBS信号.本文采用的冷却阈值法相较于原先的迭代阈值，多了一层外循环的热启动以增强稀疏性，从而加快了收敛的速度和稳定性.分别选取曲波和小波作为变换基，在相同的迭代次数下，对模拟数据进行去噪对比，发现小波变换收敛略快于曲波变换，但是曲波变换得到的结果信噪比更高.由于OBS整体的数据量较小，收敛速度不是制约预处理的主要因素，高的信噪比才是OBS预处理的真正目的，因此基于曲波变换的冷却阈值滤波方法更有优势.最终通过OBS实际数据去噪测试，证明了基于Curvelet的冷却阈值迭代去噪法能够增强OBS信号中地震震相的连续性并丰富震相信息.