0 引言

由于受地质环境与采集条件的影响，地震数据易出现缺失道或空间采样不足等问题，而后续地震数据的处理与解释又对数据的完整性与规则性有很高的要求，因此地震缺失道数据的恢复重建是非常关键的问题.传统地震数据重建方法(Ronen, 1987; Spitz, 1991)受奈奎斯特(Nyquist)采样定理的限制，采样频率至少是信号最高频率的两倍，否则极易出现假频现象而影响重建质量.压缩感知(Compressed Sensing，CS)理论(Candes et al., 2006; Donoho, 2006)依靠信号本身的稀疏性和观测方式的非相关性，只要信号在某个变换域上是稀疏或可压缩的，就可以用一个与稀疏表示基不相关的测量矩阵对信号进行降维观测，并在接收端通过求解最优化问题，以较高的概率根据其低维观测值，高精度重建原始信号，为欠采样率下地震数据的重建问题提供了新思路.目前基于压缩感知的地震数据重建方法逐渐引起国内外学者的关注(Herrman and Hennenfent, 2008; 刘成明等, 2016)，其重点研究的问题集中于：(1) 地震数据不相关的观测方法(王本锋等, 2015)；(2) 地震数据最优的稀疏表示(Ma and Plonka, 2010)；(3) 重建模型最优化的求解算法(白兰淑等, 2014)，上述一系列问题中地震数据最优的稀疏表示方法是提高重建效果的关键.传统地震数据稀疏表示的方法有Fourier变换(陈双全和李向阳, 2015)、Radon变换(陈潜等，2015)、Wavelet变换(王晓凯等，2016)等，Fourier变换是一种全局变换，没有局部化分析功能，基于Fourier变换的重建方法容易产生吉布斯现象；Radon变换能够较好地捕捉地震数据中的直线特征，但地震波前在空间上通常呈曲线状；Wavelet是一种时频域的局域性变换，因而能有效表示各向同性的点奇异特征，虽然在地震数据处理方面表现出一定能力，但在边界和波前奇异线条附近存在不光滑和模糊现象.近年来发展的多尺度几何分析技术由于具有良好的多尺度性、多方向性和各向异性，更加适合稀疏表示地震波前特征，在地震数据处理中得到了广泛应用，如Ridgelets变换(包乾宗等, 2007)、Curvelet变换(张岩等, 2015；吴蔚等, 2016; Liu et al., 2016)等方向性变换也逐渐应用于地震数据重建等处理，并取得了较好的效果.虽然Ridgelets变换对直线奇异函数具有良好的逼近性能，但对于地震数据包含的曲线奇异函数，不具有最优的逼近能力.Curvelet变换(Candès and Donoho, 2004)由于具有良好的方向性、局部性与各向异性, 可以捕获到各个方向上地震数据的同相轴信息，进而很好地稀疏表示，在地震数据处理领域逐渐得到重视，但是Curvelet变换对于简单的纹理模型，如地震波前，不具有最优的稀疏表能力.

波原子变换(Wave Atoms Transform) (Demanet and Ying, 2007)是一种特殊的二维波包的变形，与Curvelet相比，其基的支撑区间是各向同性的，而其每个波包的振动周期和支撑尺寸满足抛物尺度关系，即波长约等于支撑尺寸的平方，在这个意义下可以简单地将波原子理解为方向小波和Gabor原子的插值.相对Wavelet、Gabor和Curvelet而言，波原子对于简单纹理模型具有最优的稀疏表示，对于给定的精度，只需要O(N)个波原子系数就可以表示，但需要O(N^3/2)个Curvelet系数，O(N²)个Wavelet系数和Gabor系数才能达到同样的精度.陈书贞等(2010)将波原子变换引入图像的压缩感知重构处理中，由于波原子可以有效稀疏表示图像的纹理区域，因此得到了较好的重构效果.杨宁等(2011)利用波原子具有准确稀疏表示地震波前曲线的特点，通过系数相关性的方法对叠前地震数据进行去噪处理，具有较高的信噪比.

本文提出基于波原子域的地震数据压缩感知重建算法，在波原子域建立基于压缩感知框架的地震数据重建正则化模型，通过Landweber迭代算法反演求解L1范数最小优化问题，结合指数阈值收缩模型逐步促进编码系数的稀疏程度，保留地震数据主要特征，采用循环平移技术抑制重建数据中的噪声，实验结果表明本文算法重建数据反射同相轴清晰、连续性好，有利于后续地震资料的处理和解释.

1 波原子变换 1.1 波原子变换原理

波原子(Wave atoms)是一种特殊二维波包的变形，实现了多尺度、方向性、和局域性较好的权衡.每个小波包的振动周期和支撑尺度满足抛物线函数，即波长约等于支撑尺寸的平方，相对于Wavelet，Gabor和Curvelet而言，波原子对于振荡函数(简单的纹理模型)具有最优的稀疏表示，因此具有较好的纹理保持性能.

1.1.1 波原子定义

在二维空间考虑波原子形式的定义，记x=(x₁, x₂)，Fourier变换为

(1)

设波原子为φ_μ，下标μ=(j, m, n)=(j, m₁, m₂, n₁, n₂)，其中j、m₁、m₂、n₁、n₂均为整数，定位相空间中的点(x_μ, ω_μ)，其中，分别为大于零的常数; x_μ、ω_μ分别表征波原子在空间域和频域中的中心.此外，二维波原子在相空间(x_μ, ω_μ)附近还必须满足局部化条件.假定波包{φ_μ}的框架元为波原子，若对任意M>0，φ_μ在空间域和频域中的局部化条件分别为

(2)

(3)

以上两个条件是对波原子在时频局部化约束下的定性描述.

1.1.2 波原子构造

在一维空间中，假设g为实值无穷可微的冲击函数，支撑区间为[-7π/6, 5π/6].当|ω|≤π/3时，满足：

(4)

利用g定义：

(5)

其中，

选择适当的m，使，则{φ_m⁰(t-n)}，n∈Z，构成L²(R)空间的一组正交基.于是频域轴的多尺度覆盖可以通过二进缩放实现.引入尺度指标为j，记基函数为

(6)

虽然{φ_{m, n}^j(x)}也称作波包(Villemoes，2002)，但它与标准小波包有着本质的差异，其在空域与频域都具有一致有界局部化性质，这一点对于波原子的构造及其性质起着至关重要的作用.对于一维信号f(x)，其波原子变换可视为

(7)

进一步，由Plancherel公式，有：

(8)

在二维空间中，波原子的构造可以通过一维波原子实现.设一维波原子为φ_m^j(x-2^-jn)，μ=(j, m, n)，m=(m₁, m₂)，n=(n₁, n₂)，令：

(9)

(10)

其中H{·}表示Hilbert变换.则ψ_μ⁺(x₁, x₂)和ψ_μ^-(x₁, x₂)都是规范正交基(事实上，它们是一对张量形式的波原子正交基)，于是其组合和，构成了冗余度为2的波原子紧框架，它们在频域中有对称的两部分，即与一维波原子有相似的结构与性质，由此可见波原子变换的基函数具备完备性与反变换的唯一性.

1.2 波原子域稀疏表示地震数据

地震数据中的同相轴构成丰富的曲线状纹理，波原子变换对地震数据具有较好的稀疏逼近特性.如图 1a所示的地震数据包含256个采样点、128个记录道，分别采用小波变换、波原子变换、Curvelet变换3种不同的稀疏表示方法得到的系数数目如表 1所示，Curvelet得到的系数最多，小波得到的系数最少.分别取其中2000个最大的系数重构数据，效果分别如图 1b、c、d所示，可见波原子重构效果最好，Curvelet的曲线特征保持能力优于小波，而小波变换重构的整体SNR高于Curvelet. 图 2给出了三种不同稀疏表示方法重构的SNR随系数数目变化的曲线，可见在系数数目相同的条件下，波原子重构得到的SNR均高于小波与曲波变换.从而说明波原子变换得到的变换系数冗余程度相对低，并且能够更准确地稀疏表示地震数据.

2 波原子域地震数据压缩感知重建 2.1 波原子域压缩感知重建模型

设完整的二维地震数据为f(N_t×N_r)，N_r是检波器坐标样本数，N_t是时间轴坐标样本数，设V是向量化算子，f向量化后得x=V(f)，x为(N_tN_r)向量.欠采样地震资料设为q(N_t×M_r)，M_r < N_r，向量化后得y=V(q)，y为(N_tM_r)维向量，则欠采样地震数据可以表示为

(11)

R是一个二元采样矩阵，从x采样得到y，相当于从完整的地震资料f中N_r道抽取M_r道数据.因为M_r < N_r，显然式(11) 是欠定的方程组，无法根据y求解x，设地震资料x在Wave atoms变换域W下系数为α，有x=W^-1α，W^-1是W的反变换，则:

(12)

在具备地震数据在Wave atoms域变换系数是稀疏的先验知识，和感知矩阵A^CS满足有限等距性质(Restricted Isometric Property, RIP)(Candès, 2008)的前提下，为通过y恢复x提供了理论保证.求解欠定方程组(11) 的问题转化为L1范数优化问题，公式为

(13)

求解式(13) 是极不稳定而且是NP难的非凸优化问题.在R和W不相关的条件下，可以转为求解一个等价的、更加简单的L1范数优化问题，公式为

(14)

使问题转化为求解一个凸优化问题，进而可解.由于α是稀疏的，可以通过不断促进α的稀疏性通过式(15) 进行反演求解.公式为

(15)

其中λ是平衡因子，是α的估计值.在α的稀疏性促进过程，采用Landweber下降法的冷却阈值迭代技术，第k次迭代过程公式为

(16)

其中T为阈值算子，公式为

(17)

冷却阈值法每次迭代缩减阈值Threshold，可以找到离超平面y=Rx距离最近的向量，对应式(15) 中项，然后投影到L1球上，得到是α的估计值，利用反变换后得重建后完整地震数据向量，反向量化得重建地震数据.

2.2 循环平移技术

本文提出的地震数据重建过程中，利用波原子域冷却阈值迭代法反演求解L1最优化问题，由于波原子变换缺乏平移不变性，易在地震数据缺失道邻域产生伪吉布斯现象，导致重建数据失真，该失真与地震数据缺失道的位置密切相关.在图像处理领域Coifman和Donoho(1995)提出的循环平移技术(Cycle Spinning)，通过对含噪声图像进行循环平移、阈值去噪、逆向循环平移，每次平移后的图像进行阈值去噪会使伪吉布斯现象出现在不同的地方，针对图像行和列方向上的每组平移量都会得到一个不同的去噪结果，对所有去噪结果进行线性平均，抑制噪声.循环平移技术对地震数据可分别在时间轴和检波器轴上进行平移，设S_{i, j}表示平移量为i、j(i=0, 1, 2…K₁，j=0, 1, 2…K₂，K₁和K₂分别为时间轴和检波器轴最大平移量)的循环平移算子，为逆循环平移算子，地震数据重建迭代过程式(16) 在平移量为i和j条件下的可以表示为

(18)

其中，为循环平移原始地震数据向量，R_{i, j}=S_{i, j}(R)为循环平移观测矩阵，y_{i, j}=V(S_{i, j}(V^-1(y)))为循环平移的观测地震数据向量.

最大平移量K₁和K₂下(K₁=0, 1, 2…N_t，K₂=0, 1, 2…N_r)，对应的重建地震数据向量可以表示为：

(19)

2.3 指数阈值收缩模型

地震数据重建过程中，阈值收缩迭代参数的选择方式对重建效率有重要影响.由于欠采样地震数据曲线状的波前信息被随机缺失的地震道所打断，在波原子域表现为噪声特征，且迭代初期该噪声的强度相对较大，阈值收缩的速度需要相对加快，加速噪声的去除；在迭代接近结束时，噪声相对较小，阈值收缩速度应该适当放缓，加强地震数据的细节信息的保留.线性变化的阈值收缩模型无法满足上述需求，Gao等(2010)在凸集投影(Projection Onto Convex Sets，POCS)框架下，利用指数型阈值集合模型在Fourier域进行地震数据重建，相对于线性模型阈值可以有效减少迭代次数，提高重建效率，因此本文提出波原子域的指数阈值收缩模型，指数阈值函数表示为

(20)

其中k为迭代次数，Threshold⁽¹⁾为初始阈值，δ为阈值收缩因子.

3 波原子域地震数据压缩感知重建算法

基于波原子域稀疏表示的地震数据压缩感知重建算法的总体思想是：通过Landweber迭代算法反演求解提出的正则框架下L1范数最优化问题，每次迭代结合循环平移技术抑制重建数据中的噪声，采用指数阈值收缩模型逐步促进编码系数的稀疏程度，逼近原始地震数据.具体算法描述如下：

【算法1】：基于波原子域的地震数据压缩感知重建算法

算法步骤：

输入：给定迭代停止参数τ，最大迭代次数L，阈值收缩因子δ，地震数据观测采样矩阵R，地震欠采样观测数据q.

输出：重建地震数据

(1) 初始化：迭代计数器k=1，初始阈值Threshold⁽¹⁾，阈值收缩因子δ，最大循环平移量K₁和K₂，Wave Atoms变换算子为W，向量化算子V，y=V(q)，重建初值为x⁽¹⁾=R^T(y).

(2) 重建处理过程：

4 实验结果与分析

实验的硬件平台采用双核CPU主频3.3 G的Intel I5微机，内存容量为4 G.系统软件为32位Windows7操作系统，仿真实验软件使用Matlab R2013b.实验用例分别采用合成地震数据、marmousi模型数据，以及实际地震数据.地震数据重建效果的衡量指标采用信噪比(SNR)与峰值信噪比(PSNR)分别为

(21)

(22)

其中MSE为原始完整地震数据与重建后地震数据均方误差.

4.1 合成地震数据实验

图 3a是截取原始合成单炮地震数据部分记录，包括128个记录道，128个采样点，图 3b为随机去除50%记录道后得到的欠采样数据，图 3c为本文算法迭代200次的效果，重建后PSNR约提高11.9 dB左右，SNR约提高6.6 dB左右，且地震数据的同向轴方向性与连续性保持较好，图 3d给出了原始地震数据与重建地震数据的差值图，可见本文算法重建得到的误差较小，图 4a是SNR随迭代次数的变化曲线，随着迭代次数的增加，重建效果不断增强，最后算法取得收敛，SNR基本稳定在最大值，从而说明本文算法具有较好的收敛性与稳定性.图 4b给出了阈值随迭代次数增加的变化曲线，迭代初期阈值缩减相对剧烈，以加速去除缺失地震道而产生的噪声，迭代后期阈值收缩相对缓慢，有利于地震数据细节的保留.

4.2 Marmousi地震模型数据实验

在Marmousi模型数据中随机抽取单炮数据，截取原始地震数据200道，每道256个采样点，如图 5a所示，去除抽取40%记录道后得到的欠采样数据如图 5b所示，图 5c给出了本文算法重建效果.由于双树复小波(Dual-tree Complex Wavelet Transform, DTCWT)和Curvelet变换具有较好的方向性，二者是地震数据处理领域应用比较广泛的稀疏表示方法，接下来给出与DTCWT、Curvelet两种变换的重建效果对比，图 5d是DTCWT稀疏表示后阈值收缩迭代重建效果，可见重建后数据质量得到一定的改善，但在同相轴密集区域误差较大.图 5e是Curvelet稀疏表示后阈值收缩迭代重建效果，重建地震数据中同向轴比较连续，但在连续缺失多道处仍存在较大的误差.

表 2分别给出了DTCWT重建、Curvelet重建、未采用循环平移技术的波原子域线性阈值收缩与指数阈值收缩，及本文算法在不同的地震道采样率下的重建效果对比，可以得出，随着采样率的增加，各个算法的重建效果不断增强，本文算法在各个采样率下均能取得相对好的重建效果.

表 2所述各种重建算法在一次迭代过程中，DTWCT重建算法的运行时间为0.06 s，Curvelet重建算法的运行时间为0.13 s，未采用循环平移技术的波原子域线性阈值收缩与指数阈值收缩的运行时间为0.20 s，本文算法为0.81 s.不同重建算法计算量的主要区别在于进行变换与反变换运算的次数及复杂度.本文算法由于采用循环平移技术，在一次迭代中需要进行多次变换与反变换，而前4种算法在一次迭代过程仅需要一次变换与反变换，因此运行时间少于本文算法.由于DTWCT变换、Curvelet变换，波原子的计算复杂度依次增加，故运行的时间也依次增加.考虑地震道的缺失数据仅存在于检波器轴，而时间轴数据完整，本文仅在检波器轴进行循环平移，每次迭代采用4种方式平移，每种方式均需要进行一次变换与反变换，因此本文算法运行时间是未采用循环平移技术的4倍左右.

4.3 实际地震数据实验

图 6a为某区海洋实际叠后地震数据，选取500道，每道400个采样点，采样间隔为1 ms，该数据中既包含强能量、连续性好的同相轴，又包含能量弱的、连续性差的信号，及部分噪声，重建处理难度较大.图 6b为随机去除50%的地震道后的数据，图 6c为本文算法重建结果，从图 6c中可见各缺失道均得到了很好的恢复，地震同相轴较为光滑、连续，其中强振幅同相轴恢复效果尤为明显.

5 结论

为提高缺失地震数据的重建效果，本文在压缩感知框架下研究了基于波原子稀疏表示的地震数据重建.通过理论分析与模型实验可以得到以下结论：

(1) 利用波原子域有效稀疏表示地震数据的特点，建立的压缩感知重建求解模型是有效的.

(2) 结合循环平移技术，可以克服波原子缺乏平移不变性而引起重建失真的问题，抑制重建数据中的噪声.

(3) 地震数据重建迭代过程中，指数阈值收缩模型更适合促进编码系数的稀疏程度，去除波原子域由缺失道引起的噪声，反演逼近原始地震数据.由于本文算法在每次循环平移过程均需要进行波原子变换与反变换操作，从而增加了算法的运行时间，因此提高算法的运行效率是我们下一步研究的重点.

致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！