CN102831588B - 一种三维地震图像的降噪处理方法 - Google Patents

Abstract

本发明公开了一种三维地震图像的降噪处理方法，充分利用了混合范数迭代算法的优越性和地震图像连续性因子的准确性，使得地震图像的降噪效果得到了明显的提升。本发明的积极效果是：通过混合范数的引入，使得在地震图像降噪的过程中，何时用1-1范数对图像保边处理，何时用1-2范数对图像平滑处理，有更好的自适应性；通过连续性因子的引入，使得本发明提出的算法在降噪的过程中的自适应性更加明显和准确，在准确的保留原始图像中有用的纹理信息方面，比传统的降噪方法有着更好的效果；从试验结果可以看出，混合范数迭代算法求解总变差模型收敛速度快，计算时间少，且具有很好的鲁棒性。

Description

一种三维地震图像的降噪处理方法

技术领域

本发明涉及一种三维地震图像的降噪处理方法。

背景技术

当今，各国在地震解释技术中都在不断地寻求所采用技术的改进和提高，特别是当地震资料信号信噪比较低时，传统的地震图像降噪技术及处理技术并不能快速准确地给出地质结构和储量信息，所以对地震资料降噪提高信噪比的问题变得越来越急迫。目前，已经展开了通过地震图像处理技术提高地震数据的信噪比及清晰度，以有利于地震解释的工作。而且，经过仔细观察，发现地震图像中的噪声大致可以分为相干噪声和随机噪声。

相干噪声包括折射、面波、多次波、鸣震等，如何处理这些相干噪声，地球物理界一直争论不休。有的学者认为，相干噪声的存在严重影响子波处理、速度分析与动校正，应在叠前加以消除；有的学者则认为，目前还没有一种精确的信噪分离的方法，因而，在消除相干噪声的过程中，有效波一定受到不可逆转的损失，另外叠加是最有效的去噪方法，故在叠前不进行消除相干噪声处理为好。

随机噪声包括环境噪声、地面微震等，它是由各种不可预知的因素综合作用而成，有统一的规律，它在整张记录上随机出现，频带很宽，视速度不确定，无一定的传播方向，比较难以去除。根据随机噪声的特点和产生机制，结合地震勘探的实际情况，可将随机噪声划分为环境噪声、次生噪声、***噪声三大类型。随着地震采集设备的不断更新改进、***噪声对地震勘探的影响已非常微弱，因此，我们主要是针对环境噪声和次生噪声进行分析研究。

目前关于地震信号去噪方法的研究主要集中在以下两个方面：一是研究如何在去噪过程中，不仅能够使水平方向同相轴得到增强，而且能够增强地震记录剖面上倾角较大的倾斜同相轴或弯曲同相轴。由于大多噪声衰减方法都是基于地震记录剖面在水平方向的相关性，这往往使地震记录剖面上的水平同相轴得到加强，而具有较大倾角的倾斜同相轴或弯曲同相轴会被削弱，因此必须研究新的去噪方法，或对现有的去噪方法进行必要的改进，以既能有效地压制噪声又能同时增强水平同相轴和倾斜同相轴及弯曲同相轴的信噪比；其二，是研究对地震记录进行分频去噪的方法。由于地震记录的不同频率成分有不同的信噪比，而不同频率成分的信噪比改进对提高分辨率的作用不同，且地震有效信号的高频成分随深度的增加而损失增大，使其在高频端的信噪比很低。因此，地震信号去噪应该对于不同的频率成分能够区别对待，这就要求所研究方法应具有分频处理的特性。

与本发明相关的现有技术包括：

目前地震图像处理中所采用的降噪方法有很多种，每种方法都有条件，需要根据噪声的类型来选用合适的去噪方法，这样才能取得预期选择合适的方法，对于提高地震图像的信噪比以及分辨率是至关重要的。

由于相干噪声在时间上的出现具有规律性，有明的运动学特征，和有效波存在着频谱差异、视速度差异或到达时间的差异等，因此可用波、方向特性、相干性等进行去除。本发明重点研究的是地震图像中随机噪声的降噪处理。

现阶段常用的随机噪声消除方法有F-X域去噪、多项式拟合、径向预测滤波、矢量分解去噪及中值约束下的矢量分解去噪等。这些方法的设计思路大体分为两类：增强信号或压制噪声，目的是尽可能地提高信噪比。信号增强类在生产实践中经常使用的F-X域反褶积、多项式拟合、径向预测滤波；中值约束下的矢量分解方法则属于压制噪声类。

f-x域滤波是消除随机噪声的一种较好办法。根据线性预测理论和随机噪声无法预测的原理，对叠后剖面上的线性同相轴进行预测，分离信号与噪声，压制随机噪声，并进行有效信号增强。但是该方法的去噪效果受到地震记录信噪比高低的限制，尤其是该方法需要假定同相轴在空间方向上是线性的。所以，当地下构造比较复杂时，使用此方法将产生一系列复杂的问题。

KL变换实质上是一种加权叠加(组合)，因此在剖面上使用时，要求各道信号时间一致，这样导致倾斜界面效果较差，而且，由于该方法的准则是最大叠加能量，显然能量越大的道次的权系数也就越大，假如能量大的原始道主要是噪声，使用此方法后极有可能降低地震信号的信噪比。所以，此方法在实际的应用中通常采用先拉平再变换的方法提取噪声或者有效波。

τ-p变换法是从地震波的视速度角度出发，将地震资料分解到截距-斜率域中，采用切除方式将干扰部分切除，然后利用反变换得到有效信号记录。但由于在实际地震资料中，变换域内的有效信号和干扰往往没有明显的界限，所以该方法在去噪的同时也会损伤一部分有效信号。

时频峰值滤波技术是一项比较新的时频域降噪技术。时频峰值滤波是在时频分析的基础上，利用时频平面的峰值作为瞬时频率估计器，然后将叠加了随机噪声的信号编码为瞬时频率进行估计。这样估计得到的信号消减了随机噪声，使得信噪比得到有效恢复。时频滤波方法避免了对于待估计信号波形的假设，原理十分简单，并且能够增强低信噪比含噪声信号中的带限确定性信号，提高信噪比。

中值滤波方法的本质是中位数法，它是以误差绝对值之和达到最小为判定条件，从而确定滤波器输出的滤波方法，由于它的输出是随着序列的变化而变化的，所以它是一种能够有效压制随机噪声的非线性处理技术。在有效波同相轴水平的情况下，中值滤波能够达到保护同相轴细节特征和最大限度地剔除噪声的目的，但是，当地震剖面上的有效波同相轴非水平时，而且倾斜角度发生剧烈变化时，中值滤波就完全失去了滤波作用，这是中值滤波致命的缺陷。

混沌振子检测技术是利用不同类型信息引起由一个非线性动力学方程所构成混沌振子滤波***相态的不同表现来判别有用信息的存在。此滤波***对强随机噪声具有极强的“免疫力”，即使在低信噪比、频谱重叠、断层结构、低有效覆盖次数等条件下都可以较为准确地确定湮没在强随机噪声下的地震勘探同相轴的时空位置，然后再利用各种普通的去噪方法使湮没在噪声中的同相轴显现出来。

这些方法在生产中得到了广泛的应用，也取得了理想的应用效果。不过，它们的共同前提是有效波具有一定相似性，只有满足一定信噪比的地震记录才能使用这些技术。由于方法的局限性，随机噪声消除技术在生产实践中经常达不到最佳去噪效果。

地震图像的降噪的问题，在过去的几十年里一直被研究着。由于地震图像无可避免地受到噪声的影响，图像降噪过程无论是理论分析或是数值计算都有一定的困难。地震图像降噪处理最基本的任务是在去除由降质***引入的噪声的同时，不丢失原始数据的细节信息，然而抑制噪声和保持细节往往是一对矛盾，也是地震图像处理中至今尚未很好解决的一个问题。

总变差正则化方法正是在这种背景下提出并发展起来的。本发明在经典的Tikhonov正则化方法求解反问题的理论框架下，将传统的总变差正则化方法中的范数进行推广,提出了广义总变差正则化模型，并根据加权迭代最小二乘方法的基本思想，通过加权矩阵将一般的l-p范数转化为l-2范数，从而应用标准的二次优化方法进行迭代求解，这就是本发明重点介绍的算法——混合范数迭代算法。从数值仿真试验中对污染图像的降噪效果来看，与经典降噪算法相比，混合范数迭代（IRMN）算法在计算时间和降噪效果上都具有显著优势。

发明内容

为了克服现有技术的上述缺点，本发明提供了一种三维地震图像的降噪处理方法。

本发明解决其技术问题所采用的技术方案是：一种三维地震图像的降噪处理方法，包括如下步骤：

步骤一、参数初始化，输入三维地震图像数据；

步骤二、建立退化模型，计算出退化矩阵H；

步骤三、计算结构张量S和连续性因子λ；

步骤四、计算正则化变量α；

步骤五、得到目标表达式M(x)，并进行迭代计算；

步骤六、返回步骤四，更新正则化变量α，再次进行运算，直到到达最大迭代次数，完成三维地震图像的降噪处理。

与现有技术相比，本发明的积极效果是：提出了基于正则化混合范数迭代算法的三维地震图像的降噪处理方法，充分利用了混合范数迭代算法的优越性和地震图像连续性因子的准确性，使得地震图像的降噪效果得到了明显的提升。具体优点如下：

（1）混合范数的引入，使得在地震图像降噪的过程中，何时用1-1范数对图像保边处理，何时用1-2范数对图像平滑处理，有更好的自适应性。

（2）连续性因子的引入，使得本发明提出的算法在降噪的过程中的自适应性更加明显和准确，在准确的保留原始图像中有用的纹理信息方面，比传统的降噪方法有着更好的效果。

（3）在经典降噪方法失效的情形下，总变差正则化模型仍然能得到较好的降噪效果。同时从试验结果可以看出，混合范数迭代算法求解总变差模型收敛速度快，计算时间少，且具有很好的鲁棒性。

附图说明

本发明将通过例子并参照附图的方式说明，其中：

图1是本发明方法的流程示意图。

具体实施方式

一种三维地震图像的降噪处理方法，如图1所示，包括如下步骤：

步骤一、参数初始化（包括迭代次数的设定，Iter=0），输入三维地震图像数据；

步骤二、建立退化模型，计算出退化矩阵H：

地震图像的质量下降叫做退化，退化的形式有模糊、失真、有噪声等。图像的退化过程可以看成是噪声的污染过程，退化模型如下：

根据以上退化模型，得到退化后的图像为：

g(x,y,z)=H·f(x,y,z)+n(x,y,z)（2-1）

其中，f(x,y,z)是原始图像，g(x,y,z)是退化后的图像，H[]为综合所有退化因素的函数，n(x,y,z)是噪声。

我们对H有如下假设：

①H是线性的。

H[k₁f₁(x,y,z)+k₂f₂(x,y,z)]＝k₁Hf(x,y,z)+k₂Hf₂(x,y,z)(2-2）

②H是空间（或移位）不变的。对任一个f(x,y)和任意常数α和β都有：

H·f(x-α,y-β，z-γ)＝g(x-α,y-β，z-γ)（2-3）

也就是说图像上任一点的运算结果只取决于该点的输入值，而与坐标位置无关。

三维图像降噪的过程可以看成是一个离散卷积的过程，计算公式如下：

$g(x,y,z)=\underset{i=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{k=0}{\overset{L-1}{\mathrm{\Σ}}}f(x,y,z)h(x-i,y-j,z-k)+n(x,y,z)---(2-4)$

用矩阵表示次计算过程如下：

$g=\mathrm{Hf}+n=\left[\begin{array}{cccc}{H}_0& H_{M-1}& ...& H_1\\ H_1& H_0& ...& H_2\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ H_{M-1}& H_{M-2}& ...& H_0\end{array}\right]\left[\begin{array}{c}{f}_e\left(\mathrm{0,0}\right)\\ f_e\left(\mathrm{0,1}\right)\\ .\\ .\\ .\\ f_e(M-1,N-1)\end{array}\right]+\left[\begin{array}{c}{n}_e\left(\mathrm{0,0}\right)\\ n_e\left(\mathrm{0,1}\right)\\ .\\ .\\ .\\ n_e(M-1,N-1)\end{array}\right]---(2-5)$

其中，对于一个M×N×L像素的三维地震图像来说，H的大小是MNL×MNL，它是一个分块轮换矩阵，而H_i是由扩展矩阵h(x,y)的第i行循环构成：

$H_i=\left[\begin{array}{cccc}{h}_e(i,0)& h_e(i,N-1)& ...& h_e(i,1)\\ h_e(i,1)& h_e(i,0)& ...& h_e(i,2)\\ .& .& ...& .\\ .& .& ...& .\\ .& .& ...& .\\ h_e(i,N-1)& h_e(i,N-2)& ...& h_e(i,0)\end{array}\right]---(2-6)$

可以看出H_i也是一个右移轮换矩阵，根据搜索（Hunt）的结果，它的这种性质便于矩阵H的对角化，即H=WDW^-1，其中W是由MN×MN个L×L的块组成，而D是一个由H矩阵的特征值组成的对角矩阵，再对式(2-5)进行二维傅里叶变换得到式(2-7):

g=WDW^-1f+n

G(u,v)=H(u，v)F(u,v)+N(u,v)(2-7)

这样整个求解过程就得到了简化，因为***是线性时不变的，在空间域中建立的退化模型可通过分块轮换矩阵对角化导出频域中的恢复滤波器，将庞大的空域运算转化为较少的频域运算，从而很容易地计算出退化矩阵H。

步骤三、计算结构张量S和连续性因子λ：

连续性因子λ的引入，主要是控制何时使用最小四阶矩范数（LMF），何时使用最小均方值范数(LMS）,或者两者同时混合使用所占的比重。下面介绍连续性因子的计算方法和步骤：

首先，计算三维地震数据的结构张量S：

$S=\left(\begin{array}{ccc}\frac{\∂u\∂u}{\∂x\∂x}& \frac{\∂u\∂u}{\∂x\∂y}& \frac{\∂u\∂u}{\∂x\∂z}\\ \frac{\∂u\∂u}{\∂x\∂y}& \frac{\∂u\∂u}{\∂y\∂y}& \frac{\∂u\∂u}{\∂y\∂z}\\ \frac{\∂u\∂u}{\∂x\∂z}& \frac{\∂u\∂u}{\∂y\∂z}& \frac{\∂u\∂u}{\∂z\∂z}\end{array}\right)---(2-8)$

其中，u(x,y,z)为输入的初始地震图像，表示x,y,z三个方向的偏导数。

然后，可根据下式计算连续性因子：

$\mathrm{\λ}=\frac{\mathrm{diag}\left(S^{0}S\right)}{\mathrm{diag}\left(S^0\right)\mathrm{diag}\left(S\right)}---(2-9)$

其中，S⁰为初始化结构张量，S为当前迭代时间的结构张量，连续性因子λ值域为[0,1]，图像平滑均质区域趋近于0，图像边缘区域趋近于1。

步骤四、计算正则化变量α：

正则化变量主要控制混合范数函数和平滑滤波函数之间的关系，它影响图像数据的可靠性和平滑性，它可以通过式(2-10)计算得到：

$\mathrm{\α}\left(x\right)=\frac{\mathrm{\τJ}\left(x\right)}{1-\mathrm{\τ}{\left|\right|\mathrm{Cx}\left|\right|}_2^2}=\frac{[1-\mathrm{\λ}]{\left|\right|y-\mathrm{Hx}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|y-\mathrm{Hx}\left|\right|}_4^4}{1/\mathrm{\τ}-{\left|\right|\mathrm{Cx}\left|\right|}_2^2}---(2-10)$

其中，J为混合范数函数，λ为连续性因子，τ表示凸性因子，C是一个高通滤波器，当α(x)是收敛的。

步骤五、得到目标表达式M(x)，并进行迭代计算：

综合以上变量得到目标表达式如下：

$M\left(x\right)=[1-\mathrm{\λ}]{\left|\right|y-\mathrm{Hz}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|y-\mathrm{Hx}\left|\right|}_4^4+\mathrm{\α}\left(x\right){\left|\right|\mathrm{Cx}\left|\right|}_2^2$

$=J\left(x\right)+\mathrm{\α}\left(x\right){\left|\right|\mathrm{Cx}\left|\right|}_2^2---(2-11)$

其中，J(x)为最小均方混合范数函数（LMMN），α(x)为正则化变量，C是一个高通滤波器，在本发明中，用3-D的拉普拉斯算子来表示这个高通滤波器。下面给出了2-D拉普拉斯算子表达式：

${\mathrm{\δ}}_{i,j}\left(\begin{array}{ccc}0& -0.25& 0\\ -0.25& 1& -0.25\\ 0& -0.25& 0\end{array}\right)---(2-12)$

其中，δ_i，j=∑_m∑_nh_i，j(m,n)，h_i，j是PSF散射函数。

由以上步骤我们得到了对地震图像降噪的目标表达式，然后本发明中通过运用连续逼近的迭代算法来最小化式(2-11)的误差问题。M(x)在x方向的梯度向量可表示为：

$\▿M\left(x\right)=-2[1-\mathrm{\λ}]H^T(y-\mathrm{Hx})-4\mathrm{\λ}{H}^T{(y-\mathrm{Hx})}^3+2\mathrm{\α}\left(x\right)C^T\mathrm{Cx}+{\▿}_x\mathrm{\α}\left(x\right){\left|\right|\mathrm{Cx}\left|\right|}_2^2---(2-13)$

其中，式(2-13)满足当时，正因为如此，连续逼近迭代算法可以用下面的式子表达：

x_k＋1＝x_k+β{[1-λ]H^T(y-Hx)+2λH^T(y-Hx)³-α(x)C^TCx}(2-14)

其中，β是松弛变量，它用来控制迭代过程中计算式的收敛性。在每一次的迭代时，α(x)都由新的α(x_k)替换，每次迭代后，梯度向量也由新计算出的替换。通过式(2-14)的迭代算法，完成对三维地震图像的降噪处理。

步骤六、返回步骤四，更新参数α和梯度向量，再次进行运算，Iter=Iter+1，直到到达最大迭代次数IterMax，完成三维地震图像的降噪处理。

Claims (2)

1.一种三维地震图像的降噪处理方法，其特征在于：包括如下步骤：

步骤一、参数初始化，包括迭代次数的设定，输入三维地震图像数据；

步骤二、建立退化模型，g(x,y,z)＝H·f(x,y,z)+n(x,y,z)，

其中，f(x,y,z)是原始图像，g(x,y,z)是退化后的图像，H[]可以理解为综合所有退化因素的函数，n(x,y,z)是噪声；

计算出退化矩阵H；

步骤三、计算结构张量S和连续性因子λ；

根据下式计算连续性因子：

$\mathrm{\λ}=\frac{\mathrm{diag}\left(S^{0}S\right)}{\mathrm{diag}\left(S^0\right)\mathrm{diag}\left(S\right)}---(2-9)$

其中，S⁰为初始化结构张量，S为当前迭代时间的结构张量，

所述S为当前迭代时间的结构张量，根据下式计算得到：

$S=\left(\begin{array}{ccc}\frac{\∂u\∂u}{\∂x\∂x}& \frac{\∂\∂}{\∂x\∂y}& \frac{\∂u\∂u}{\∂x\∂z}\\ \frac{\∂u\∂u}{\∂x\∂y}& \frac{\∂u\∂u}{\∂y\∂y}& \frac{\∂u\∂u}{\∂y\∂z}\\ \frac{\∂u\∂u}{\∂x\∂z}& \frac{\∂u\∂u}{\∂y\∂z}& \frac{\∂u\∂u}{\∂z\∂z}\end{array}\right)---(2-8)$

其中，u(x,y,z)为输入的初始地震图像，表示x,y,z三个方向的偏导数；

步骤四、计算正则化变量α；

步骤五、得到目标表达式M(x)，并进行迭代计算；

步骤六、返回步骤四，更新正则化变量α，再次进行运算，直到到达最大迭代次数，完成三维地震图像的降噪处理。

2.根据权利要求1所述的一种三维地震图像的降噪处理方法，其特征在于：所述正则化变量α的计算公式为：

$\mathrm{\α}\left(x\right)=\frac{\mathrm{\τJ}\left(x\right)}{1-\mathrm{\τ}{\left|\right|\mathrm{Cx}\left|\right|}_2^2}=\frac{[1-]\mathrm{\λ}{\left|\right|y-\mathrm{Hx}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|y-\mathrm{Hx}\left|\right|}_4^4}{1/\mathrm{\τ}-{\left|\right|\mathrm{Cx}\left|\right|}_2^2}$

其中，J为混合范数函数，λ为连续性因子，τ表示凸性因子，C是一个高通滤波器。