形态分量分析框架下基于DCT与曲波字典组合的地震信号重建

周亚同,刘志峰,张志伟

(1.河北工业大学电子信息工程学院,天津300401;2.中国科学院地质与地球物理研究所,北京100029)

形态分量分析框架下基于DCT与曲波字典组合的地震信号重建

周亚同1,2,刘志峰1,张志伟1

(1.河北工业大学电子信息工程学院,天津300401;2.中国科学院地质与地球物理研究所,北京100029)

针对单一型数学变换或字典不能有效刻画地震信号的形态特征多样性这一问题，在形态分量分析(MCA)框架下，提出了一种基于离散余弦变换(DCT)与曲波字典组合的地震信号重建方法。该方法首先建立MCA框架下的地震信号重建模型，依托模型将信号分解成局部奇异形态分量以及平滑与线状形态分量。然后采用DCT字典表示局部奇异分量，采用曲波字典表示平滑与线状分量。再以迭代求解方式逐一重建各分量，最后将重建后的分量合并。人工合成地震信号和二维叠前及叠后实际地震信号重建实验结果表明，该方法能很好完成信号重建，重建精度不仅要高于非抽样小波变换(UDWT)与曲波字典组合、曲波与曲波字典组合、余弦调制滤波器组与曲波字典组合，而且更高于DCT，UDWT，或曲波等单一型字典。

形态分量分析;地震信号重建;离散余弦变换;曲波变换;字典

地震信号重建在地震资料处理中具有重要意义[1-5]。现有地震信号重建方法大致分为以下3类：基于波场算子、基于预测滤波、基于数学变换的方法。在基于数学变换这类重建方法中，早年采用的是傅里叶变换[1]和Randon变换[6]等单一型经典变换。近年随着信号的多尺度几何分析理论的发展，涌现出一些新的变换，如曲波变换[7]、Seislet变换[8]和Dreamlet变换[9]等并被用于地震信号重建。另外随着稀疏信号表示及压缩感知理论的发展，又有学者尝试用超完备学习型字典来重建地震信号[10-11]。字典作为由若干个基函数所组成的超完备冗余框架，与单一型变换相比拥有更强的表达能力。但无论是叠前还是叠后地震剖面，通常都是由多种形态特征明显的元素组合而成的。以叠后剖面为例，既包含断点或尖灭点等局部奇异元素，也包含弱反射区或连续同相轴等平滑与线状元素。上述不同的元素对应着地震信号的不同内蕴特征，因此仅靠单一型变换或单一字典未必能有效刻画这种多样性的特征。

李海山等曾提出一种基于形态分量分析(morphological component analysis,MCA)的地震信号重建方法[12]：首先将地震信号分解成多个形态分量，然后从字典组合中挑选合适的字典分别重建各分量，最后将重建结果合并。这种分而治之的策略受数字图像修复启发[13]，用字典组合来表达地震信号中的不同内蕴特征，探索了地震信号重建的新途径。MCA框架下的信号分解与重建，面临的一个关键问题是字典组合的选择。字典组合选择不仅要匹配被表达信号的不同内蕴特征，而且要保证信号分解和重建可以数值化实现。李海山等在重建地震信号时选用的是非抽样小波变换(UDWT)与曲波字典组合，即用UDWT字典表示局部奇异分量，用曲波字典表示平滑与线状分量，取得了良好的重建效果。

本文在李海山等工作的基础上，尝试在MCA框架下用多种不同的字典组合重建地震信号，并对重建结果进行定量评价，以此为依据挑选出最优的字典组合。合成地震信号、叠前及叠后实际地震信号重建结果表明，本文选用的离散余弦变换(DCT)与曲波字典组合具有最优的重建精度。这可能是因为DCT字典依靠局域离散余弦变换，由多个固定大小的子块构成冗余字典。字典中稀疏的高频部分包含的是奇异性较强的局部细节特征信息，适合表示断点或尖灭点等局部奇异分量。而曲波字典具有很强的方向性，适合处理形态各异的同相轴纹理，即能对剖面的平滑与线状分量提供稳定、高效和近似最优表示。

1 形态分量分析框架下的地震信号重建模型

1.1 地震信号重建问题

二维地震信号重建问题可通过(1)式描述：

Y=MF

(1)

式中：F代表完整的二维原始地震采样信号;Y代表二维缺失地震采样信号;M为由元素0和1组成的采样矩阵，且元素满足：

(2)

式中：Ω代表地震信号中的缺损待重建区域。地震信号重建就是由缺失信号Y和采样矩阵M恢复出完整信号F的过程。

1.2 形态分量分析(MCA)框架基本原理

MCA是一种基于多种基函数的稀疏信号分解框架[14]，可以看作是经典基追踪算法(BP)[15]和匹配追踪算法(MP)[16]的结合。MCA利用信号组成分量的形态多样性，将信号分解成多个分量的线性组合，不同分量用不同字典稀疏表示。每个分量对应的字典仅能稀疏表示该分量，而对其它分量不能稀疏表示。现假设二维信号F可表示成K个形态分量Fk之和：

(3)

式中：αk为分解系数;Φk为过完备字典。可以通过求解(4)式约束规划得到系数αk：

(4)

式中：符号“‖‖1”代表L1范数。(4)式也可改写成以下形式：

(5)

式中：λ为给定阈值。上述求解系数{α1,…,αK}也可以转化为求解信号分量{F1,…,FK}问题：

(6)

1.3 MCA框架下的地震信号重建模型

根据MCA的基本原理，若在该框架下由缺失地震信号Y和采样矩阵M重建完整信号F，考虑到(1)式和(3)式，可以通过改写(5)式来实现重建：

(7)

MCA框架下的地震信号重建可用图1直观阐释：采用分而治之策略，首先将地震信号分解成多个形态分量Fk(k=1,2,…,K)，然后用不同的字典Φk分别表示各分量，最后将重建结果合并。图1中字典Φk对重建质量有重要影响。理论上可以根据一些预设准则自适应构造字典。但这样的构造方法通常比较复杂，在工程中往往根据经验选用一些现有变换基作为字典。如曲波字典、局部脊波字典、DCT字典、小波字典、小波包字典、Gabor字典等。

图1 MCA框架下的地震信号重建模型

MCA框架下的地震信号分解与重建，面临的一个关键问题是字典组合{Φ1,…,Φk,…,ΦK}的选择。本文尝试使用多种不同的字典组合进行地震信号重建，并对重建结果进行定量评价，以此为依据挑选出最优的字典组合。为计算方便，假设(7)式中形态分量数目K=2，所选用的4种字典组合分别为：DCT与曲波字典组合、UDWT与曲波字典组合、曲波与曲波字典组合、余弦调制滤波器组(CMFB)与曲波字典组合。其中，CMFB通过对线性相位低通原型滤波器进行余弦调制来实现。另外为便于对比，重建中还采用了单一型字典如DCT字典、UDWT字典和曲波字典等，此时形态分量数目K=1。

2 MCA框架下基于DCT与曲波字典组合的地震信号重建

根据MCA框架下的地震信号重建模型，在选定字典组合以后，首先需要构造出各字典，然后基于(7)式计算各字典对应的分解系数，再基于(3)式获得重建信号。目前大致存在两种构造过完备字典的方式：第一种是直接基于现有线性变换构造;第二种是基于学习方式，在给定样本集上学习出能够稀疏表示样本的字典[11]。由于第一种方式能借助快速变换和反变换实施数值计算，而且产生的字典往往高度结构化，因此本文采用第一种方式构造字典。

2.1 采用DCT字典表示局部奇异分量

DCT是数字信号处理领域常用的一种正交变换，有很强的能量集中特性，算法复杂度也比较低。设某信号块F(x,y)经离散化后共有M道，每道有N个采样点，对应的二维离散余弦变换定义如下：

(8)

其对应的反变换为：

(9)

DCT字典基于DCT变换构造。对于DCT变换后所获得的完备字典，采用分数频率法将其扩展成为过完备字典，具体做法是将得到的字典在其频率上做更精细的遍历和抽样，从而获得一个新的过完备字典。考虑到DCT字典自身特点，本文用其表示地震信号的局部奇异分量。

2.2 采用曲波字典表示平滑与线状分量

曲波变换是一种建立在脊波变换基础上的广泛使用的多尺度几何分析手段[17]。该变换将任意平方可积空间L2(R)中的函数映射为系数序列，实质上是通过把曲线无限分割，近似的把每段看作是直线段，然后对每个直线段进行脊波变换。曲波变换具有多尺度、多方向性、各向异性等特点，能有效表示和检测到信号中包含的多方向的线状变化特征[18-19]，因此十分适合用于表示地震信号中的线状同相轴。由于曲波变换在时间域内对信号的定义域进行正方形块分割，会造成边缘处出现不连续性，因此，Candès等[20]又提出了第二代曲波变换，把频率域划分为一组具有不同尺度的楔形区域，然后在每个区域上定义一个基函数。

曲波变换采用基函数和信号内积的形式实现信号稀疏表示，曲波变换可表示如下：

(10)

式中：符号“〈·〉”代表内积;cj,l,k为系数序列;φj,l,k为曲波函数，包含尺度j,方向l和位置k等参数。曲波变换紧框架通过标准正交基分解，可以表示一系列曲波函数的加权叠加，在L2范数意义下，对应的重建公式可以表示为：

(11)

(12)

曲波超完备字典基于曲波紧框架构造。首先对该紧框架中的尺度参数j,方向参数l以及位置参数k进行离散化，然后再平移滑动窗对最小尺度下不同方向的曲波原子进行截取，最终获得曲波超完备字典。由于曲波紧框架对信号的平滑与线状部分提供稳定、高效和近似最优表示，因此基于该框架得到的超完备字典，能够很好地匹配地震信号中的平滑与线状分量，能处理形态各异的同相轴纹理。鉴于曲波超完备字典的上述特性，本文用其表示二维地震信号中的平滑与线状分量。

2.3 MCA迭代重建的数值实现

在构造出各个字典后，需要基于(7)式计算各系数αk(k=1,2,…,K)，进而基于(3)式重建出完整的地震信号。直接求解(7)式十分困难，下面以形态分量数目K=2为例，分步阐述以迭代求解方式实现重建的全过程。

1) 根据地震信号缺失情况，确定(1)式中的采样矩阵M，同时根据地震信号的形态特征选择合适的字典组合{Φ1,Φ2}。

3) 初始化总残差R(0)=Y，其中Y为(1)式中的已知缺失地震信号。

4) 由硬阈值法得到起始阈值λ0，同时设置迭代停止阈值λstop。

5) 设置最大迭代次数Imax，在本文中最大迭代次数一般设置为50次。

8) 计算重建地震信号和原始地震信号间的各种误差及相似度，以此评价重建质量。

3 实验结果及分析

为了检验重建方法的效果，特设计人工合成地震信号和二维叠前、叠后实际地震信号重建3个实验。实验在CPU 2.53GHz，内存2GB，预装Windows 7旗舰版32位操作系统的个人计算机上进行，软件编程环境为MATLAB 7.10.0(R2010a)。实验中采用如下5个指标定量评价地震信号的重建精度：均方误差(MSE)、均方根误差(RMSE)、平均绝对误差(MAE)、峰值信噪比(PSNR)和结构相似度(SSIM)。

(13)

(14)

(15)

(16)

另外为度量地震信号重建前、后在整体结构形态方面的差异，定义二者间的结构相似度[21]为：

(17)

3.1 人工合成地震信号的重建过程展示

为清晰展示重建过程，选择如图2a所示某简单楔形合成地震信号，共有96道，单道160个采样点，道间距为10m，时间采样率为0.001s。图2b为随机缺失48道后的信号(缺失率达50%)。现采用本文算法对之进行重建，阈值类型设置为硬阈值，形态分量数目设为2，迭代起始阈值λ0设为1.3178，停止阈值λstop为10-6。图2c和图2d是在MCA框架下，分别采用UDWT与曲波字典组合、DCT与曲波字典组合分别进行重建的结果。从图2c和图2d可以看出，两种字典组合均能很好地完成重建。表1给出了采用两种字典组合重建时，合成地震信号的6个重建指标。从表1可以看出，DCT与曲波字典组合比UDWT与曲波字典组合的重建性能更好，重建时间也短一些(因为DCT的计算复杂度要比UDWT的计算复杂度低)。

图2 合成地震信号以及MCA框架下两种字典组合方式的重建结果

图3给出了在MCA框架下采用各种字典组合重建时，峰值信噪比(PSNR)和均方根误差(RMSE)随迭代次数的变化曲线。从图3中可以看出：①在MCA框架下采用字典组合方式，其重建效果比采用单一字典要好;②在4种字典组合方式中，DCT与曲波字典组合的重建效果最好，其次是UDWT与曲波字典组合，再次是曲波与曲波字典组合，CMFB与曲波字典组合的重建效果最差;③MCA框架下的3种单一字典重建方式的效果大致相当。

表1 MCA框架下采用两种字典组合时的合成地震信号重建质量评价

图3 各种字典组合重建的峰值信噪比(a)和均方根误差(b)随迭代次数的变化曲线

图4a给出了缺失地震信号的采样矩阵，即(1)式中的矩阵M。图4b给出了采用DCT与曲波字典组合进行迭代重建时，阈值λi随迭代次数的变化曲线。图4c为重建前、后地震信号的采样值统计直方图，其中左图是重建前原始合成地震信号(图2a)的采样值统计直方图，右图是由DCT与曲波字典组合重建后的地震信号(图2d)的统计直方图。对比发现二图非常接近，说明DCT与曲波字典组合确实能很好地完成重建。

图4 缺失地震信号的采样矩阵(50%道缺失)(a)、迭代重建阈值曲线(b)以及重建前、后地震信号采样值统计直方图(c)

3.2 二维叠前实际海洋地震信号重建结果比较

本实验选取某二维叠前实际海洋地震信号进行重建。图5a为原始地震信号，共有64道，每道含256个采样点，道间距为10m，时间采样率为0.001s。图5b为随机缺失19道后的信号。在MCA框架下，采用实验一中提及的各种字典组合方式对其进行重建，重建结果分别如图5c至图5i所示。对于DCT与曲波字典组合，所采用的重建参数为：阈值类型为硬阈值，形态分量数目为2，迭代起始阈值λ0为824.9917，停止阈值λstop为10-6。从图5 中可以看出，采用各种字典组合都能恢复缺失道，而且能保持与相邻道的振幅关系。但仔细对比各图发现，每种字典组合的重建结果均不相同。

图5 二维叠前实际海洋地震信号以及MCA框架下各种字典组合方式的重建结果

为了定量比较评价各种字典组合的重建效果，逐一计算6个重建指标并列于表2中。从表2可以看出：①在MCA框架下采用字典组合，其重建效果比采用单一字典要好;②在所有单一字典或字典组合中，DCT与曲波字典组合具有最优重建精度;③在相同的计算条件下，字典组合所需重建时间要比采用单一字典长。在实际计算中可通过减少迭代次数、选择合适的停止阈值等方式来减少计算时间。

3.3 二维叠后实际地震信号重建结果比较

叠后地震信号与叠前地震信号相比，通常同相轴的形态更加多样化。而且部分同相轴可能存在弯曲并伴随断点等，因此包含的元素类型更多，重建难度更大。选取如图6a所示某叠后实际地震剖面片段，该剖面片段共128道，单道含256个采样点，道间距为10m，时间采样率为0.001s。图6b显示的是缺道叠后剖面，在MCA框架下采用不同的字典组合对其进行重建。图6c和图6d 分别给出了UDWT与曲波字典组合、DCT与曲波字典组合的重建结果。从图6可以看出，两种字典组合均能完成重建，但重建结果仍有差异。

表2 MCA框架下采用不同字典组合时的叠前实际海洋地震信号重建质量评价

表3列出了采用不同字典组合时，叠后实际地震剖面的6个重建指标。无论是从MSE，RMSE或MAE，还是从PSNR或SSIM等指标来看，DCT与曲波字典组合均具有最优的重建精度。另外采用字典组合方式的重建效果要比采用单一字典好，但所需重建时间更长。上述结论与叠前地震信号重建一致。

图6 叠后实际地震信号以及MCA框架下两种字典组合方式的重建结果

表3 MCA框架下采用不同字典组合时的叠后实际地震剖面重建质量评价

字典组合类型MSERMSEMAEPSNRSSIM重建时间/sDCT与曲波字典组合3．42251．84980．801526．86500．9531104．7915UDWT与曲波字典组合5．47182．33921．119924．82610．9252124．5500曲波与曲波字典组合6．15642．48121．290624．31420．9159131．3001CMBF与曲波字典组合5．67312．49631．256124．58680．9209221．1349单一DCT字典6．67942．58451．331823．96010．909810．9905单一UDWT字典6．29752．50951．415824．21580．908426．6257单一曲波字典6．21202．48231．313024．29260．913524．3287

4 结束语

在MCA框架下，提出了一种基于DCT与曲波字典组合的地震信号重建方法。3个重建实验从不同角度验证了该方法的有效性，实验结果表明，DCT与曲波字典组合的重建精度不仅要高于UDWT与曲波字典组合、曲波与曲波字典组合、CMBF与曲波字典组合，而且要高于DCT，UDWT，或曲波等单一型字典。本文提出的重建方法在借助MCA框架发掘地震信号的内蕴特征、进而对信号进行更精细、灵活刻画方面提供了启示。下一步研究工作包括MCA框架下的字典自适应选择、地震信号分量分解质量的精细评价、迭代重建数值计算效率提升等。另外地震信号的反假频重建在实际资料处理中非常重要，下一步研究将考虑在本文的重建方法中引入反假频环节。

[1] Duijndam A J W,Schonewille M A,Hindriks C O H.Reconstruction of band-limited signals,irregularly sampled along one spatial direction[J].Geophysics,1999,64(2):524-538

[2] Zhang Y,Zhang G,Bleistein N.True amplitude wave equation migration arising from true amplitude one-way wave equations[J].Inverse Problems,2003,19(5):1113-1138

[3] Hindriks K,Duijndam A J W.Reconstruction of 3-D seismic signals irregularly sampled along two spatial coordinates[J].Geophysics,1999,64(1):253-263

[4] Zwartjes P,Gisolf A.Fourier reconstruction of marine-streamer data in four spatial coordinates[J].Geophysics,2006,71(6):171-186

[5] Schonewille M A,Romijn R,Duijndam A J W,et al.A general reconstruction scheme for dominant azimuth 3D seismic data[J].Geophysics,2003,68(6):2092-2105

[6] Kabir M M N,Verschuur D J.Restoration of missing offsets by parabolic Radon transform[J].Geophysical Prospecting,1995,43(3):347-368

[7] 张华,陈小宏.基于jitter采样和曲波变换的三维地震数据重建[J].地球物理学报,2013,56(5):1637-1649 Zhang H,Chen X H.Seismic data reconstruction based on jittered sampling and curvelet transform[J].Chinese Journal of Geophysics,2013,56(5):1637-1649

[8] Fomel S,Liu Y.Seislet transform and seislet frame[J].Geophysics,2010,75(3):25-38

[9] Wu R S,Geng Y,Ye L L.Preliminary study on dreamlet based compressive sensing data recovery[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013,3585-3590

[10] 唐刚.基于压缩感知和稀疏表示的地震数据重建与去噪[D].北京：清华大学,2010 Tang G.Seismic data reconstruction and denoising based on compressive sensing and sparse representation[D].Beijing:Tsinghua University,2010

[11] 周亚同,王丽莉,蒲青山.压缩感知框架下基于K-奇异值分解字典学习的地震数据重建[J].石油地球物理勘探,2014,49(4):652-660 Zhou Y T,Wang L L,Pu Q S.Seismic data reconstruction based on K-SVD dictionary leaning under compressive sensing framework[J].Oil Geophysical Prospecting,2014,49(4):652-660

[12] 李海山,吴国忱,印兴耀.形态分量分析在地震数据重建中的应用[J].石油地球物理勘探,2012,47(2):236-243 Li H S,Wu G C,Yin X Y.Morphological component analysis in seismic data reconstruction[J].Oil Geophysical Prospecting,2012,47(2):236-243

[13] 张涛,洪文学.基于自适应字典选择的MCA图像修复方法[J].光学技术,2010,36(5):672-676 Zhang T,Hong W X.Image inpainting based on MCA featured adaptive dictionary selection[J]. Op-tical Technique,2010,36(5):672-676

[14] Starck J L,Moudden Y,Bobin J,et al.Morphological component analysis[J].Proceedings of SPIE，2005(5914):1-15

[15] Chen S S.Basis pursuit[D].California:Stanford University,1995

[16] Mallat S,Zhang Z.Matching pursuit with time frequency dictionaries[J].IEEE Transactions on Signal Processing,1993,41(12):3397-3415

[17] Starck J L,Candes E J,Donoho D L.The curvelet transform for image denoising[J].IEEE Transactions on Image Processing,2002,11(6):670-684

[18] Candès E J,Donoho D L.Curvelets a surprisingly effective nonadaptive representation for objects with edges[M].Vanderbilt:Vanderbilt University Press,1999:71-89

[19] Candès E J,Demanet L,Donoho D L,et al.Fast discrete curvelet transforms[J].Multiscale Modeling and Simulation，2006，5(3):861-899

[20] Candès E J,Donoho D L.New tight frames of curvelets and optimal representations of objects with smooth singularities[R].California:Department of Statistics,Stanford University,2002

[21] Wang Z,Bovik A C,Sheikh H R,et a1.Image quality assessment:from error visibility to structural similarity[J].IEEE Transactions on Image Processing,2004,13(4):600-612

(编辑：陈 杰)

Seismic signal reconstruction under the morphological component analysis framework combined with discrete cosine transform (DCT) and curvelet dictionary

Zhou Yatong1,2,Liu Zhifeng1,Zhang Zhiwei1

(1.SchoolofElectronicandInformationEngineering,HebeiUniversityofTechnology,Tianjin300401,China;2.InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China)

Aiming at the problem that the mathematic transforms and dictionaries can not effectively depict the morphological features diversity of seismic signals,we propose a seismic signal reconstruction method under the morphological component analysis (MCA) framework combined with discrete cosine transform (DCT) and curvelet dictionary.Firstly the seismic signal reconstruction model is built under the MCA framework.Then the signal is decomposed into local singular and smooth linear component based on the model.Following that,the local singular component is represented by DCT dictionary,and the smooth linear component is represented by curvelet dictionary.We combine two kinds of components together after iterative reconstruction.The experiments on synthetic and real seismic signals illustrates that the proposed method can be used to reconstruct signals very well.The reconstruction precision of the method is not only higher than some dictionary combinations such as UDWT & curvelet dictionary combination,curvelet & curvelet dictionary combination,CMBF & curvelet dictionary combination,but also higher than some single dictionaries such as DCT,UDWT or curvelet etc.

morphological component analysis,seismic signal reconstruction,discrete cosine transform,curvelet transform,dictionary

2014-12-02;改回日期：2015-05-30。

周亚同(1973—)，男，博士，教授，主要研究方向为地震信号处理、模式识别与机器学习等。

国家自然科学基金项目(60972106，51475136)、河北省自然科学基金项目(F2013202254)和中国博士后科学基金项目(2014M56053)共同资助。

P631

A

1000-1441(2015)05-0560-09

10.3969/j.issn.1000-1441.2015.05.009