唐国维， 程鑫华， 张岩

(东北石油大学 计算机与信息技术学院， 大庆 163318)

0 引言

现在的高分辨率地震勘探技术对地震数据的规则性以及完整性提出了更高的要求。然而由于客观地质条件、勘探成本以及其他环境因素，在地震勘探时，无效的地震道采样、坏道、以及各种不规则的干扰数据通常都会对地震数据进行一定程度的干扰，使得野外勘探得到的数据往往是不规则和不完整的。若不对缺失的地震数据进行处理，会影响多次波压制、偏移等处理结果，并最终影响地震资料的解释和油气藏的评价。在实际勘探过程中，处理不规则采样数据的一些相对简单的方法，如拷贝或线性插值出邻近道、忽略空缺道等做法均有不足之处。因此，寻求一个快速而有效的地震数据重建方法，对于地震数据的分析以及后续工作具有重大意义。

地震数据重建，即利用一定的方法和手段对已有的道缺失地震数据进行插值重建处理，生成完整的或者采样率更高的地震数据。地震数据的重建方法大致一共有3类：第一类是基于滤波器的策略：将欠采样地震数据与某种插值滤波器做卷积，常用的一种是预测误差滤波器，但这类方法通常将随机采样数据当作规则数据来处理，并使用高斯窗进行插值操作，这种操作有着非常大的不确定性，会造成很大的误差；第二类是基于波动方程的方法：通过DMO或AMO正、反演算子迭代进行插值处理，此种方法能够处理随机采样的地震数据，但是计算量巨大，并且需要有地下结构的先验信息，对于使用较粗网格的欠采样地震数据的重建效果并不理想；第三类是基于变换函数的方法：利用特定的变换函数对欠采样地震数据进行变换，在变换域内进行插值等处理之后再进行逆变换，最终得到重建后的地震数据。

在过去几十年里，奈奎斯特(Nyquist)采样定理在采样、存储、传输和处理等信号处理的传统领域中起到了非常关键的作用。最近发展出的压缩感知(Compressed Sensing,CS)理论表明[1,2]，若待处理的数据是稀疏的，通过一定的采样方法，就可以突破奈奎斯特定理，利用更少的数据来重建出满足一定精度要求的重建数据。基于以上特点，将压缩感知理论应用于地震数据重建[3,4]:国内已有白兰淑提出了一种基于曲波变换的联合迭代重建算法[5]，周亚同提出了一种基于K-奇异值分解的字典训练的地震数据重建方法[6]，唐刚则使用泊松碟采样的方式对地震数据进行重建等。因此压缩感知的稀疏表示部分通常是先将欠采样的地震数据进行某种变换来满足压缩感知理论中对于数据稀疏性的要求。通常使用的变换方法有离散余弦变换(DCT)、傅立叶变换、离散小波变换等。离散小波变换作为信号和数据处理工具取得了巨大的成果，离散小波变换由于其时频局部化特性、多分辨率特征性、边缘检测特性等使得其在地震数据重建方面有着重要应用。但是传统的离散小波具有的自身局限性，主要体现在：缺乏平移不变性；缺乏方向选择性；震荡性；频谱混叠性。普通的复小波变换可以解决上述局限性的问题。但是由于超过一层分解的复小波变换的输入都是复数形式，所以很难找到与之对应的完全重构的滤波器。为了解决这个问题，Nick Kingsbury等提出了双树复小波变换(DTCWT)[7]，它按照一定的规则采用双树滤波的形式设计，既保留了一般复小波的优点，又可以完全重建，能够更有效的处理地震数据。

2 地震数据稀疏表示

2.1 压缩感知理论

压缩感知(Compressed Sensing,CS)理论表明，如果采集到的原始数据具有稀疏属性，同时使用一定的采样方法，就可以通过少量的欠采样数据来重建出具有一定精度的接近原数据的地震数据。从以上描述可以得知，压缩感知理论已在地震数据重建领域占有了一席之地。但对于绝大部分的地震图像来说，稀疏性不是一个共通的属性，所以需要对地震图像的稀疏表示，即进行某种正交变换，使其满足压缩感知理论中的稀疏特性。欠采样的地震数据重建过程可以表示为模型为式(1)。

y=Rf

(1)

来表述。式中:N维向量f表示原始的完整地震数据；R是一个M×N阶的采样矩阵，未采集的点的位置为0，其余采集点的值为1；N维向量y表示向量化后的不完整数据。由于采集数据的不完整性，R是一个欠定矩阵，整个地震数据重建的过程，就是由采样矩阵R和不完整的地震数据y重建出完整的、达到一定精度的地震数据f的过程。

根据压缩感知理论，如使用某种变换D∈RNXL使得向量a=DHf是稀疏的，则稀疏过程可以表示为式(2)。

y=RDa=Aa

(2)

其中A=RD为感知矩阵。因为向量a本身具有稀疏性，则原始重建问题转化为求解最小l0范式问题为式(3)。

min‖a‖0y=RDa

(3)

由于实际操作中允许误差的存在，则求解式 的最优化问题可以用一个最简单的近似求解形式代替为式(4)。

min‖a‖0‖y-RDa‖≤ε

(4)

其中ε为一个非常小的常量。l0范数最小使得稀疏非零元素的个数最少，而K

2.2 双树复小波变换

普通的复小波变换可以解决离散小波的局限性问题。但是由于一层分解以上的复小波变换的输入都是复数形式，所以很难找到与之对应的完全重建的滤波器。Nick Kingsbury等提出了双树复小波(DTCWT)解决了这个问题，它使用了双树滤波的形式进行设计，既保留了一般复小波的优点，又可以完全重建[6]。复小波可以表示为式(5)。

ψ(t) =ψr(t) +jψi(t)

(5)

ψr(t)，ψi(t)分别表示复小波的实部和虚部，他们都是实函数，这样双树复小波变换可以表示为两个独立的实小波变换，一个给出实部，一个给出虚部，如图1所示。

图1 双树复小波变换示意图

树a和树b为一对平行树，树a给出双树复小波变换的实部，树b给出双树复小波变换的虚部。双树复小波变换的优良性能来源于其解析性，为使ψ(t)满足解析性，则需要ψr(t)与ψi(t)组成Hilbert变换对，两个正小波函数组成Hilbert对的充要条件是：两低通滤波器满足半采样延迟条件，确保树b中的第一层向下采样取到树a中因隔点采样而舍弃的，未保留采样值。所以双树复小波具备了频谱单边性的优良特性，同时在二抽样条件下，让双树复小波具有其自身优势：

1) 平移不变性：双树复小波变换具有平移不变性，即信号的微小平移不会导致在各尺度上能量的变化。

2) 多方向选择性：双树复小波不单融合了离散小波所具有的良好视频特性，同时还有更好的方向分析手段。

3) 具有完全重构特性，地震图像数据在分解后可以保证完全重建，保证了数据的重建效果。

4) 数据冗余较为有限。对于一维信号其冗余为2∶1；对于二维信号冗余为4∶1。

5) 较少的计算量。分解重建过程的计算量比非抽象离散小波变换都要少很多。

基于以上所述的双树复小波所具有的优良特性，能够在地震数据重建领域有着很好的应用。

3 重建算法

双树复小波变换后，利用阈值软迭代算法(IST)[7]，进行地震数据的重建，同时考虑到相同尺度内的子带维数均相同，同位置的小波系数在同一尺度内的各个子带的位置是固定的，这样有利于分析同一位置小波系数间的统计特性。用合适的分布拟合小波系数分布，并通过分析系数间的相关性特征，建立合适的统计模型，可提高重建算法的精度。本文利用地震信号的双树复小波变换系数中同一方向子带与当前子带的父带小波系数之间的相关性，利用双阈值软迭代法[8,11]进行地震数据的最后重建步骤。

用a1、a2表示对采样前完整地震数据进行稀疏变换得来的小波系数及其父系数，u1、u2表示有缺失道地震数据稀疏变换后的小波系数及其父系数，ε1、ε2表示由于道缺失对地震数据造成的假频影响的小波变换系数及其父系数，则有式(6)。

u1=a1+ε1

(6)

双阈值法考虑的是当前系数与其父系数之间的关系，则有式(7)。

u2=a2+ε2

(7)

则上述式子可以改写为式(8)。

u=a+ε

(8)

(9)

根据条件概率公式进行推导，上式可变为式(10)。

(10)

在迭代重建过程中，注意到地震数据的道缺失具有随机性，则对于道缺失对于双树复小波变换后地震数据的影响ε，此处设ε服从均值为0，方差为σε高斯随机分布，即式(11)。

(11)

设原地震数据在同一方向的小波系数与其父系数的双变量联合概率密度函数为式(12)。

(12)

其中ɑ是待定参数，σ2表示被估计信号小波系数的方差。

则将上述式子代入，通过计算得到当前系数ɑ1的MAP估计为式(13)。

(13)

在以上计算推导过程中，为了计算得到当前系数的估计值 ，必须先得到模型中的未知参数 与被估计数据的方差 。在实际应用中，由于地震数据的却是造成的信号干扰的方差一般是未知的，通常需要进行估计。由于道缺失对于原始地震数据造成的影响主要集中在高频子带，且在各个高频子带中不尽相同且无法准确估计，此处，我们可以利用小波系数来进行中值估计 (Donoho,1994)[12]为式(14)。

(14)

其中表示第i个高频子带的双树复小波分解系数。

对于道缺失的地震数据的小波系数的方差 、被估计小波系数的方差 和假频小波系数方差 ，假设这三者之间满足关系式为式(15)。

(15)

对于方差 的估计值计算公式为式(16)

(16)

其中N是邻域U的大小。利用上式可以得到 的估计值为式(17)。

(17)

其中max表示原地震图像数据中的最大值，对于灰度图像的地震数据即255、MSE表示原图像与处理图像之间的均方误差。PSNR值越大，表示与原地震数据的差异越小。

双阈值收缩法插值重建由以下几个步骤组成：

1) 对欠采样地震数据进行双树复小波变换。利用正交小波变换的快速算法获得低分辨率下的尺度系数以及各个分辨率下的小波系数，其中尺度系数和小波系数共N个。

2) 对小波系数进行非线性阈值处理。此处利用不同频率间子带与父带的相关性，使用以上描述的双阈值算法进行重建，对分解过程中的低频系数不做处理。硬阈值虽然可以保留更多的信号特征，但是在平滑方面还有所欠缺。根据本文情况选用双阈值软收缩法。

3) 进行逆向双树复小波变换。由所有低频尺度系数、以及经由阈值处理后的小波系数做逆变换进行重建，提取欠采样的地震道，用插值法与欠采样地震数据插值重建为较为接近原地震数据的结果。

本文使用6层分解的双树复小波变换，需要逐层处理高频分解系数，每次插值对每层小波系数都进行双阈值处理，则整体的算法复杂度为o(N)。

4 实验结果及分析

实验的硬件平台采用双核CPU主频3.3G的Intel I5微机，内存容量为4G。系统软件为32位Windows7操作系统，仿真实验软件使用Matlab R2013b。实验数据采用广泛使用的marmousi模型，此模型由法国石油研究院于1998年做出，这种模型及其声波有限差分合成数据已被全世界成百上千的科研人员用于许多地球物理科研项目，直到今天仍然是出版最多的地球物理资料集之一。地震数据重建效果的衡量指标采用峰值信噪比(PSNR)，如式(18)。

(18)

其中MSE为不含噪声原始地震数据与去噪后地震数据的均方误差，定义为式(19)。

(19)

为了对比本文算法的性能，对采样地震数据分别进行传统的迭代阈值法与文中上述双阈值插值重建法进行重建。整个过程经过200次迭代，如图2—图7所示。

图2 迭代算法示意

图3 原始地震数据

图4 添加随机道缺失之后的地震数据，采样率为60%：(PSNR=20.829 8)

5 总结

由上述结果可知：本文算法相对于单阈值软迭代算法更为优秀，PSNR提升约为1.01，相对于欠采样的地震数据PSNR提升为3.6525。且从上述结果图中可以较为直接的看出，本文算法在地震数据重建中能够更好的还原地震数据的纹理与细节。需要指出的是，由于使用的双树复小波变换，相对传统离散小波变换其滤波器的构造更为复杂，且双阈值迭代每次迭代都需要与当前子层的父层一起进行迭代运算，需要更多的运算时间。

图6 使用双树复小波双阈值软迭代进行重建后的地震数据：(PSNR=24.482 3)

图7 双阈值迭代过程中PSNR的变化曲线