变换域噪声压制是指利用数学变换找出某种基函数来稀疏表达有效信号, 使有效信号能量主要分布在少数较大的系数上, 而噪声分布在较小的系数上, 从而实现信噪分离。如20世纪80年代初发展起来的小波变换域噪声压制方法^[1-10], 利用“正方形”基函数来逼近信号, 虽然容易检测出边缘奇异点, 但无法准确表示边缘点间的连续性, 且缺乏方向性, 不能很好地表达图像中的曲线或面奇异。DO等^[11]在小波变换的基础上发展了Contourlet变换, 用“长条形”基函数来表达信号, 其支撑区间的宽度与长度近似符合平方律, 具有优越的非线性逼近能力, 在图像增强、图像去噪、模式识别、统计分析等领域得到了广泛应用。彭才等^[12]将Contourlet变换应用于随机噪声衰减, 采用硬阈值法修改Contourlet系数, 反变换重建原信号。张恒磊等^[13]采用相关叠加和非线性阈值相结合的方法对Contourlet系数进行修改, 以压制随机噪声。王建花等^[14]利用同一方向子带内相邻Contourlet系数之间的相关性, 对每一个系数计算邻域均方根振幅, 定义该系数与邻域均方根振幅比值的平方为邻域相关系数, 根据邻域相关系数定义阈值函数自适应地衰减随机噪声。

K-L变换最初应用在遥测多谱资料处理中。HEMON等^[15]1978年首次将K-L变换应用到地震资料处理领域, 基于随机噪声与有效信号的相干性差异去除随机噪声和相干噪声。其后, 傅才芳等^[16]应用K-L变换去除地震资料中的随机噪声。杨文采^[17]用非线性K-L变换去除海洋地震资料中的随机噪声。张广智等^[18]先将地震数据分块并拉平相干同相轴, 再利用K-L变换进行随机噪声衰减。LIU^[19]先选出面波区域并将面波拉平, 再用K-L变换提取面波并从原始信号中减去面波, 取得较好的去噪效果。

地震资料经Contourlet变换后, 有效反射信号的Contourlet系数较大并具有较强的相关性, 随机噪声则均匀分布于变换域, 系数较小且相关性弱, 对Contourlet系数的处理方式决定了Contourlet变换的去噪效果。本文充分利用有效信号的Contourlet系数在同一方向子带内及不同方向子带之间的相关性和非高斯性, 结合K-L变换在分类特征提取上的优势, 对Contourlet系数进行K-L变换以进一步分离有效信号和随机噪声。利用最大似然估计法^[20-21]和多尺度噪声估计法对有效信号和随机噪声的Contourlet系数进行方差估计, 并定义K-L变换域能量百分比阈值函数, 自适应地确定K-L反变换的特征向量, 修改Contourlet变换的系数, 再通过Contourlet反变换压制随机噪声。

1 方法原理 1.1 Contourlet变换

Contourlet变换也可称作小轮廓变换, 采用分段光滑的基函数对原始信号进行逼近。Contourlet变换可分为两个独立的步骤:①根据信号的特点, 采用拉普拉斯金字塔滤波器^{[11, 22]}(Laplacian Pyramid, LP)对二维信号进行多尺度分解, 以捕获不同尺度上信号的边缘奇异点; ②通过方向滤波器组^[23-25](directional filter bank, DFB)对LP分解后的每一尺度高频信号进行方向分解, 并将同方向上的奇异点连成线, 合并为同一系数, 即Contourlet系数。

1.1.1 拉普拉斯金字塔分解

拉普拉斯金字塔分解源于影像的高斯金字塔分解, 即将平滑滤波与下采样操作结合获得原始图像的近似分量, 对此近似分量重复进行平滑滤波与下采样操作, 则可得到一系列不同尺度的高斯图像金字塔, 相邻两层高斯金字塔图像相减则为拉普拉斯金字塔图像。故拉普拉斯金字塔分解过程如下:

1) 通过低通滤波和下采样操作得到原始信号x的近似分量, 也就是原始信号的低频子带信号;

2) 对低频子带信号进行上采样和高通滤波, 生成原始信号的预测信号d;

3) 将原始信号x与预测信号d相减, 得到原始信号的高频子带信号c。

此后利用每一步分解所产生的低频子带信号, 迭代进行LP分解, 生成一个低频信号和一系列高频子带信号。分解及重构过程如图 1所示, 其中, H为(低通)分解滤波器, G为(高通)合成滤波器, M为采样矩阵。拉普拉斯金字塔滤波实现了信号从精细尺度到粗略尺度的逼近。

1.1.2 方向滤波器组

方向滤波器组DFB可以通过l级树型分解来实现。在每一级中, 频域被划分为2^l个楔型的方向子带^[25]。DFB的实现分以下两步:①将扇型滤波器与梅花型采样矩阵相结合, 实现对原始信号二维频谱的切分, 即第一与第二级分解; ②从第三级分解开始, 将扇型滤波器与幺模采样矩阵相结合, 实现局部方向信息的提取。图 2为前两级DFB分解的频率分割示意图^[16]。第一级和第二级分解采用梅花型采样矩阵Q₀, Q₁:

$ {\mathit{\boldsymbol{Q}}_0} = \left[ {\begin{array}{*{20}{c}} 1&{ - 1}\\ 1&1 \end{array}} \right] $

(1)

$ {\mathit{\boldsymbol{Q}}_1} = \left[ {\begin{array}{*{20}{c}} 1&1\\ { - 1}&1 \end{array}} \right] $

(2)

第三级及以后各步采用幺模采样矩阵, 分别为:

$ {\mathit{\boldsymbol{R}}_0} = \left[ {\begin{array}{*{20}{c}} 1&1\\ 0&1 \end{array}} \right] $

(3)

$ {\mathit{\boldsymbol{R}}_1} = \left[ {\begin{array}{*{20}{c}} 1&{ - 1}\\ 0&1 \end{array}} \right] $

(4)

$ {\mathit{\boldsymbol{R}}_2} = \left[ {\begin{array}{*{20}{c}} 1&0\\ 1&1 \end{array}} \right] $

(5)

$ {\mathit{\boldsymbol{R}}_3} = \left[ {\begin{array}{*{20}{c}} 1&0\\ { - 1}&1 \end{array}} \right] $

(6)

DFB的第一级将原始信号经扇形滤波器分解为水平和垂直两个方向(如图 2a), 第二级将扇形滤波器与梅花型采样矩阵Q₀交换形成如图 2b所示的象限滤波器, 将第一级与第二级滤波器组合后进行Q₁采样, 则可获得如图 2c所示的原始信号4个方向子带。从第三级开始所用的滤波器组是扇形滤波器与幺模采样矩阵形成的平行滤波器组, 实现信号的2^l级分解。方向滤波器组会将信号的低频分量泄漏于几个方向子带中, 在应用之前需要将信号的低频部分移除, 因此将LP分解模块与DFB分解模块结合, 形成了双层滤波器组结构的Contourlet变换。图 3为Contourlet变换基本流程^[26], 图 4为Contourlet变换4级分解系数。

1.2 基于Contourlet变换的K-L变换随机噪声衰减方法

图像经Contourlet变换稀疏表达后, 采用边缘统计模型及联合统计模型对变换系数进行研究表明:Contourlet变换系数具有聚类性和持续性, 以及非高斯、高峰度和长拖尾等特点^[27]。这保证了图像中有效信息主要集中于变换域的有限区域内, 且系数较大, 而高斯性的随机噪声经Contourlet变换后仍为白噪声, 其均匀分布于变换域且幅值较小^[28]。同时Contourlet系数之间存在一定的相关性, 即如果一个Contourlet系数幅值大, 则与其相邻的Contourlet系数很有可能是“大”系数。

地震资料中的有效信号具有丰富的纹理和边缘特征, 而随机噪声符合高斯白噪声的特点。有效信号的Contourlet变换域系数较大且相互关联, 而随机噪声的Contourlet变换域系数较小且均匀分布于各子带上, 无明显相关性。利用这些特性修改Contourlet系数可实现对随机噪声的压制。然而, 任何一种变换都不能达到信号与噪声的完全分离, 有必要进一步扩大有效信号和随机噪声的分离度。

K-L变换具有分类特征提取的优势, 它对输入信号按其协方差矩阵的归一化特征向量进行正交分解。经K-L变换后, 存在相干性的有效信号集中在前面几个特征值对应的特征向量上, 随机噪声则分布在最后几个特征值对应的特征向量上, 反变换时去除代表随机噪声的特征向量, 则可将相干性好的有效信号保存下来, 压制随机噪声, 提高地震资料的信噪比。本文对Contourlet系数进行K-L变换, 以最大限度地将有效信号与随机噪声分离。

K-L变换原理^[29]如下。

对于第i层第l个方向子带的Contourlet系数矩阵

$ \mathit{\boldsymbol{C = }}\left[ {\begin{array}{*{20}{c}} {{c_{11}}}&{{c_{12}}}& \cdots &{{c_{1n}}}\\ {{c_{21}}}&{{c_{22}}}& \cdots &{{c_{2n}}}\\ \vdots&\vdots &{}& \vdots \\ {{c_{m1}}}&{{c_{m2}}}& \cdots &{{c_{mn}}} \end{array}} \right] $

(7)

其协方差定义为:

$ {\rm{Cov}}\left( \mathit{\boldsymbol{C}} \right) = E\left\{ {\left[ {\mathit{\boldsymbol{C}} - E\left( \mathit{\boldsymbol{C}} \right)} \right]{{\left[ {\mathit{\boldsymbol{C}} - E\left( \mathit{\boldsymbol{C}} \right)} \right]}^{\rm{T}}}} \right\} $

(8)

式中:E(C)表示C的数学期望, [C－E(C)]^T为[C－E(C)]的转置。Cov(C)为m×m的实对称正定矩阵, 存在一个正交矩阵U, 使得

$ {\mathit{\boldsymbol{U}}^{\rm{T}}}{\rm{Cov}}\left( \mathit{\boldsymbol{C}} \right)\mathit{\boldsymbol{U}} = \left[ {\begin{array}{*{20}{c}} {{\lambda _1}}&0& \cdots &0\\ 0&{{\lambda _2}}& \cdots &0\\ \vdots&\vdots &{}& \vdots \\ 0&0& \cdots &{{\lambda _m}} \end{array}} \right] = \mathit{\boldsymbol{ \boldsymbol{\varLambda} }} $

(9)

式中:λ_i是Cov(C)矩阵的特征值; Λ是按照特征值递减顺序排列的对角矩阵; U为特征向量矩阵。

$ \mathit{\boldsymbol{U = }}\left[ {\begin{array}{*{20}{c}} {{u_{11}}}&{{u_{12}}}& \cdots &{{u_{1m}}}\\ {{u_{21}}}&{{u_{22}}}& \cdots &{{u_{2m}}}\\ \vdots&\vdots &{}& \vdots \\ {{u_{m1}}}&{{u_{m2}}}& \cdots &{{u_{mm}}} \end{array}} \right] $

(10)

对二维矩阵C做正交变换可得K-L变换公式:

$ \mathit{\boldsymbol{Y}} = {\mathit{\boldsymbol{U}}^{\rm{T}}}\mathit{\boldsymbol{C}} $

(11)

式中:U^T为U的转置。重构信号公式为:

$ \mathit{\boldsymbol{\hat C}} = \mathit{\boldsymbol{UY}} $

(12)

式中: $\mathit{\boldsymbol{\hat C}}$ 为经K-L滤波后的Contourlet系数矩阵。只要取前p(1≤p≤m)个特征值所对应的特征向量对信号进行重构, 就可以达到衰减随机噪声的目的。

在K-L变换域, Contourlet系数矩阵C的能量可由其协方差矩阵的特征值来度量, 即输入数据的总能量为 $\sum\limits_{i = 1}^m {{\lambda _i}} $ 。使用前p个主分量重构Contourlet系数矩阵 $\mathit{\boldsymbol{\hat C}}$ , 再进行Contourlet反变换则可得到去噪后的地震记录。其中用于K-L反变换重构的特征向量个数p依据能量百分比δ来确定:

$ \delta = \frac{{\sum\limits_{i = 1}^p {{\lambda _i}} }}{{\sum\limits_{i = 1}^m {{\lambda _i}} }} $

(13)

不同尺度的噪声引起的Contourlet系数方差沿着分解层次近似为指数分布, 而同一个尺度内各方向的噪声方差基本相等^[30]。因而不同尺度的噪声方差σ_m可表示为:

$ \sigma _{\rm{m}}^2 = \sigma _1^2{{\rm{e}}^{1 - {i^{1.2}}}} $

(14)

式中:i代表分解层次(尺度); σ₁为最小尺度的噪声方差。

σ₁可采用经典的鲁棒中值估计法求取^[31]:

$ {\sigma _1} = \frac{{{\rm{median}}\left( {\left| {{c_{il}}} \right|} \right)}}{{0.6745}} $

(15)

式中:c_il为第i层、l方向子带系数; median表示求系数中值。

定义第i层、第l个方向子带中大小为m×n, 代表有效信号的Contourlet系数方差为σ, 利用最大似然估计法来求取^[25]:

$ {\sigma ^2} = \max \left[ {\frac{1}{{m \times n}}\sum\limits_{i = 1}^m {\sum\limits_{l = 1}^n {c_{il}^2 - \sigma _{\rm{m}}^2} ,0} } \right] $

(16)

在Contourlet域, 有效信号与噪声的Contourlet系数是相互独立的, 所以输入信号(即含噪信号)的Contourlet系数方差σ_c满足

$ \sigma _{\rm{c}}^2 = {\sigma ^2} + \sigma _{\rm{m}}^2 $

(17)

考虑到实际资料的复杂性, 引入调整因子θ∈(0, 1), 定义K-L变换域能量百分比阈值函数为f(θ), 则:

$ f\left( \theta \right) = 1 - \theta \frac{{\sigma _{\rm{m}}^2}}{{\sigma _{\rm{c}}^2}} $

(18)

当p满足:

$ \delta = \frac{{\sum\limits_{i = 1}^p {{\lambda _i}} }}{{\sum\limits_{i = 1}^m {{\lambda _i}} }} \ge f\left( \theta \right) $

(19)

时, 确定p为用于信号重构的最大特征向量个数。

不同尺度、不同方向子带的Contourlet系数c_il不同, 导致σ_c², σ_m²和λ_i也不同, 公式(18)中的f(θ)和公式(19)中的p受信号驱动而变化, 所以基于Contourlet变换的K-L变换随机噪声衰减方法是一种自适应随机噪声压制方法。

图 5为基于Contourlet变换的K-L变换去噪方法流程, 具体实现步骤如下:

1) 对输入数据进行Contourlet变换, 得到不同尺度、不同方向子带的Contourlet系数矩阵C;

2) 对每个系数矩阵C做K-L变换, 得到对应的特征值和特征向量;

3) 利用(14)~(18)式求取K-L变换域能量百分比阈值函数f(θ);

4) 利用公式(19)求取p值, 确定K-L反变换时的特征向量个数;

5) 利用公式(12)进行K-L反变换, 得到去除随机噪声后的Contourlet系数矩阵 $\mathit{\boldsymbol{\hat C}}$ ;

6) 对去噪后的Contourlet系数矩阵 $\mathit{\boldsymbol{\hat C}}$ 进行Contourlet反变换得到去噪后的地震数据。

2 数值模拟

对小波变换去噪、常规Contourlet变换阈值去噪和本文基于Contourlet变换的K-L变换去噪方法进行了数值模拟对比试验(图 6)。图 6a为数值模拟的原始单炮记录; 图 6b为加入随机噪声的单炮记录, 其信噪比为0.55;图 6c为小波变换阈值去噪结果, 其信噪比为1.02;图 6d为Contourlet变换阈值去噪结果, 其信噪比为1.28;图 6e为本文方法的去噪结果, 其信噪比为2.01, 优于小波变换阈值去噪结果和Contourlet变换阈值去噪结果; 图 6f为本文方法去除的噪声。图 7为本文方法去噪后单炮记录频谱分析结果, 去噪后单炮记录频谱与原始模拟单炮记录频谱相似。说明本文方法不仅对随机噪声有良好的压制作用, 而且对有效信号有很高的保真度。

3 实际资料处理效果分析

采用本文基于Contourlet变换的K-L变换去噪方法对我国南方某工区实际地震资料进行了处理, 如图 8所示。该区单炮记录噪声比较严重(图 8a), 1000~2000ms的有效反射信号被噪声淹没。利用本文方法进行去噪处理后, 单炮记录信噪比得到明显提高, 1000~2000ms有效反射波同相轴清晰可见(图 8b), 去除的噪声中肉眼未见有效反射信号(图 8c), 说明本文方法不仅能够很好地去除地震资料中的随机噪声, 而且不损伤有效信号, 具有良好的保真性。

利用本文方法对整个工区的单炮记录进行叠前去噪, 并分别对叠前去噪处理前、后的单炮记录进行叠加, 结果如图 9所示。

对比可见, 利用本文方法进行叠前去噪处理后, 叠加剖面信噪比得到了显著提高, 有效反射波同相轴连续性增强, 波组特征明显改善。

4 结束语

本文对Contourlet系数进行K-L变换, 利用最大似然估计法和多尺度噪声估计法来估算K-L变换域能量百分比函数, 根据K-L变换域能量百分比函数自适应地确定K-L反变换所用的特征向量, 修改Contourlet系数, 去除随机噪声, 提高地震资料信噪比。数据模拟和实际资料处理结果表明, 本文方法不仅能很好地去除地震资料中的随机噪声, 提高资料信噪比, 而且不损伤有效信号, 具有良好的保真性。