空间自适应EPLL低频地震随机噪声降噪方法

林红波，马 阳

(吉林大学 通信工程学院，吉林 长春 130012)

在利用地震勘探对地下结构进行成像与分析时，提取地震图像蕴含的地下信息是关键。然而，在接收地下反射信号时，不可避免地会接收到人文活动和自然环境信息，如风的运动引起的随机噪声。受随机噪声污染的实际地震勘探图像Y=X+V，其中随机噪声V的不可预测性和复杂的噪声性质严重影响了信号X的提取精度[1]。在沙漠地区，随机噪声受复杂地质和地表条件影响，其波形、空间结构与反射信号有一定的相似性，是一种低频色噪声，这种性质使得低频反射信号的辨识和提取更为困难。此外，地震勘探图像中非平稳地震勘探信号的强度和频率在时间和空间上变化剧烈。首先到达的信号能量强，在图像中占主导地位，来自深层反射的信号能量弱，其在强信号和弱相似随机噪声的干扰下很难分辨。如何压制非平稳地震勘探图像中的低频弱相似随机噪声，且不损失非平稳地震信号，一直是沙漠地震勘探研究热点和难点。

国内外地球物理学者和信号处理学者在地震勘探去噪方面做了大量研究，并不断提出新的去噪方法以适应地震勘探技术的发展和地震勘探目标的新要求。从最初经典的傅里叶变换发展到适应非平稳信号的短时傅里叶变换[2]，从小波变换到Curvelet 变换[3]、脊波变换[4]，从时频峰值滤波[5]到径向道时频峰值滤波[6]等一系列方法抑制地震图像中的噪声，并取得了显著的去噪效果。然而，这些降噪方法未能充分挖掘地震勘探信号复杂的结构特征，导致在压制噪声的同时有效信号也存在一定程度丢失。

近年来，基于字典学习方法有效改进了图像降噪效果。基于字典学习的去噪模型通过更新字典学习图像中信号的结构特征，能够更好地抑制图像中的随机噪声。AHARON等[7]提出一种K-奇异值分解(K-Singular Value Decomposition，K-SVD)的字典学习方法来稀疏地描述图像特征。YANG等[8]结合自定进度的学习实现基于非负矩阵分解的非负系数稀疏字典学习方法，在时空域中抑制地震图像中的随机噪声。王勇等[9]结合字典学习思想，提出一种冗余自适应字典表示的稀疏贝叶斯学习算法。从待重建图像的迭代中间图像块中学习字典，并以该字典作为基获得图像的稀疏解。刘帅等[10]在分层字典学习框架下进行统计建模，通过引入高斯噪声项和稀疏噪声项，有效地表达图像噪声特性；同时，应用吉布斯采样求解统计模型，以实现降噪的目的。字典学习的方法通常利用图像块学习字典，减小字典学习的复杂度，然而却忽略了图像的全局特征。针对此问题，ZORAN等[11]提出了块期望对数似然降噪方法(Expected Patch Log Likelihood，EPLL)，利用高斯混合模型(Gaussian Mixture Model，GMM)学习目标信号的结构特征作为先验约束去噪过程。该算法既能精准描述图像块统计特征，又能从强噪声中无失真地复原自然图像。PAPYAN等[12]进一步提出多尺度EPLL方法，改善算法的降噪性能。SULAM等[13]利用稀疏先验替代高斯混合模型，提出稀疏期望块对数似然方法来恢复自然图像。基于字典学习的去噪方法也在地震图像噪声抑制中表现出优秀的地震信号恢复能力。XI等[14]结合块分类和高斯混合模型提出基于块分类的期望块对数似然算法(Patch Classification EPLL，PC-EPLL)，实现非平稳地震随机噪声干扰下的信号恢复。基于EPLL对地震特征良好的恢复能力，笔者进一步将其拓展到非平稳地震勘探图像，压制低频有色随机噪声。

针对地震勘探信号非平稳特性，笔者提出空间自适应块期望对数似然(Spatially Adaptive EPLL，SA-EPLL)低频地震随机噪声降噪方法。在EPLL图像降噪方法中，正则化参数通常设为与随机噪声方差有关的一组常数。然而，当地震勘探图像中信号强度随空间非平稳变化时，具有相同噪声水平的图像块具有不同的信噪比，所需要的去噪强度不同。因此，EPLL算法正则化参数仅随噪声方差变化而无法实现对不同信噪比图像块施加不同的去噪水平，导致EPLL算法对非平稳地震勘探信号去噪效果不佳。为了同时表征地震勘探信号的非平稳特征和噪声特性，笔者提出利用更能反映非平稳特征的块信噪比自适应地控制正则化参数，构建基于块信噪比的加权图像重构方案。合成图像和实际地震勘探图像验证了空间自适应EPLL算法有效性，SA-EPLL能够有效压制沙漠地震勘探非平稳图像中低频弱相似随机噪声，缓解了非平稳地震信号的损失问题。

1 期望块对数似然降噪方法算法

EPLL算法通过学习样本图像块的结构特征作为先验，约束图像去噪反问题。对于含噪图像Y=X+V，假设V是均值为0、方差为σ2的高斯白噪声。利用EPLL方法对Y去噪，即通过求解如下反问题估计干净图像X，代价函数表示为

(1)

其中，logP(RiX)表示第i个重叠图像块的对数似然先验，Ri为抽取第i个重叠图像块的算子，重叠图像块大小为q×q。参数λ为常数，通常设置为q2/σ2。在对数似然先验约束下的优化问题不能直接求解，ZORAN和WEISS采用半二次分裂算法[11]，通过引入辅助变量zi，将式(1)的优化问题改写为

(2)

(3)

给定X，求解式(2)关于zi的最小化问题。由于对数似然先验约束下没有闭合解，在实际应用中，利用已经训练好的GMM先验模型获得近似解。从GMM模型中搜索与含噪图像块zi最似然的GMM分量

(4)

(5)

基于EPLL降噪原理，正则化参数β通常根据整体噪声方差σ2设定，由于在迭代过程中图像的噪声方差已经改变，在图像去噪应用中通常设定β为一组常数序列[1，2，4，8，…]/σ2。

将EPLL应用到沙漠地震勘探图像时，正则化参数β的设置需要考虑地震勘探信号的非平稳特性。地震图像是由多个检波器接收到的地下反射信号构成的。受地下介质吸收影响，地震图像中信号强度在空间和时间上变化。随着时间增加，地震信号的幅度迅速降低，强信号与弱信号的强度差异巨大；由于检波器间隔距离远，不同检波器接收到的信号强度也不均衡。因此，地震图像下方的弱信号在上方强信号的掩盖下很难辨识。另一方面，沙漠地震勘探图像中的随机噪声为低频色噪声，噪声特性较高斯白噪声更为复杂，甚至部分沙漠勘探随机噪声与地震信号具有相似的波形和形态。这种复杂的随机噪声进一步增加了沙漠勘探图像去噪的难度。当EPLL算法应用于沙漠地震图像时，地震信号的非平稳性对降噪效果会产生较大的影响。

以图1(a)为例说明对非平稳地震图像的影响。图像块中信号幅度分别为[0.2，0.8，1，10，1000]，含噪图像块(图1(b))中沙漠噪声方差为0.08；可见，沙漠噪声结构与地震信号相似。根据EPLL原理将正则化参数β设为常数序列，去噪结果如图1(c)所示，信号幅度为0.2的图像块噪声压制效果不明显，幅度为10和1 000的图像块噪声压制良好，但左下角信号出现变形。幅度为0.8与1的图像块降噪效果好。上述结果表明，正则化参数仅与噪声方差有关时，对于非平稳信号图像块滤波效果不同。含有大幅度信号的图像块信噪比高，噪声方差相对于信号来说相当小，若EPLL采用噪声方差设置正则化参数，则对高信噪比图像块相当于施加了过强的去噪，导致信号失真；对于小幅度信号图像块，信噪比极低，噪声方差相对于信号来说相当大，此时利用噪声方差设置正则化参数相当于降低了EPLL降噪能力，图像块去噪不理想。由此可见，针对非平稳地震图像处理，不能简单采用噪声方差调节EPLL参数。为了衡量信号与噪声的相对信息，笔者提取图像块的信号方差，利用信号方差与噪声方差的比值设置正则化参数。图1(d)所示EPLL去噪结果表明，不论是弱信号图像块还是强信号图像块，去噪结果和信号保持程度均得到改善。上述分析表明，对于非平稳地震图像，采用反映块信噪比的特征设置正则化参数更为合理。

图1 不同正则化参数的非平稳图像块去噪效果分析

2 基于空间自适应EPLL的沙漠地震图像去噪

在对非平稳沙漠地震勘探图像去噪时，若仅考虑噪声方差水平设置EPLL正则化参数，要保证强信号图像块不失真，则会导致弱信号块噪声压制不彻底，很难同时兼顾强信号保真和弱信号恢复。这就需要EPLL正则化参数根据每个图像块特征调节。笔者在EPLL框架下，将能反映信号非平稳特性的块信号方差与反应噪声水平的参数β相结合，提出空间自适应块期望对数似然(SA-EPLL)去噪算法。该算法利用块信号方差调节参数β得到每个块的信噪比，并通过块信噪比实现对强信号和弱信号图像块施加不同的正则化参数，自适应控制图像块去噪水平。同时利用块信噪比加权重构图像，平衡局部先验和全局特征，改进非平稳图像去噪效果。

2.1 空间自适应EPLL

(6)

(7)

(8)

(9)

(10)

SA-EPLL流程如算法1所示。

算法1：空间自适应EPLL算法(SA-EPLL)。

输入：含噪地震图像Y、学习到的GMM先验。

设t为当前迭代次数，开始迭代。

② 求解平稳化图像去噪结果。

2.2 训练集构造

在使用SA-EPLL算法对沙漠地震图像去噪时，需要利用GMM模型学习地震图像先验，因此，构建包含实际地震图像结构特征的训练集对于恢复地震图像尤为重要。笔者采用雷克子波、零相位波、混合相位波[15]模拟实际地震图像中的信号特征，生成训练集。为了防止在学习过程中过拟合大的数据特征，笔者构建均值为0，方差为1的训练集，使学习到的模型更加鲁棒。利用上述构建的训练集，笔者用EM算法[16]学习高斯混合模型参数，对地震信号结构特征建模。所学习的混合GMM模型包含200个高斯分量，在保证特征丰富的同时，又不会增加计算量。

在先验字典学习中，图像块的大小是影响EPLL降噪效果的重要因素。在地震勘探图像处理时，地震勘探同相轴结构信息对于精准探查地下结构至关重要。要在图像块中体现出倾斜同相轴的空间相干性，需要较大的图像块，但采用太大的图像块尺寸，会导致学习时间迅速增加。为了尽可能多地保留地震勘探信号的结构特征，避免学习时间过长，笔者分析大量实验降噪效果寻找最优值，将图像块大小设置为30×30。

3 实验及分析

为了验证SA-EPLL对沙漠图像中的随机噪声抑制能力，笔者将其应用于合成地震勘探图像和实际地震图像，并与BM3D、EPLL和多尺度EPLL去噪方法进行比较。

3.1 合成沙漠地震图像去噪

笔者利用合成沙漠地震勘探图像验证空间自适应EPLL算法压制沙漠勘探中随机噪声的能力。所采用的合成沙漠地震图像如图2(a)所示，包含5个同相轴，其幅度从上到下分别设为[1.8，1.6，1.4，1.2，1，0.8]。通常假设地震勘探系统为线性系统，噪声与信号独立，随机噪声为加性噪声。在合成沙漠图像中加入模拟的沙漠随机噪声，得到的含噪图像(图2(b))的信噪比为-4.2 dB。图2(c)为从含噪图像中计算出的块信噪比，可见块信噪比反映了地震信号强度变换和结构特征。利用块信噪比调节SA-EPLL正则化参数能够自适应控制非平稳图像去噪强度。

图2 合成沙漠地震图像

笔者分别采用BM3D、EPLL、多尺度EPLL、SA-EPLL方法处理合成沙漠地震图像。EPLL、多尺度EPLL和 SA-EPLL算法均采用步长为15重叠分块，图像块大小为30×30。在实验中设置ε为0.001，经多次实验，表明迭代次数为6次时具有好的去噪效果，增加迭代次数会导致过度去噪，引起信号失真。在实验中，当相邻两次迭代结果的谱范数小于ε或者迭代次数达到6次时，停止迭代。

4种方法的滤波结果如图3所示。分析滤波结果可知，在经过BM3D方法、EPLL方法与多尺度EPLL处理之后，仍然残留着大量的噪声，尤其在EPLL和多尺度EPLL去噪结果中噪声残留十分明显。且从滤波结果和干净图像的残差图中可以看出，经过BM3D方法、EPLL方法与多尺度EPLL方法去噪后，残差地震图像中仍然可以看到同相轴，说明这3种方法在滤波同时造成有效信号的大量丢失。笔者提出的方法背景更为干净，沙漠随机噪声去除的更彻底，同相轴也变得更清晰，在残差图像中几乎没有有效信号残留。由此可知，笔者提出的方法在有效压制噪声的同时，尽可能地保留了非平稳地震信号。

图3 不同方法去噪结果及对应残差图

为了进一步评估去噪方法，笔者改变模拟噪声强度加入合成地震勘探图像，分别使用BM3D、EPLL、多尺度EPLL以及SA-EPLL算法对其去噪，计算滤波前后的信噪比(SNR)，由r表示，和均方误差(MSE)，由e表示

(11)

(12)

式中，X为干净的合成信号，Q是地震图像中的采样点数。表1给出了4种方法的信噪比和均方误差。从表1可以看出，笔者提出的方法获得最高的信噪比及最小的均方误差，降噪效果优于BM3D、EPLL和多尺度EPLL算法。

表1 不同方法去噪结果信噪比与均方误差比较：SNR(dB)/MSE

3.2 实际沙漠地震图像去噪

验证笔者提出的方法对实际地震勘探图像的去噪能力，并与BM3D方法、EPLL以及多尺度EPLL方法进行比较。所采用的实际地震勘探图像来源于中国塔里木盆地，包含200道记录，抽样频率为1 000 Hz，道间距20 m，信噪比为-3.29 dB，如图4所示。从实际沙漠地震勘探图像中可以看出，信号能量随着时间的增加而衰减，在直达波区域信号占主导地位，在反射波区域噪声逐渐占主导地位。

图4 实际沙漠地震图像

笔者分别采用4种滤波方法对实际沙漠地震图像处理。由图5可知，BM3D方法对沙漠噪声的压制能力很弱，信号没有被恢复，信噪比为 2.46 dB。EPLL与多尺度EPLL方法能够滤除部分沙漠噪声，但噪声残留明显，信噪比分别为3.75 dB与4.77 dB。笔者提出的SA-EPLL方法能够抑制实际沙漠地震图像中的大部分沙漠噪声，信噪比提高到9.18 dB。对比图5中的白框区域发现，笔者提出的方法在去噪图像中背景更干净，有效同相轴变得更光滑，明显优于对比的方法。

图5 实际沙漠地震图像及不同方法去噪结果

4 结束语

笔者提出空间自适应EPLL低频地震随机噪声降噪方法。该方法针对地震图像中信号非平稳问题，利用图像块的统计特征，控制EPLL算法中的正则化参数，构建了空间自适应EPLL算法。该算法实现了根据块信噪比自适应控制每个块的正则化参数，调节对每个块的去噪强度，进而提升块去噪效果。此外，在重构过程中以块信噪比对图像块加权平均，以避免非平稳地震信号的损失。合成和实际地震图像去噪结果表明，笔者提出的空间自适应EPLL方法在保留非平稳地震信号和抑制低频随机噪声方面优于其他竞争性降噪方法。