FK和Shearlet域联合压缩感知数据重构技术

闫海洋 周 辉 刘海波 徐朝红 孙赞东 刘 昭

(①中国石油大学(北京)油气资源与探测国家重点实验室，北京 102249； ②中国石油大学(北京)CNPC物探重点实验室，北京 102249； ③中国石油大学(北京)地球物理学院，北京 102249； ④东方地球物理公司海洋物探处，天津 300457)

0 引言

稀疏炮检点采集或野外采集因素造成地震数据的不规则影响后续地震资料成像质量。基于压缩感知理论的重构方法能够有效恢复地震数据。压缩感知数据重构技术是将重构问题转化为稀疏优化问题。由于空间随机缺失的地震道引起的空间假频在稀疏域呈白噪分布，有效信号在稀疏域集中分布，可采用反演迭代的方法逐步恢复缺失的地震数据。

Ma等[1]实现了基于Fourier变换的快速迭代阈值收缩算法，在Fourier域重构地震数据； Zhang等[2]在FK域数据重构过程中引入了凸集投影(POCS)算法，在迭代过程中引入了指数衰减阈值参数。Fourier变换是全局变换，不能描述地震数据的方向特征，且稀疏度不够，数据重构精度有待提高，于是能够描述地震信号方向特征的曲波变换应用到压缩感知数据重构中。刘国昌等[3]和Yang等[4]实现了曲波域的POCS方法地震数据重构，采用指数阈值提高了迭代收敛的速度； Zhang等[5]提出了一种基于曲波变换的三维数据同步重构和噪声压制方法； Zhang等[6]提出了一种基于非等距快速离散曲波变换的三维重构方法，采用线性Bregman迭代提高反演效率； Cai等[7]在曲波域基于L1范数的谱梯度投影(SPG-L1)算法实现了数据重构； Liu等[8]提出了局部随机采样的曲波域数据重构； Wang等[9]提出了一种基于曲波变换和高分辨率Radon变换的多域稀疏约束高精度地震数据恢复方法。

曲波变换采用网格旋转方式实现方向表征，无法实现对线奇异在连续系统和离散系统的统一，而Shearlet变换克服了这个缺点，采用剪切方式代替曲波变换的旋转方式，能更加稀疏表达地震信号。冯飞等[10]提出了采用Jitter欠采样的Shearlet变换稀疏约束地震数据重构方法； 张良等[11]实现了基于压缩感知的Shearlet变换地震数据重构，其重构效果优于Fourier变换、离散余弦变换、小波变换和曲波变换； 刘成明等[12]采用Shearlet变换与基于Landweber加速下降迭代法进行插值，在提高计算效率的同时保证了插值精度； 杨冠雨等[13]采用基于Shearlet变换的双正则化方法进行地震数据重构，同时兼顾了信号的稀疏性和地质结构的复杂性； 王常波[14]在Shearlet域利用快速凸集投影(FPOCS)算法和指数阈值进行压缩感知地震数据重构，与Fourier变换、曲波变换的重构结果对比表明，Shearlet变换重构精度最高。

与固定稀疏基不同，形态分量分析(Morphological Component Analysis，MCA)框架采用不同基函数线性组合方式对地震数据稀疏表达，通过迭代方式逐一求解并重建各形态分量，最后通过合并重建后的形态分量完成数据的重构。张凯等[15-16]提出了一种Shearlet字典和离散余弦变换(DCT)字典组合的地震数据重构方法； 周亚同等[17]提出在MCA框架下的基于DCT与曲波字典组合的地震信号重构方法。MCA框架的压缩感知数据重构需要迭代求解、重建各形态分量，计算效率有待提升。

本文采用Shearlet变换对地震数据稀疏表示，与常规对有缺失的地震数据进行Sheartlet变换不同，是将时空域的地震道重构转化为FK域的随机噪声压制问题。首先对有缺失的地震数据做FK变换，然后对FK域数据做Shearlet变换，通过阈值迭代法消除FK域由于地震道缺失引起的空间假频，从而实现地震道的时空域重构。与形态分量分析(MCA)框架下的地震数据重构不同，MCA认为信号由多个不同形态的分量线性组合而成，用不同的稀疏字典组合方式进行地震数据重构，本文采用的 FK变换加Shearlet变换，可以看作一种新的稀疏基变换，不需要迭代求解各形态分量。

稀疏采样方法是基于压缩感知观测系统设计的重要内容。全局随机采样空间假频呈白噪分布，分段随机采样空间假频具有蓝色噪声的频谱特征，在保持采样随机性的同时，有效解决了重采样结果疏密不匀问题，提高了地震数据的重构效果。最后通过Shearlet域和全局随机采样、Shearlet域和分段随机采样、FK+Shearlet域和全局随机采样、FK+Shearlet域和分段随机采样四种数据重构结果对比，表明FK+Shearlet域和分段随机采样结合重构效果最好，信噪比和保真度最高。

1 方法原理

1.1 FK+Shearlet域联合重构方法原理

考虑到传统采样定理的不足及信号的稀疏特性，Candes等[18]、Donoho[19]提出了压缩感知理论。压缩感知数据重构包括三个方面：①信号的稀疏表述，信号是稀疏的或在变换域内是稀疏的； ②测量矩阵的设计，测量矩阵要与稀疏基完全不相干； ③重构算法的构造，信号的有效恢复依赖高效高精度重构优化算法。

基于压缩感知的地震数据重构可表示为求解方程

d=Mf

(1)

式中：f∈Rp为完整数据；d∈Rq为观测数据；M为对角采样矩阵，对角线上的元素与d对应，采样前、后数据的维数分别为p和q，且p>q。当d中地震道数据未缺失时，对应的M对角线上的元素为1； 地震道数据缺失时，对应的M对角线上的元素为0。M一般选取高斯随机矩阵。

时空域的地震数据不是稀疏的，但可以通过数学变换满足稀疏特性。张良等[11]通过实验认为基于Shearlet变换的数据重构精度高于Fourier变换、DCT变换、小波变换和曲波变换。

本文利用FK+Shearlet变换对数据进行稀疏表示，将时空域不连续的地震数据重构问题转为连续FK域数据随机噪声压制问题，稀疏基为Shearlet基函数

x=φΨf

(2)

式中：Ψ为FK正变换；φ为Shearlet正变换；x为f的FK+Shearlet变换结果。

将式(2)代入式(1)，可得

d=MΨ-1φ-1x

(3)

式中：Ψ-1为FK逆变换；φ-1为Shearlet逆变换。令

θ=MΨ-1φ-1

(4)

为感知矩阵。

式(3)为欠定问题，无法直接求解。感知矩阵θ只有满足有限等距性质(RIP)时，则可以求取唯一解[18-20]。高斯随机矩阵与稀疏基不相干，θ满足RIP条件。地震数据重构问题可以表示为

(5)

式(5)在一定的条件下，L0范数最小化求解问题可以转化为L1范数最小化问题[22]，即

(6)

设ε为重构误差，式(6)可写为

(7)

关于高精度重构方法，国内外学者对Bregman迭代[6,21]、POCS[14，22-23]、SPG-L1[7，24-25]、正交匹配追踪法[26-27]、交替乘子方向算法(ADMM)[13，28]等方法进行研究或改进。温睿等[29]通过实验对比分析了POCS、IHT(Iterative Hard Thresholding)、Bregman和JRSI(Joint Reconstruction by Sparsity-promoting Inversion)四种迭代算法，认为POCS和IHT算法较为稳定，经过多次迭代能够达到较高的信噪比。本文采用POCS进行压缩感知数据重构，算法如下

(8)

(9)

1.2 稀疏采样

稀疏采样方法是压缩感知地震采集研究的重要内容，决定了由于地震数据缺失引起的空间假频在稀疏域的随机分布特征。相比全局随机采样，Jitter稀疏采样在保持随机性的同时，分布更均匀，减少了假频和有效信号的混叠，重建效果优于全局随机采样[10，12，30]。本文采用王汉闯等[31]提出的分段随机采样方法，该方法不仅具有更为优越的蓝谱特征，还突破了Jitter采样对段内样点为奇数的限制，具有更强的实用性。

以一道为例分析重采样因子对数据频谱(波数谱)的改造。对完整数据进行稀疏采样，就是对重采样因子与完整数据进行点乘，即

d(t)=r(t)f(t)

(10)

时间(空间)域的乘积对应频率(波数)域的褶积，即

D(ω)=R(ω)*F(ω)

(11)

式中：f为单道完整地震数据；r为重采样因子；d为采样结果；F、R和D分别为f、r和d的Fourier变换。完全采样时，r(t)≡1。

图1为规则、全局随机、分段随机1/2采样示意图。规则采样时采样因子r中0和1规则交替出现； 1/2全局随机采样的采样因子r中0和1全局随机出现； 1/2分段随机采样中2个采样点为1段，每段采样因子r中0和1随机出现。假定原始数据采样间隔为2ms，采样点数为3000，则奈奎斯特频率为250Hz。图2a为完全采样因子及其频谱，0频为有效信号频谱，频谱值为3000，不存在假频。图2b为1/2规则采样因子及其频谱，0频处为有效信号频谱，频谱值降为原来的1/2(1500)，假频在-250Hz处，与0频处有效信号频谱值相等。图2c为1/2全局随机采样因子及其频谱，0频处谱值降为原来的1/2(1500)，假频的频谱呈白噪特征，假频和有效信号混叠。图2d为1/2分段随机采样因子及其频谱，0频处谱值降为原来的1/2(1500)，假频的频谱呈蓝谱特征，假频和有效信号混叠。

图1 规则(a)、全局随机(b)、分段随机(c)1/2采样示意图实心表示采样； 空心表示未采样

图3a和图3b分别为完全采样和1/2规则采样的数据及其频谱。规则采样后，如果信号频率和假频混叠，则无法保真恢复原始信号。图3c和图3d分别为全局1/2随机采样和分段1/2随机采样的数据和频谱。随机采样后，假频振幅值相对较小，有效数据的频谱大于假频频谱，方便后续通过反演迭代进行全频带数据恢复，理论上可以突破奈奎斯特频率的限制。但分段随机采样数据的蓝谱特征明显，由于地震数据为带限数据，分段随机更有利于后期数据重构。

图2 不同稀疏采样因子(左)及其频谱(右)(a)完全采样因子； (b)1/2规则采样因子； (c)1/2全局随机采样因子； (d)1/2分段随机采样因子

1.3 FK+Shearlet变换

信号的稀疏表示就是将信号投影到正交变换基时，有效信号对应的变换系数集中分布，系数相对较大，噪声(对地震数据重构来说，主要指由于数据缺失造成的空间假频)对应的变换系数的绝对值很小。Shearlet变换[32-34]可以有效地表示地震数据各向异性特征。由剪切方式代替了曲波变换的旋转方式，能够在离散系统和连续系统做统一处理。

对于任意二维平方可积函数b(τ)，其Shearlet变换可表示为

SHφb{j,l,m}=〈b,φj,l,m〉

(12)

式中

(13)

(14)

其中：A为各向异性膨胀矩阵，控制尺度；B为剪切矩阵，控制方向，A和B均为2阶可逆矩阵；j为尺度参数；l为剪切参数；m为平移参数。

Shearlet变换是一类多尺度几何分析方法，通过对基函数的缩放、剪切和平移等仿射变换生成具有

图3 不同稀疏采样因子采样后的数据(左)及其频谱(右)(a)完全采样因子； (b)1/2规则采样因子； (c)全局1/2随机采样因子； (d)分段1/2随机采样因子

不同特征的剪切波函数，具有较好的地震数据的方向性描述和稀疏表达能力。

应用模拟单炮数据(采样间隔为2ms，512个样点，512道)对比Shearlet变换与FK+Shearlet变换效果。采用1尺度的Shearlet分解，将模拟单炮分解成1个不带方向特征的低频分量和8个带方向特征的高频分量。图4为模拟单炮及FK谱，图5为对模拟单炮进行Shearlet分解的9个系数，由于分解的8个高频系数(图5a～图5h)能量较弱，对其振幅放大10倍显示。

图6为图5的9个Shearlet域系数对应的时空域数据。图7为对模拟单炮进行FK+Shearlet分解的9个系数谱，图8为对模拟单炮进行FK+Shear-let分解的9个系数对应的时空域数据。表1对比了模拟单炮不同稀疏基变换分解系数能量占比。对模拟单炮直接进行Shearlet变换，不带方向特征的低频分量系数能量占比为98.75%，8个带方向特征的高频分量系数能量占比仅为1.25%； 对模拟单炮进行FK+Shearlet变换，不带方向特征的低频能量占比为8.81%，8个带方向特征的高频分量系数能量占比为91.19%，FK+Shearlet域的带有方向特征的高频系数相对Shearlet域能量更强，能更好地观测地震数据的方向(线奇异)特征。

图4 模拟单炮(a)及其FK谱(b)

图5 模拟单炮Shearlet分解的9个系数(a)～(h)带方向特征的8个高频系数谱，振幅放大10倍显示； (i)不带方向特征的低频系数谱

图6 图5的9个系数对应的时空域数据(a)～(h)带方向特征的8个高频系数谱，振幅放大10倍显示； (i)不带方向特征的低频系数谱

图7 模拟单炮分解的FK+Shearlet域 9个系数谱(a)～(h)带方向特征的8个高频系数谱； (i)不带方向特征的低频系数谱

图8 图7的9个系数谱对应的时空域数据(a)～(h)带方向特征的8个高频系数谱； (i)不带方向特征的低频系数谱

表1 模拟单炮不同变换分解的系数能量占比

2 模拟数据实验

采用重构信噪比评价重构效果

(15)

式中f*表示重构的地震数据。

应用上述模拟单炮记录对比全局随机和分段随机采样250道(采样率48.8%)在Shearlet域与FK+Shearlet域的数据重构效果。表2 为模拟单炮重构信噪比统计。

图9为模拟单炮记录的全局随机采样结果及其FK谱。图10为Shearlet域和FK+Shearlet域数据重构结果及误差。Shearlet域重构数据(图10a)的信噪比为13.39dB，FK+Shearlet域重构数据(图10c)的信噪比为21.05 dB； Shearlet域数据重构的误差(图10b)明显大于FK+Shearlet域数据重构误差(图10d)，表明FK+Shearlet域重构效果优于Shearlet域。全局随机采样数据FK谱(图9b)的空间假频呈白噪分布，FK+Shearlet域重构数据有效信号的FK谱(图11b)能量的恢复及空间假频的压制效果明显优于Shearlet域重构数据的FK谱(图11a)。

图9 模拟单炮全局随机采样结果(a)及其FK谱(b)

图10 模拟单炮全局随机采样Shearlet域与FK+Shearlet域模拟单炮重构效果分析(a)Shearlet域重构结果； (b)Shearlet域重构误差； (c)FK+Shearlet域重构结果； (d)FK+Shearlet域重构误差

图11 模拟单炮全局随机采样不同域重构结果的FK谱(a)Shearlet域； (b)FK+Shearlet域

图12为模拟单炮记录的分段随机采样结果及其FK谱。图13为Shearlet域和FK+Shearlet域数据重构结果及误差。Shearlet域重构数据(图13a)的信噪比为13.89dB，FK+Shearlet域重构数据(图13c)的信噪比为23.52dB，Shearlet域数据重构的误差(图13b)大于FK+Shearlet域数据重构误差(图13d)，FK+Shearlet域重构效果优于Shearlet域。分段随机采样数据的FK谱(图12b)的空间假频呈蓝谱特征，空间假频与有效信号混叠相对全局随机(图9b)要弱。FK+Shearlet域重构数据有效信号的FK谱(图14b)的能量恢复及空间假频的压制效果优于Shearlet域重构数据的FK谱(图14a)。从数据重构信噪比(表2)可见，Shearlet域分段随机重构效果(图13a)优于全局随机重构(图10a)，FK+Shearlet域分段随机重构效果(图13c)优于全局随机重构(图10c)。

通过不同稀疏采样比例进一步验证FK+Shearlet的重构效果。图15为分段随机采样下不同稀疏采样比例和不同稀疏基变换模拟单炮重构信噪比对比。整体上，FK+Shearlet比Shearlet域重构信噪比高。当稀疏采样比例低于60%时，FK+Shearlet比Shearlet域重构信噪比高6dB以上； 当采样比例高于70%时，随着采样比例的增加，两种稀疏基重构效果逐渐接近。

图12 模拟单炮分段随机采样结果(a)及其FK谱(b)

图13 模拟单炮分段随机采样Shearlet域与FK+Shearlet域重构效果分析(a)Shearlet域重构结果； (b)Shearlet域重构误差； (c)FK+Shearlet域重构结果； (d)FK+Shearlet域重构误差

图14 模拟单炮分段随机采样不同域重构结果的FK谱(a)Shearlet域； (b)FK+Shearlet域

表2 模拟单炮不同域重构信噪比统计

图15 不同稀疏采样比例FK+Shearlet(红色)与Shearlet域(蓝色)重构数据信噪比对比

3 实际数据应用

选择海域地震资料(图16a)进行压缩感知数据重构。由于海洋资料频带较宽，同时原始数据也存在空间采样不足引起的空间假频，给数据重构带来一定难度，以70%全局随机采样(图16b)和分段随机采样(图16c)对比Shearlet域与FK+Shearlet域数据重构效果。图17为原始、全局随机采样、分段随机采样数据的FK谱，分段随机采样数据FK谱(图17c)的空间假频呈蓝谱特征，空间假频与有效信号混叠相对全局随机(图17b)要弱。表3为不同稀疏采样和稀疏基变换实际数据重构信噪比对比。

图16 实际原始地震数据(a)及其全局(b)和分段随机采样结果(c)

图17 实际原始地震数据(a)及其全局(b)和分段随机采样结果的FK谱(c)

在全局随机采样情况下，Shearlet域重构数据(图18a)的信噪比为10.08dB，FK+Shearlet域重构数据(图18b)的信噪比为15.66dB，Shearlet域数据重构的误差(图18c)明显大于FK+Shearlet域数据重构误差(图18d)，FK+Shearlet域重构效果优于Shearlet域重构效果。在分段随机采样情况下，Shearlet域重构数据(图19a)的信噪比为10.47dB，FK+Shearlet域重构数据(图19b)的信噪比为16.52dB，Shearlet域数据重构的误差(图19c)大于FK+Shearlet域数据重构误差(图19d)，FK+Shearlet域重构效果优于Shearlet域重构。

在全局随机采样情况下，FK+Shearlet域重构数据有效信号的FK谱(图20b)能量的恢复及空间假频的压制效果明显优于Shearlet域重构数据的FK谱(图20a)。在分段随机采样情况下，FK+Shearlet域重构数据有效信号的FK谱(图20d)能量恢复及空间假频的压制效果优于Shearlet域重构数据(图20c)。

图18 实际数据全局随机采样下不同稀疏域重构结果对比(a)Shearlet域重构结果； (b)FK+Shearlet域重构结果； (c)Shearlet域重构误差； (d)FK+Shearlet域重构误差

图19 实际数据分段随机采样下不同稀疏域重构结果对比(a)Shearlet域重构结果； (b)FK+Shearlet域重构结果； (c)Shearlet域重构误差； (d)FK+Shearlet域重构误差

图20 实际数据两种采样方式、不同域重构结果的FK谱对比(a)全局随机，Shearlet域； (b)全局随机，FK+Shearlet域； (c)分段随机，Shearlet域； (d)分段随机，FK+Shearlet域

由表3的数据重构信噪比可见，Shearlet域分段随机重构效果(图19a)优于全局随机重构(图18a)，FK+Shearlet域分段随机重构效果(图19b)优于全局随机重构(图18b)。

表3 实际数据的重构信噪比统计

4 结论

压缩感知数据重构技术要求地震数据随机稀疏分布，由于数据缺失引起的假频在稀疏域表现为噪声，有效信号集中分布，假频和有效信号混叠，但是假频相对较弱，通过迭代的方法能进行数据重构。分段随机采样方法具有更为优越的“蓝色噪声”频谱分布，分段随机采样在保持采样随机性的同时，有效解决了重采样结果疏密不均问题。本文采用FK+Shearlet变换作为稀疏基进行分段随机采样数据重构。实验表明，FK+Shearlet域比Shearlet域能更好地描述地震数据的方向(线奇异)特征，分段随机采样FK+Shearlet域重构精度高于全局随机采样Shearlet域重构、分段随机采样Shearlet域重构和全局随机采样FK+Shearlet域重构。