CN104730576A - 基于Curvelet变换的地震信号去噪方法 - Google Patents

Abstract

本发明涉及一种基于Curvelet变换的地震信号去噪方法，建立一个不同子波长度的楔形模型，模型中只存在一个薄层，不同子波分辨薄层和楔形的能力也不同，选用时间分辨率低的地震记录作为匹配道，时间分辨率高的地震记录作为参考道，通过Curvelet变换的匹配滤波方法对匹配道地震记录以及掺加高斯白噪声的地震记录进行匹配处理，并与参考道分辨率高的地震记录进行对比。经实际对比试验结果表明，通过FITSA算法同时处理两方面的影响，匹配后的结果波形一致性高，并且噪声得到了好的压制，明显的提高了地震数据处理的分辨率、信噪比。在改善信噪比的同时极大地提高了地震资料的处理效率。

Description

基于Curvelet变换的地震信号去噪方法

技术领域

本发明涉及一种地震数据的处理方法，尤其是基于Curvelet变换的地震信号去噪方法。

技术背景

地震勘探中，不同震源的地震数据受匹配滤波方法的限制，匹配滤波后的地震数据仍存在信噪比低、分辨率低等一系列差异；或新老地震资料进行匹配处理后出现分层不清晰，使得波阻抗反演的精度难以保证和地震属性参数的提取困难；或不同区块地震资料的拼接也存在同相轴不连续、信噪比低等，给地震数据匹配处理带来了困难，这一系列问题都严重制约着地震数据处理的分辨率、信噪比。

吉林大学2008年硕士论文公开了《基于Curvelet变换的地震数据去噪方法研究》，介绍了一种基于Curvelet变换的地震数据处理技术研究；将第二代Curvelet变换引入地震数据处理,通过选取恰当的阈值,对地震数据进行Curvelet变换。计算了模拟数据的Curvelet变换,对比了不同信噪比下的去噪效果。并且通过Curvelet变换对实际地震数据进行去噪处理，实验结果得到了比较好的效果。

CN103543469A公开了《一种基于小波变换的小尺度阈值去噪方法》，首先，小尺度地对地震数据进行扫描，得到一个小时窗内的相关系数值，其次，设置一个阈值对这个小尺度时窗内的数据进行判断，是以地震信号为主还是以噪声信号为主，然后采用合适的分频数及合适的小波阈值进行去噪，最后将去噪后的小波尺度进行小波重构，从而得到去噪后信噪比较高的地震道集。

CN103399348A公开了一种《基于Shearlet变换的地震信号去噪方法》，读入二维地震剖面数据，将二维地震数据S扩展为一长宽为奇数的方阵S1，构造频域方向滤波器组，将各个变换矩阵分别与S1信号向量相乘，再分别做二维傅里叶反变换,得到各个方向和尺度的Shearlet系数C_i,j，阈值处理，将经过阈值处理的Shearlet变换系数做Shearlet反变换获得去噪之后的地震数据。本发明将含有噪声的地震信号做拉普拉斯分解然后利用剪切波函数进行滤波处理，得到对应剪切波系数，通过阈值处理滤除噪声信号，再通过逆非下采样剪切波变换恢复去噪的信号，获得较好的去噪效果，具有较好的实用价值。

上述现有技术相比虽然在一定程度上能够有效的压制噪声，但过程复杂，计算速度较慢，信噪比低，地震数据处理效果不是很好。

发明内容

本发明的目的就是针对上述现有技术的不足，提出了将采集到的地震数据进行匹配滤波与噪声压制，通过FITSA算法同时处理两方面的影响，在改善信噪比的同时极大地提高地震资料处理效率。

本发明的目的是通过以下方式实现的：

建立一个不同子波长度的楔形模型，模型中只存在一个薄层，不同子波分辨薄层和楔形的能力也不同，选用时间分辨率低的地震记录作为匹配道，时间分辨率高的地震记录作为参考道，通过Curvelet变换的匹配滤波方法对匹配道地震记录以及掺加高斯白噪声的地震记录进行匹配处理，并与参考道分辨率高的地震记录进行对比。

基于Curvelet变换的地震信号去噪方法，包括以下步骤：

A、输入含噪声信号道集；

B、利用Curvelet变换，由原始地震数据得到匹配算子A；

C、利用匹配算子A构建等式Ax＝b+n……(1)，

在匹配滤波过程中：A为匹配算子，x为待匹配数据，b为匹配道数据，n为随机噪声；

在压制随机噪声过程中：A表示为产生随机噪声的算子矩阵，x是不含随机噪声的数据，b为伴随随机噪声的地震记录；以上两种情况的n都是未知的噪声数据；

D、将等式(1)中数据转换到Curvelet域变为稀疏数据；

E、将等式转化为目标函数： $J={\left|\right|\mathrm{Dm}-b\left|\right|}_2^2+\mathrm{\λ}{\left|\right|m\left|\right|}_1\·\·\·\left(2\right),$

式中：m是x在Curvelet域中稀疏系数，D＝AR，R是x的Curvelet逆变换，即x＝Rm，λ是正则化参数；

F、令目标函数j最小，利用FISTA算法求解x；

G、输出滤波后的地震数据。

有益效果：建立一个不同子波长度的楔形模型，对无噪声和对含高斯白噪声的地震数据使用FISTA算法进行Curvelet域精细尺度下地震数据匹配,经实际对比试验，匹配后的结果波形一致性高，并且噪声得到了好的压制，明显的提高了地震数据处理的分辨率、信噪比。通过实验结果对比可以看出，匹配滤波效果好，可以有效地压制噪声，提高了地震数据信噪比、分辨率，增强了地震勘探效果。

附图说明：

图1为基于Curvelet变换的地震信号去噪流程图。

图2(a)为模拟数据1波形图，由最小相位子波合成的褶积模型构成。

图2(b)是对模拟数据1进行频谱分析得到的振幅谱。

图2(c)是对模拟数据1进行频谱分析得到的相位谱。

图3(a)为模拟数据2波形图，由零相位子波合成的褶积模型构成。

图3(b)是对模拟数据2进行频谱分析得到的振幅谱。

图3(c)是对模拟数据2进行频谱分析得到的相位谱。

图4(a)为模拟数据1在Curvelet域经过FISTA计算匹配得到的地震波形图。

图4(b)是对图4(a)波形图进行频谱分析得到振幅谱图。

图4(c)是对图4(a)波形图进行频谱分析得到相位谱图。

图5(a)为模拟数据3波形图，由模拟数据1加入一定量的高斯白噪声构成，信噪比为-7.45dB。

图5(b)为模拟数据3进行频谱分析得到的振幅谱图。

图5(c)为模拟数据3进行频谱分析得到的相位谱图。

图6(a)为模拟数据4波形图，由模拟数据2加入一定量的高斯白噪声构成，信噪比为1.09dB。

图6(b)为模拟数据4进行频谱分析得到的振幅谱图。

图6(c)为模拟数据4进行频谱分析得到的相位谱图。

图7(a)是对加入噪声的模拟数据3进行Curvelet域匹配滤波的结果。

图7(b)是对图7(a)波形图进行频谱分析得到振幅谱图。

图7(c)是对图7(a)波形图进行频谱分析得到相位谱图。

具体实施方式：

下面结合附图和实施例对本发明作进一步详细的说明：

如图1是基于Curvelet变换的地震信号去噪流程图，

基于Curvelet变换的地震信号去噪方法，包括以下步骤：

A、输入含噪声信号道集；

B、利用Curvelet变换，由原始地震数据得到匹配算子A；

C、利用匹配算子A构建等式Ax＝b+n………(1)，

在匹配滤波过程中：A为匹配算子，x为待匹配数据，b为匹配道数据，n为随机噪声；

在压制随机噪声过程中：A表示为产生随机噪声的算子矩阵，x是不含随机噪声的数据，b为伴随随机噪声的地震记录；以上两种情况的n都是未知的噪声数据；

D、将等式(1)中数据转换到Curvelet域变为稀疏数据；

E、将等式转化为目标函数： $J={\left|\right|\mathrm{Dm}-b\left|\right|}_2^2+\mathrm{\λ}{\left|\right|m\left|\right|}_1\·\·\·\left(2\right),$

式中：m是x在Curvelet域中稀疏系数，D＝AR，R是x的Curvelet逆变换，即x＝Rm，λ是正则化参数；

F、令目标函数j最小，利用FISTA算法求解x；

G、输出滤波后的地震数据。

通过Curvelet变换，由原始地震数据得到匹配算子A。

再利用匹配算子A构建等式Ax＝b+n，在匹配滤波过程中：A匹配滤波算子矩阵，x为待匹配数据，b为匹配道数据；在压制随机噪声过程中：A表示为产生随机噪声的算子矩阵，x是不含随机噪声的数据，b为伴随随机噪声的地震记录；以上两种情况下。n都是未知的噪声数据；

对等式(1)中数据进行Curvelet变换，转换到Curvelet域变为稀疏数据；

离散Curvelet变换原理如下：

定义一对窗函数：径向窗函数W(r)和角度窗函数V(t)，其中，r∈(1/2,2),t∈[-1,1]。它们都满足可允许条件：

$\underset{j=-\∞}{\overset{\∞}{\mathrm{\Σ}}}{W}^2\left(2^{j}r\right)=1,r\∈(3/\mathrm{4,3}/2)\·\·\·\left(3\right)$

$\underset{j=-\∞}{\overset{\∞}{\mathrm{\Σ}}}{V}^2(t-l)=1,t\∈(-1/\mathrm{2,1}/2)\·\·\·\left(4\right)$

对于每一个j≥j₀，在频域中定义频窗U_j如下：

$U_j(r,\mathrm{\θ})=2^{-3j/4}W\left(2^{j}r\right)V\left(\frac{2^{[j/2]}\mathrm{\θ}}{2\mathrm{\π}}\right)\·\·\·\left(5\right)$

其中[j/2]是j/2的整数部分。U_j的支撑区间受W和V支撑区间限制所获得的楔形区域，楔形区域符合各向异性尺度的特征。假设原子φ_j的频率表达式为φ_j(ω)＝U_j(ω)，若知道了尺度j上的φ_j，则其它2^-j尺度上的Curvelet均可以由φ_j通过旋转和平移获得。设等间隔的旋转角序列θ_l＝2π·2^-[j/2]·l，l为非负整数，0≤θ≤2π；平移参数为：k＝(k₁,k₂)∈Z²。

综合以上概念，定义尺度在2^-j，方向θ_l，平移参数(k₁,k₂)处的Curvelet为：

${\mathrm{\φ}}_{j,l,k}\left(x\right)={\mathrm{\φ}}_j\left(R_{{\mathrm{\θ}}_l}(x-x_k^{(j,k)})\right)\·\·\·\left(6\right)$

其中 $x_k^{(j,l)}=R_{\mathrm{\θl}}^{-1}(k_1\·2^{-j},k_2\·2^{-j/2}),$ R_θ和θ_l旋转获得。Curvelet变换可以表示为：

$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)}\mathrm{dx}\&CenterDot;\&CenterDot;\&CenterDot;\left(7\right)$

通过笛卡尔坐标系下的f[t₁,t₂],0≤t₁,t₂＜n为信号的输入，则Curvelet变换的离散形式为：

$c[j,l,k]=\underset{0\&le;t_1,t_2<n}{\mathrm{\&Sigma;}}f[t_1,t_2]\stackrel{\&OverBar;}{{\mathrm{\&phi;}}_{j,l,k}[t_1,t_2]}\&CenterDot;\&CenterDot;\&CenterDot;\left(8\right)$

其中φ_j,l,k[t₁,t₂]为离散Curvelet波形。

通常情况下受勘探条件以及噪声的影响，会使得线性问题求解存在多解性，需要引入正则化方法对解进行约束。从信息论的角度来看，最优匹配处理时一个欠定问题，Chen(2001)等则在理论上证明了在一定条件下，用稀疏建模具有最佳的重构效果，因此，为降低多解性，采用稀疏约束反演法

将数据转换到Curvelet域变为稀疏数据后，将方程式(1)转化为下式，只需将目标函数J最小，便可求得解：

$J={\left|\right|\mathrm{Dm}-b\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|m\left|\right|}_1\&CenterDot;\&CenterDot;\&CenterDot;\left(9\right)$

其中，m是公式(1)中x在Curvelet域中的稀疏系数，D＝AR，R是x的Curvelet逆变换，即，x＝Rm。λ是正则化参数，其控制L1范数的约束项与L2范数的误差项之间的权重。方程式(9)是curvelet域中待求解的稀疏系数。由于方程式(9)的D可以用AR表示，可以将方程式写成如下式子：

$J={\left|\right|\mathrm{ARm}-b\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|m\left|\right|}_1\&CenterDot;\&CenterDot;\&CenterDot;\left(10\right)$

令目标函数j最小，利用FISTA算法对方程式(2)进行迭代运算求解，得到待匹配数据x。

采用FISTA(Fast Iterative Shrinkage-Thresholding Algorithm)(Berkand Teboulle,2009)求解方程式(10)。FISTA作为一种有效的解决线性反演问题的算法，具有收敛速度快，迭代次数少的优点，大大提高计算效率。FISTA求解过程如下：

第1步：b₁＝x₀和t₁＝0

第k步：

$\begin{array}{c}{x}_k=\mathrm{soft}(x_{k-1}+\frac{1}{\mathrm{\&alpha;}}{A}^H(b-{\mathrm{Ax}}_{k-1}),\frac{1}{\mathrm{\&alpha;}}),\\ t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2},\\ b_{k+1}=x_k+\frac{t_{k-1}}{t_{k+1}}(x_k-x_{k-1}).\end{array}\&CenterDot;\&CenterDot;\&CenterDot;\left(11\right)$

其中，soft(x,t)定义为软阈值，当输入x为复数时，令x＝ze^iω。软阈值可表示为：

$\mathrm{soft}(x,t)=\left\{\begin{array}{ccc}(z-t)e^{\mathrm{i\&omega;}},& \mathrm{if}& z>t\\ 0,& \mathrm{if}& z\&le;t\end{array}\right.\&CenterDot;\&CenterDot;\&CenterDot;\left(12\right)$

FISTA算法具有实现简单，收敛速度快，并能将复数的运算也包括在内，扩大了使用范围。

不同种类地震子波合成的楔形模型代表不同震源或者年代的地震记录。子波波长较长，时间分辨率较低，反之亦然。不同子波分辨薄层和楔形的能力也不同，文中给出两种不同类型子波的褶积模型。为了更好地模拟真实情况以证明本方法的有效性，本文建立一个不同子波长度的楔形模型，且模型中存在一个薄层。本发明中选用时间分辨率低的地震记录作为匹配道，时间分辨率高的地震记录作为参考道。

模型共有45道，采样时间2ms，零相位和最小相位子波长度均为140ms，模型厚度是4000ms，薄层厚度是120ms。对模型进行Curvelet域多尺度分解，而只分解了三个尺度的Curvelet变换。

如图2是最小相位子波合成的褶积模型，图2(b)和图2(c)分别是对最小相位地震子波的褶积模型频谱分析得到的振幅谱和相位谱，其主频范围10-80Hz。由于子波类型影响，褶积模型中分辨楔形和薄层的能力大大降低，给后续处理与解释带来了不便，因此将最小相位子波获得褶积模型作为匹配道。图3(a)为模拟数据2波形图，由零相位子波合成的褶积模型构成。对模拟数据2频谱分析得到振幅谱图3(b)和相位谱图3(c)。由频谱分析可知，地震波主频能量也集中在10-80Hz左右，相位谱3(c)是一个非线性的谱线，与相位谱图2(c)相比相位存在较大的差别。模拟数据2具有较高的时间分辨率，为地震数据的后续处理和解释提供了便利，因而将模拟数据2作为参考地震数据。图4(a)是Curvelet域经过FISTA计算匹配得到的地震波形图；对匹配滤波结果进行频谱分析得到振幅谱图4(b)和相位谱图4(c)。与参考道地震数据2进行波形、振幅谱及相位谱对比分析，匹配后地震数据波形、振幅及相位一致性很高，说明Curvelet域最优约束范数匹配具有很好的效果。

在实际地震勘探中，地震数据往往存在噪声。为更好的模拟实际情况，本发明对模拟数据1和模拟数据2加入一定量的高斯白噪声。对模拟数据1加入噪声称为模拟数据3，图5(a)是模拟数据3的波形图，信噪比为-7.45dB。对模拟数据3进行频谱分析得到振幅谱图5(b)和相位谱图5(c)。对模拟数据2加入高斯白噪声得到模拟数据4，信噪比为1.09dB，如图6(a)所示。对加入噪声的模拟数据进行频谱分析得到振幅谱图6(b)和相位谱图6(c)。由于信噪比高和零相位子波类型，其分辨楔形和薄层能力较强，因而将模拟数据4作为参考道地震数据。图7(a)是对加入噪声的模拟数据3进行Curvelet域匹配滤波的结果，对匹配后地震数据经过频谱分析得到振幅谱7(b)和相位谱7(c)。与模拟数据4对比，波形具有很高的相似性，信噪比提高到26dB。对比振幅谱和相位谱,它们也具很较高的一致性。

Curvelet变换都具有很高的稀疏性，对地震数据进行多尺度分解，为信号和噪声的分离奠定了基础，因而利用稀疏约束求解匹配问题更加准确。本发明方法取得了更好的效果，为匹配处理提供了新的思路,将最优范数和Curvelet分解结合在一起，为分辨薄层和新老地震数据提供了有价值的信息。FISTA算法具有较快的计算效率，且具收敛性能好等特点，易于计算机编程等特点,能够很好地应用于数据拟合的测量实践中。Curvelet域匹配滤波与噪声压制的结果具有同相轴连续性好、波形、振幅谱和相位谱一致性高的特点。

Claims (2)

1.一种基于Curvelet变换的地震信号去噪方法，其特征在于，建立一个不同子波长度的楔形模型，模型中只存在一个薄层，不同子波分辨薄层和楔形的能力也不同，选用时间分辨率低的地震记录作为匹配道，时间分辨率高的地震记录作为参考道，通过Curvelet变换的匹配滤波方法对匹配道地震记录以及掺加高斯白噪声的地震记录进行匹配处理，并与参考道分辨率高的地震记录进行对比。

2.按照权利要求1所述的基于Curvelet变换的地震信号去噪方法，其特征在于，包括以下步骤：

A、输入含噪声信号道集；

B、利用Curvelet变换，由原始地震数据得到匹配算子A；

C、利用匹配算子A构建等式Ax＝b+n…(1)，

在匹配滤波过程中：A为匹配算子，x为待匹配数据，b为匹配道数据，n为随机噪声；

在压制随机噪声过程中：A表示为产生随机噪声的算子矩阵，x是不含随机噪声的数据，b为伴随随机噪声的地震记录；以上两种情况的n都是未知的噪声数据；

D、将等式(1)中数据转换到Curvelet域变为稀疏数据；

E、将等式转化为目标函数： $J={\left|\right|\mathrm{Dm}-b\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|m\left|\right|}_1...\left(2\right),$

式中：m是x在Curvelet域中稀疏系数，D＝AR，R是x的Curvelet逆变换，即x＝Rm，λ是正则化参数；

F、令目标函数j最小，利用FISTA算法求解x；

G、输出滤波后的地震数据。