CN105676292A - 一种基于二维曲波变换的三维地震数据去噪方法 - Google Patents

Abstract

本发明提出一种基于二维曲波变换的三维地震数据去噪方法，其特征在于，首先抽取含噪三维地震数据的时间切片，对其进行多尺度多方向二维曲波变换得到曲波域系数，然后在曲波域采用局部阈值法去噪方法，对曲波变换后的每一个尺度都选取一个软阈值算子，通过阈值处理，得到各尺度下的有效波曲波系数，最后将提取出来的有效波曲波系数进行逆变换重构出地震信号，从而达到去噪目的。本发明提出了采用二维曲波变换进行三维地震数据去噪，并且提出局部阈值法去噪方法，采用软阈值算子，对每个时间切片进行单独去噪，克服了传统全局阈值法去噪效果不理想的缺点，同时降低了对计算机内存的要求，大幅度提高了计算效率，节省了运算时间。

Description

一种基于二维曲波变换的三维地震数据去噪方法

技术领域

本发明涉及的是地震数据去噪的方法，具体是一种基于二维曲波变换的三维地震数据去噪方法。

技术背景

随着我国地震勘探进程的不断向前发展，山区，沙漠，厚黄土、砂土等复杂地区的地震勘探项目逐渐增多，而这些地区表层激发条件都不够理想，再加上野外各种外界环境的干扰，所采集到的地震数据包含各种严重的噪声干扰，掩盖了有效波信息，使得有效信号同相轴模糊不清，能量相对较弱，信噪比低。尽管在野外可以采取相应的一些抗噪措施来提高信噪比。但是，在多风季节或受其它复杂条件限制时，仅靠野外采集阶段完全压制随机噪声是不可能，这就需要在室内进行有效地去噪工作，提高叠前地震资料的信噪比，以便后续资料的处理。

目前，基于数学变换的去噪方法较多，这些方法充分利用地震数据在数学变换域内稀疏特点，通过对变换域中的系数进行阈值处理，从而达到去噪的目的。而为了达到所需要的去噪效果，这就要求稀疏表示所采用的基函数能够捕捉到地震波前，只要保留少数较大的稀疏系数就可以表示出原始数据的主要特征，而滤除的大量小系数不影响数据的主要特征。

假设采用多尺度多方向的二维曲波变换进行去噪，由于二维曲波变换由各向异性的曲线状基元所构成，可以更加稀疏地表示地震波前特征，从而能够弥补了其他数学变换方法的不足。然而在去噪领域中，常规方法常被用于压制二维地震数据的随机噪声，但随机噪声存在于三维地震数据空间中，导致以往二维去噪方法无法获得满意的去噪结果。而对于曲波变换如果想真正压制三维地震数据随机噪声，理论上讲应该需要使用三维曲波变换，但三维曲波变换运算速度较慢，处理时间过长，远满足不了海量数据处理的要求。

因此，本发明进行折中考虑，将二维曲波变换应用于三维地震资料去噪处理，提出在每次去噪过程中逐次对时间切片进行多尺度多方向二维曲波正变换，从而得到不同尺度层的曲波系数，同时根据有效信号系数和噪声系数的分布不一样，改变以往采用单一全局阈值会损伤部分有效波的做法，提出采用局部阈值方法，也即每一尺度都采用一个阈值参数，这样可以分别提取出每一尺度的有效波系数，最后再进行二维曲波反变换，从而完成整个处理流程。

发明内容

本发明的目的是为了能够快速高精度去除地震勘探数据中的噪声干扰，而提供了一种基于二维曲波变换的三维地震数据去噪方法。

本发明提出一种基于二维曲波变换的三维地震数据去噪方法，首先抽取含噪三维地震数据的时间切片，对其进行多尺度多方向二维曲波变换得到曲波域系数，然后在曲波域采用局部阈值法去噪方法，对二维曲波变换后的每一个尺度都选取一个软阈值算子，通过阈值处理，得到各尺度下的有效波曲波系数，最后将提取出来的有效波曲波系数进行逆变换重构出地震信号，从而达到去噪目的。

一个含噪声的三维地震信号的模型可以表示成如下的形式：

f(i,j,t)＝s(i,j,t)+k·e(i,j,t)i＝1,2,...,m,j＝1,2,..·,n,t＝1,2,...,k

式中f(i,j,t)为三维地震数据含噪信号，s(i,j,t)为三维不含噪地震信号，e(i,j,t)为噪声信号，k表示噪声水平值。去噪的过程就是从含噪信号f(i,j,t)中，提取真实信号s(i,j,t)，去除噪声干扰信号e(i,j,t)。

在数据去噪过程中，采用局部阈值法，具体去噪步骤如下：

(1)首先逐次抽取三维含噪地震数据的时间切片，然后选择分解尺度记为N，将含噪时间切片进行N尺度二维曲波变换，得到含噪的曲波系数，这些系数中少量的较大曲波系数能够代表信号本身，而大部分较小值的曲波系数则表示高频噪声干扰信号。

(2)根据各分解尺度上曲波系数的分布情况，选择与分解尺度有关的局部阈值参数，以反映曲波系数在不同尺度上的不同特征，然后对该尺度曲波系数分量进行软阈值处理，保留有效波曲波系数。

(3)根据经过软阈值处理后的各尺度有效波曲波系数分量，进行地震信号的二维曲波反变换，对反变换后的时间切片进行组合，所得到的数据体即为本发明去噪后的三维地震数据。

进一步，所述二维曲波变换的定义为：

$C(j,l,k)=<f,{\mathrm{\&phi;}}_{j,l,k}>={\&Integral;}_{R^2}f\left(x\right)\stackrel{\&OverBar;}{{\mathrm{\&phi;}}_{j,l,k}\left(x\right)}dx$

式中：φ_j,l,k表示曲波函数，j,l,k分别表示尺度，方向和位置参数，f(x)为

地震数据，其频率域定义式为：

$C\left(j,l,k\right)=\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}\&Integral;\hat{f}\left(\mathrm{\&omega;}\right)\stackrel{\&OverBar;}{{\widehat{\mathrm{\&phi;}}}_{j,l,k}\left(\mathrm{\&omega;}\right)}d\mathrm{\&omega;}=\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}\&Integral;\hat{f}\left(\mathrm{\&omega;}\right)U_j\left(R_{{\mathrm{\&theta;}}_l}\mathrm{\&omega;}\right)e^{i<x_k^{\left(j,l\right)},\mathrm{\&omega;}>}d\mathrm{\&omega;}$

经过变换后得到的曲波系数，用C{j}{l}(k₁,k₂)表示其结构，其中j表示尺度，l表示方向，(k₁,k₂)表示j尺度l方向上的矩阵系数。

进一步，局部阈值参数，其表达式如下：

λ_c＝aστ

式中λ_c表示局部阈值参数，a表示与尺度方向有关的经验参数，σ表示最高尺度层曲波系数标准差的估计值，τ代表曲波域内所有曲波系数的标准差估计值。

进一步，所述软阈值算子，其表达式如下：

$F=\left\{\begin{array}{cc}F-{\mathrm{\&lambda;}}_c,& F\&GreaterEqual;{\mathrm{\&lambda;}}_c\\ F+{\mathrm{\&lambda;}}_c,& F\&le;-{\mathrm{\&lambda;}}_c\\ 0,& F<{\mathrm{\&lambda;}}_c\end{array}\right.$

式中F表示软阈值算子，λ_c表示与分解尺度有关的局部阈值参数，不同分解尺度的阈值参数不一样。

本发明的优点:本发明采用了具有多尺度和多方向性的二维曲波变换进行三维地震资料去噪，通过逐次对时间切片进行去噪处理，实现了基于二维曲波变换的三维地震资料去噪处理方法，从而避免了计算速度慢，效率低的缺点，大幅度地提高了计算效率，节省了运算时间。同时在去噪过程中，提出局部阈值参数的去噪方法，尽可能地保护了微弱的有效波信号，从而使反射波同相轴更加连续、清晰，提高了地震数据的信噪比，降低了对计算机内存的要求。

附图说明

图1是本发明实施例中三维地震数据去噪流程图。

图2是原始地震数据及其加噪对比图。

图3是为含噪数据二维曲波6尺度分解图。

图4是局部阈值去噪结果图。

图5是去除的噪声剖面图。

具体实施方式

以下实施案例用于说明本发明，但不用来限制本发明的范围。

实施例1

实现该方法的步骤主要包括，去噪方程的构建，二维曲波变换，局部阈值法，阈值算子处理等。具体步骤如下：

步骤1：去噪方程的构建。一个含噪声的三维地震信号的模型可以表示成如下的形式：

f(i,j,t)＝s(i,j,t)+k·e(i,j,t)i＝1,2,...,m,j＝1,2,...,n,t＝1,2,...,k

式中f(i,j,t)为三维地震数据含噪信号，s(i,j,t)为三维不含噪地震信号，e(i,j,t)为噪声信号，k表示噪声水平值。去噪的过程就是从含噪信号f(i,j,t)中，提取真实信号s(i,j,t)，去除噪声干扰信号e(i,j,t)。

步骤2：首先逐次抽取三维含噪地震数据的时间切片，然后选择合适的分解尺度(记为N)，将含噪时间切片进行N尺度二维曲波变换，得到含噪的曲波系数，这些系数中少量的较大曲波系数能够代表信号本身，而大部分较小值的曲波系数则表示高频噪声干扰信号。

所述二维曲波变换的定义为：

$C(j,l,k)=<f,{\mathrm{\&phi;}}_{j,l,k}>={\&Integral;}_{R^2}f\left(x\right)\stackrel{\&OverBar;}{{\mathrm{\&phi;}}_{j,l,k}\left(x\right)}dx$

式中：φ_j,l,k表示曲波函数，j,l,k分别表示尺度，方向和位置参数，f(x)为地震数据，其频率域定义式为：

$C\left(j,l,k\right)=\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}\&Integral;\hat{f}\left(\mathrm{\&omega;}\right)\stackrel{\&OverBar;}{{\widehat{\mathrm{\&phi;}}}_{j,l,k}\left(\mathrm{\&omega;}\right)}d\mathrm{\&omega;}=\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}\&Integral;\hat{f}\left(\mathrm{\&omega;}\right)U_j\left(R_{{\mathrm{\&theta;}}_l}\mathrm{\&omega;}\right)e^{i<x_k^{\left(j,l\right)},\mathrm{\&omega;}>}d\mathrm{\&omega;}$

经过变换后得到的曲波系数，可用C{j}{l}(k₁,k₂)表示其结构，其中j表示尺度，l表示方向，(k₁,k₂)表示j尺度l方向上的矩阵系数。

步骤3：根据各分解尺度上曲波系数的分布情况，选择与分解尺度有关的局部阈值参数，以反映曲波系数在不同尺度上的不同特征，然后对该尺度曲波系数分量进行软阈值量化处理，保留有效波曲波系数。

去噪后的有效波信号可以下述方法估算出，即

s＝C^-1(F(CS))

式中C表示二维曲波变换，C^-1表示逆二维曲波变换，F表示软阈值算子，其表达式如下：

$F=\left\{\begin{array}{cc}F-{\mathrm{\&lambda;}}_c,& F\&GreaterEqual;{\mathrm{\&lambda;}}_c\\ F+{\mathrm{\&lambda;}}_c,& F\&le;-{\mathrm{\&lambda;}}_c\\ 0,& F<{\mathrm{\&lambda;}}_c\end{array}\right.$

式中λ_c表示与分解尺度有关的阈值参数，不同分解尺度的阈值参数不一样，其表达式如下：

λ_c＝aστ

式中λ_c表示局部阈值参数，a表示与尺度方向有关的经验参数，σ表示最高尺度层曲波系数标准差的估计值，τ代表曲波域内所有曲波系数的标准差估计值。

步骤4：根据经过软阈值量化处理后的各尺度有效波曲波系数分量，进行地震信号的二维曲波反变换，对反变换后的时间切片进行组合，所得到的数据体即为本发明去噪后的三维地震数据。

实现该方法具体操作为：

为了详细比较理论模型中基于二维曲波变换的三维地震资料去噪效果，本发明定义信噪比公式为SNR＝20log₁₀||x₀||₂/||x-x₀||₂，x₀表示原始不含噪声的模型(原始数据)，x表示去除噪声后的地震数据，单位为dB，信噪比越高，表示去噪效果越好。同时去噪过程中二维曲波变换的尺度数为6，最粗尺度上的角度数为8。

本发明采用声波有限差分方法，模拟出无噪声干扰的地震正演剖面，并对所得到的正演地震数据按检波器，炮点以及时间进行排列成三维数据体，其中炮距和道距都为12米，采样率4毫秒。由于不能全部显示理论的三维数据模型，本发明只能从三个不同的方向进行显示，理想原始数据如图2(a)所示(图2(a)表示原始模型数据)，其中时间切片为0.44s，共炮点对应距离为1524m(128炮)，共检波点对应的距离为1524m(128道)。然后对其加入一定的随机噪声，如图2b所示(图2(b)表示三维加噪地震数据)。由于地震数据结构复杂，阈值参数的选取也与分解的尺度和方向有关系，单一的全局阈值参数去噪结果会导致一些细节特征没有很好地保持，损失部分有效波系数。图3为图2(b)含噪时间切片(0.44s)所分解的6尺度系数图(图3(c)表示第1尺度曲波系数，图3(d)表示第2尺度曲波系数，图3(e)表示第3尺度曲波系数，图3(f)表示第4尺度曲波系数，图3(g)表示第5尺度曲波系数，图3(h)表示第6尺度曲波系数)，该尺度系数为其它尺度曲波系数置零而进行曲波逆变换而得来的。可以看出尺度1～尺度6含有不同程度的噪声，不能采用单一的全局阈值参数进行处理，因此，本发明选用与尺度有关的局部阈值参数，以反映地震数据和曲波系数在不同尺度上的不同特征。从图3可以看出，尺度1和尺度2几乎不出现随机噪声，可以不进行阈值处理，而尺度3～尺度6含有不同程度的噪声，从而采取的阈值参数应该也不同，根据各尺度曲波系数的阈值参数测试，第3尺度到第6尺度分别采用的阈值为近似保留本尺度40％、20％、10％、和5％的最大曲波系数，然后将保留下来的曲波系数进行曲波逆变换，从而得到去噪结果如图4所示，信噪比为15.63dB，可以看出该地震数据的信噪比得到了提高，有效波同相轴比较连续，从图5中也可以看出本发明方法去除噪声比较彻底，基本上不损失有效波信号，并且由于采用二维曲波变换进行处理，相对于三维曲波变换方法来讲，大幅度地缩短了计算时间。

Claims (7)

1.一种基于二维曲波变换的三维地震数据去噪方法，其特征在于，首先抽取含噪三维地震数据的时间切片，对其进行多尺度多方向二维曲波变换得到曲波域系数，然后在曲波域采用局部阈值法去噪方法，对二维曲波变换后的每一个尺度都选取一个软阈值算子，通过阈值处理，得到各尺度下的有效波曲波系数，最后将提取出来的有效波曲波系数进行逆变换重构出地震信号，从而达到去噪目的。

2.根据权利要求1所述的一种基于二维曲波变换的三维地震数据去噪方法，其特征在于，一个含噪声的三维地震信号的模型表示成如下的形式：

f(i,j,t)＝s(i,j,t)+k·e(i,j,t)i＝1,2,…,m,j＝1,2,···,n,t＝1,2,…,k

式中f(i,j,t)为三维地震数据含噪信号，s(i,j,t)为三维不含噪地震信号，e(i,j,t)为噪声信号，k表示噪声水平值。去噪的过程就是从含噪信号f(i,j,t)中，提取真实信号s(i,j,t)，去除噪声干扰信号e(i,j,t)。

3.根据权利要求1或2所述的一种基于二维曲波变换的三维地震数据去噪方法，其特征在于，在数据去噪过程中，采用局部阈值法，具体去噪步骤如下：

(1)首先逐次抽取三维含噪地震数据的时间切片，然后选择分解尺度记为N，将含噪时间切片进行N尺度二维曲波变换，得到含噪的曲波系数，这些系数中少量的较大曲波系数能够代表信号本身，而大部分较小值的曲波系数则表示高频噪声干扰信号；

(2)根据各分解尺度上曲波系数的分布情况，选择与分解尺度有关的局部阈值参数，以反映曲波系数在不同尺度上的不同特征，然后对该尺度曲波系数分量进行软阈值处理，保留有效波曲波系数；

(3)根据经过软阈值处理后的各尺度有效波曲波系数分量，进行地震信号的二维曲波反变换，对反变换后的时间切片进行组合，所得到的数据体即为去噪后的三维地震数据。

4.根据权利要求3所述的一种基于二维曲波变换的三维地震数据去噪方法，其特征在于，所述二维曲波变换的定义为：

$C(j,l,k)=<f,{\mathrm{\&phi;}}_{j,l,k}>={\&Integral;}_{R^2}f\left(x\right)\stackrel{\&OverBar;}{{\mathrm{\&phi;}}_{j,l,k}\left(x\right)}dx$

式中：φ_j,l,k表示曲波函数，j,l,k分别表示尺度，方向和位置参数，f(x)为地震数据，其频率域定义式为：

$C\left(j,l,k\right)=\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}\&Integral;\hat{f}\left(\mathrm{\&omega;}\right)\stackrel{\&OverBar;}{{\widehat{\mathrm{\&phi;}}}_{j,l,k}\left(\mathrm{\&omega;}\right)}d\mathrm{\&omega;}=\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}\&Integral;\hat{f}\left(\mathrm{\&omega;}\right)U_j\left(R_{{\mathrm{\&theta;}}_l}\mathrm{\&omega;}\right)e^{i<x_k^{\left(j,l\right)},\mathrm{\&omega;}>}d\mathrm{\&omega;}$

经过变换后得到的曲波系数，用C{j}{l}(k₁,k₂)表示其结构，其中j表示尺度，l表示方向，(k₁,k₂)表示j尺度l方向上的矩阵系数。

5.根据权利要求3所述一种基于二维曲波变换的三维地震数据去噪方法，其特征在于，所述局部阈值参数，其表达式如下：

λ_c＝αστ

式中λ_c表示局部阈值参数，a表示与尺度方向有关的经验参数，σ表示最高尺度层曲波系数标准差的估计值，τ代表曲波域内所有曲波系数的标准差估计值。

6.根据权利要求3所述一种基于二维曲波变换的三维地震数据去噪方法，其特征在于，所述软阈值算子，其表达式如下：

$F=\left\{\begin{array}{cc}F-{\mathrm{\&lambda;}}_c,& F\&GreaterEqual;{\mathrm{\&lambda;}}_c\\ F+{\mathrm{\&lambda;}}_c,& F\&le;-{\mathrm{\&lambda;}}_c\\ 0,& F<{\mathrm{\&lambda;}}_c\end{array}\right.$

式中F表示软阈值算子，λ_c表示与分解尺度有关的局部阈值参数，不同分解尺度的阈值参数不一样。

7.根据权利要求1所述一种基于二维曲波变换的三维地震数据去噪方法，其特征在于，该方法的步骤包括，去噪方程的构建，二维曲波变换，局部阈值法，阈值算子处理，具体步骤如下：

步骤1：去噪方程的构建；一个含噪声的三维地震信号模型表示成如下的形式：

f(i,j,t)＝s(i,j,t)+k·e(i,j,t)i＝1,2,…,m,j＝1,2,···,n,t＝1,2,…,k

式中f(i,j,t)为三维地震数据含噪信号，s(i,j,t)为三维不含噪地震信号，e(i,j,t)为噪声信号，k表示噪声水平值；去噪的过程就是从含噪信号f(i,j,t)中，提取真实信号s(i,j,t)，去除噪声干扰信号e(i,j,t)；

步骤2：首先逐次抽取三维含噪地震数据的时间切片，然后选择分解尺度记为N，将含噪时间切片进行N尺度二维曲波变换，得到含噪的曲波系数，这些系数中少量的较大曲波系数能够代表信号本身，而大部分较小值的曲波系数则表示高频噪声干扰信号；

所述二维曲波变换的定义为：

$C(j,l,k)=<f,{\mathrm{\&phi;}}_{j,l,k}>={\&Integral;}_{R^2}f\left(x\right)\stackrel{\&OverBar;}{{\mathrm{\&phi;}}_{j,l,k}\left(x\right)}dx$

式中：φ_j,l,k表示曲波函数，j,l,k分别表示尺度，方向和位置参数，f(x)为地震数据，其频率域定义式为：

$C\left(j,l,k\right)=\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}\&Integral;\hat{f}\left(\mathrm{\&omega;}\right)\stackrel{\&OverBar;}{{\widehat{\mathrm{\&phi;}}}_{j,l,k}\left(\mathrm{\&omega;}\right)}d\mathrm{\&omega;}=\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}\&Integral;\hat{f}\left(\mathrm{\&omega;}\right)U_j\left(R_{{\mathrm{\&theta;}}_l}\mathrm{\&omega;}\right)e^{i<x_k^{\left(j,l\right)},\mathrm{\&omega;}>}d\mathrm{\&omega;}$

经过变换后得到的曲波系数，用C{j}{l}(k₁,k₂)表示其结构，其中j表示尺度，l表示方向，(k₁,k₂)表示j尺度l方向上的矩阵系数；

步骤3：根据各分解尺度上曲波系数的分布情况，选择与分解尺度有关的局部阈值参数，以反映曲波系数在不同尺度上的不同特征，然后对该尺度曲波系数分量进行软阈值处理，保留有效波曲波系数；

去噪后的有效波信号用下述方法估算出，即

s＝C^-1(F(CS))

式中C表示二维曲波变换，C^-1表示逆二维曲波变换，F表示软阈值算子，其表达式如下：

$F=\left\{\begin{array}{cc}F-{\mathrm{\&lambda;}}_c,& F\&GreaterEqual;{\mathrm{\&lambda;}}_c\\ F+{\mathrm{\&lambda;}}_c,& F\&le;-{\mathrm{\&lambda;}}_c\\ 0,& F<{\mathrm{\&lambda;}}_c\end{array}\right.$

式中λ_c表示与分解尺度有关的阈值参数，不同分解尺度的阈值参数不一样，其表达式如下：

λ_c＝aστ

式中λ_c表示局部阈值参数，a表示与尺度方向有关的经验参数，σ表示最高尺度层曲波系数标准差的估计值，τ代表曲波域内所有曲波系数的标准差估计值；

步骤4：根据经过软阈值处理后的各尺度有效波曲波系数分量，进行地震信号的二维曲波反变换，对反变换后的时间切片进行组合，所得到的数据体即为去噪后的三维地震数据。