CN111399057A - 一种基于非凸稀疏约束的地震资料噪声压制方法 - Google Patents

Abstract

本发明涉及一种基于非凸稀疏约束的地震资料噪声压制方法，包括以下步骤：1)获得原始地震资料，并对原始地震资料进行GBW小波变换，确定其标架算子和伴随算子；2)构造基于GBW和非凸稀疏约束的随机噪声压制模型；3)求解步骤2)中的随机噪声压制模型，获得去噪后的结果；4)利用小波变换的反变换公式，将去噪后的结果转换到时域内，得到原始地震资料去噪后的时域信号。通过合成数据和实际数据对比，本发明提出的基于非凸稀疏约束的小波域随机噪声压制方法能够在保证地震资料的连续性时，最大限度的压制随机噪声。

Description

一种基于非凸稀疏约束的地震资料噪声压制方法

技术领域

本发明属于地震勘探技术领域，涉及一种叠前地震资料随机噪声压制的方法，尤其是涉及一种基于非凸稀疏约束的地震资料噪声压制方法，利用GBW(广义Beta小波)和非凸稀疏约束来压制叠前地震资料的随机噪声。

背景技术

随着采集系统和仪器的发展、数据采集道数的增加为地震资料高分辨率提供了基础，使得地震数据可以得到比旅行时和地质构造更多的信息。为了有效地利用这些信息，必须要求更高的信噪比和保真度。例如，储层属性估计和反演对叠前数据带宽、振幅和相位保真都有更高要求。由于复杂的地表和地形条件，如山地、黄土塬、崎岖海底、沙漠、戈壁和多期叠合构造、逆掩推覆、强烈褶皱等地质因素对地震信号产生各种畸变和干扰。这在地震记录上同样会表现为信噪比低、噪声类型复杂、有效信号能量弱的特点，同时将严重制约目标信号的准确成像和归位。因此，为了能更充分和准确地应用采集的地震资料，噪声衰减问题仍是地震资料处理的关键问题之一。

1984年Canales等人提出在fourier域的F-x预测滤波的方法。2013年王平和高静怀提出了在小波域利用小波系数的稀疏性来去除地震资料中的随机噪声。2016年Siahsar等人提出在同步挤压域利用SVD分解的方法来压制地震资料中的随机噪声等等。后续还有人研究在seislet以及curvelet域对随机噪声进行压制的方法。这些方法基本都是在变换域中进行噪声压制的。在变换域中，有效信号的能量分布在有限的空间内，但是随机噪声分布在整个变化域空间里。

尽管许多新的时频变换已经应用于噪声衰减方法研究中，但是小波变换依然在噪声衰减方法研究中具有重要的意义。小波变换可以构成紧标架来降低计算速度。根据标架理论，小波域的噪声衰减方法可以表示为一个带L1稀疏约束的最小二乘问题。但是，L1稀疏约束的结果容易低估高幅度值，从而导致去噪结果不准确。

综上，以上技术具有如下缺点：

1、Morlet母小波与地震子波不匹配，导致求解出来的小波变换时频分辨率较低；

2、L1稀疏约束的结果容易低估高幅度值，从而导致去噪结果不准确。

发明内容

针对上述问题，本发明的目的是提供一种基于非凸稀疏约束的地震资料噪声压制方法，该方法引入GBW母小波，从而更好的匹配地震子波，获得更加稀疏的小波系数。同时，为了获得更加稀疏的、准确的去噪结果，该方法引入非凸的稀疏约束，与带有L1正则化约束的去噪方法对比，本发明提出的去噪方法能够获得更高的SNR(信噪比)。

为实现上述目的，本发明采取以下技术方案：一种基于非凸稀疏约束的地震资料噪声压制方法，其包括以下步骤：

1)获得原始地震资料s，并对原始地震资料s进行GBW小波变换，确定其标架算子K和伴随算子K^*；

2)基于确定的标架算子K和伴随算子K^*，构造基于GBW和非凸稀疏约束的随机噪声压制模型；

3)求解步骤2)中的随机噪声压制模型，获得去噪后的结果

4)利用小波变换的反变换公式，将去噪后的结果

转换到时域内，得到原始地震资料去噪后的时域信号

进一步的，所述步骤1)中，所述标架算子K和伴随算子K^*的计算公式为：

x＝<s,ψ_m,n>＝K^*s

s＝Kx

其中，ψ_m,n表示小波族；<·,·>表示二者的内积；x表示有效信号在小波域的系数。

进一步的，所述步骤2)中，所述随机噪声压制模型为：

其中，F(x)为优化问题的目标函数；K为GBW小波变换的标架算子；x为有效信号在小波域的系数；B为给定的矩阵；λ为规则化参数；s为原始地震信号；v为中间变量；ψ_B(x)为非凸的稀疏约束。

进一步的，所述步骤3)中，求解步骤2)中的随机噪声压制模型，获得去噪后的结果

的方法，包括以下步骤：

3.1)确定输入变量，包括原始地震信号s和标架算子K；

3.2)确定迭代参数；

3.3)对有效信号在小波域的系数x和中间变量v进行初始化，并根据确定的迭代参数对二者进行更新；

3.4)计算相邻两次迭代过程中有效信号在小波域的系数的差值，判断其差值是否小于预设阈值，如果小于预设阈值或达到迭代次数则迭代结束，否则进入步骤3.5)；

3.5)对正则化参数γ进行更新，并返回步骤3.3)，直到得到有效信号在小波域的最优系数

进一步的，所述步骤3.2)中，迭代参数包括正则化参数γ，且0<γ<1；中间变量μ，且

最大迭代次数N；正则化参数的范围[λ_min,λ_max]：其中，λ_min为正则化参数的最小值，λ_max为正则化参数的最大值。

进一步的，所述步骤3.3)中，更新公式为：

x^k+1＝x^k-μK^*(K(x^k+γ(v^k-x^k))-s)

v^k+1＝v^k-μγK^*K(v^k-x^k)

x^k+1＝soft(x^k+1,λ^k+1)

v^k+1＝soft(v^k+1,λ^k+1)

其中，x^k和x^k+1分别为第k次和第k+1次迭代过程中有效信号在小波域的系数；v^k和v^k+1分别为第k次和第k+1次迭代过程的中间变量；μ为中间变量；γ为正则化参数；K和K^*分别为标架算子和伴随算子；soft(·,·)是软阈值策略，定义为：

进一步的，所述步骤3.4)中，差值计算公式为：

式中，x^k和x^k+1分别为第k次和第k+1次迭代过程中有效信号在小波域的系数；error为两次迭代过程中的系数差值，tolerance为预设阈值。

进一步的，所述步骤3.5)中，正则化参数λ的选择方法为：

其中，λ_max和λ_min分别是最大和最小正则化参数；N是最大迭代次数；λ^k是第k次迭代的正则化参数。

本发明由于采取以上技术方案，其具有以下优点：1、本发明引入GBW母小波，从而更好的匹配地震子波，获得更加稀疏的小波系数。2、本发明为了获得更加稀疏的、准确的去噪结果，该方法引入非凸的稀疏约束，与带有L1正则化约束的去噪方法对比，本研究提出的去噪方法能够获得更高的SNR。3、本发明的算法内容易于实现，可操作性好；并且能够有效的压制地震资料的随机噪声，同时可以保护有效信号的结构。因此，本发明可以广泛应用于地震勘探技术领域。

附图说明

图1(a)～图1(d)是合成记录，其中，图1(a)是无噪合成地震数据；图1(b)是含噪合成地震数据，SNR＝-4dB；图1(c)是无噪数据的GBW变换的时频谱；图1(d)是含噪数据的GBW变换的时频谱；

图2(a)～图2(d)是去噪结果：图2(a)是本发明提出的算法；图2(b)是基于GMC和morlet小波的随机噪声压制结果；图2(c)是基于IST和GBW的随机噪声压制结果；图2(d)是基于FIST和GBW的随机噪声压制结果；

图3是不同方法每次迭代输出数据的SNR值对比；

图4是实际2D地震数据，该数据含有200道和750个时间采样点；

图5(a)～图5(c)是实际2D地震数据的去噪结果：图5(a)是本发明提出的算法；图5(b)是基于FIST和GBW的随机噪声压制方法；图5(c)是小波阈值方法；

图6(a)～图6(c)是实际2D地震数据的去噪差剖面：图6(a)是本发明提出的；图6(b)基于FIST和GBW的随机噪声压制方法；图6(c)小波阈值方法。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

随着地震勘探的不断深入，勘探目标的地表条件和地质条件越来越复杂，使得所采集数据中包含了更加复杂的噪声。如何有效地衰减地震资料中的随机噪声和相干噪声，突出有效信号是地震勘探数据处理中的一个重要的研究课题。本发明提出一种基于非凸稀疏约束的地震资料噪声压制方法，包括以下步骤：

1)获取原始地震资料s，并对原始地震资料s进行GBW小波变换，确定其标架算子K和伴随算子K^*；

2)基于确定的标架算子K和伴随算子K^*，构造基于GBW和非凸稀疏约束的随机噪声压制模型；

3)利用forward-backward(FB)方法对步骤2)中构造的随机噪声压制模型进行求解，获得去噪后的结果

4)利用小波变换的反变换公式，将去噪后的结果

转换到时域内，获得原始地震资料去噪后的时域信号

上述步骤1)中，具体包括以下步骤：

1)引入GBW小波，其定义为：

其中，U(ω)为单位阶跃函数，ω为频率；α和β为GBW小波的自由参数，且α>0,β>0；a_α,β为归一化常量。GBW小波具有严格的解析性，通过调整两个自由参数，可以展示出灵活多变的特点。

2)根据GBW小波变换公式对获取的原始地震资料s进行GBW小波变换，确定其标架算子K和伴随算子K^*，其定义为：

其中，ψ_m,n为小波族；<·,·>为二者的内积；x为有效信号在小波域的系数，其中，有效信号也就是理论上原始地震资料s去噪后的真实解。

上述步骤2)中，构造基于GBW和非凸稀疏约束的随机噪声压制模型时，当选用合适的小波函数对含噪信号进行小波变换，将其投影到时间尺度域的时候，有效信号的能量将会分布在较小的子空间，噪声的能量会扩散到比较大的子空间，甚至整个时间-尺度域。换句话说，有效信号和少数的系数相对应，而噪声几乎分布在全部的系数上。如果我们确定了有效信号对应的子空间，得到有效信号对应的系数，将其它系数置零，则噪声就会在变换域被压制，信噪比得到提高。因此，可以根据反问题思想将小波域去噪表示为一个反问题求解问题，本发明基于该思想，建立的基于GBW和非凸稀疏约束的随机噪声压制模型为：

其中，F(x)为优化问题的目标函数；K为GBW小波变换的标架算子；x为有效信号在小波域的系数；B为给定的矩阵；λ为规则化参数；s为原始地震资料；v为中间变量。ψ_B(x)为非凸的稀疏约束。为了方便后续优化算法的求解，必须使得目标函数F(x)为凸的，因此矩阵B满足

上述步骤3)中，在求解基于GBW和非凸稀疏约束的随机噪声压制模型时，引入指数衰减策略选择正则化参数。

在迭代的初始阶段，求得

与真实解x差别很大，因此，目标函数中的误差项

将非常大，这时候，应该充分利用先验信息以快速得到近似的解。而先验信息正是通过规则化项体现出来，规则化参数λ给出了第一个误差项与先验约束的折中。这意味着非常有必要在迭代的初始阶段设定一个较大的权重。具体计算方法包括以下步骤：

3.1)确定输入变量，包括原始地震信号s和标架算子K；

3.2)确定迭代参数，迭代参数包括正则化参数γ，且0<γ<1；中间变量μ，且

最大迭代次数N；正则化参数的范围[λ_min,λ_max]：其中，λ_min为正则化参数的最小值，λ_max为正则化参数的最大值；

3.3)对有效信号在小波域的系数x和中间变量v进行初始化，并根据确定的迭代参数进行更新，更新公式为：

x^k+1＝x^k-μK^*(K(x^k+γ(v^k-x^k))-s) (4)

v^k+1＝v^k-μγK^*K(v^k-x^k) (5)

x^k+1＝soft(x^k+1,λ^k+1) (6)

v^k+1＝soft(v^k+1,λ^k+1) (7)

其中，x^k和x^k+1分别为第k次和第k+1次迭代过程中有效信号在小波域的系数；v^k和v^k+1分别为第k次和第k+1次迭代过程的中间变量；μ为中间变量；γ为正则化参数；K和K^*分别为标架算子和伴随算子。其中，soft(·,·)是软阈值策略，定义为：

3.4)计算相邻两次迭代过程中有效信号在小波域的系数的差值，判断其差值是否小于预设阈值，如果小于预设阈值或达到迭代次数则迭代结束，否则进入步骤3.5)；

差值计算公式为：

式中，error为两次迭代过程中的系数差值，tolerance为预设阈值。

3.5)对正则化参数γ进行更新，并返回步骤3.3)，直到得到最优系数

随着迭代过程的进行，

越来越接近真实解，此时，先验信息的重要性逐渐降低，而误差项的权重应该逐步变大，因此应该选定更小的规则化参数λ，λ选择方法为：

其中，λ_max和λ_min分别是最大和最小正则化参数；N是最大迭代次数；λ^k是第k次迭代的正则化参数。

实施例一

合成地震记录数据

如图1(a)～图1(d)所示，为了验证本发明提出方法对地震资料随机噪声压制的有效性，本实施例选取由35Hz的Ricker子波构成的合成记录。其中，图1(a)为无躁合成地震数据，图1(b)是在无躁合成数据中加入随机噪声，含噪合成数据的SNR为-4dB，图1(c)是无噪数据的GBW变换的时频谱；图1(d)是含噪数据的GBW变换的时频谱。本实施例同时与基于GMC和morlet小波的随机噪声压制结果、基于IST和GBW的随机噪声压制结果、基于FIST和GBW的随机噪声压制结果去噪方法进行对比。

如图2(a)～图2(d)所示，分别是不同去噪方法去噪后的结果。从图2(a)～图2(d)中可以发现，相对于其他五种去噪方法，本发明提出的方法能够有效的压制合成数据中的随机噪声。

如图3所示，是每次迭代过程中输出信号的SNR值。从图3中可以看出，本发明提出的去噪方法需要的迭代次数更少，并且获得的SNR值更高。

实施例二

实际地震记录数据

如图4所示，本实施例选择实际地震数据进一步验证本发明去噪方法的有效性。图4所示的实际地震数据包含200道，且每道数据有750个时间采样点，时间采样间隔为2ms。

如图5(a)～图5(c)所示，本实施例中，只将本发明提出的方法与基于FIST和GBW的随机噪声压制和小波阈值方法两个去噪方法进行对比，可以看出，本发明方法相较于其他两种方法去噪结果更好。

如图6(a)～图6(c)所示，为去噪后的差剖面。通过图5(a)～图5(c)和图6(a)～图6(c)进行对比可得出，本发明提出的方法能够有效的压制随机噪声且对有效信号的影响最小。

以上给出一种具体的实施方式，但本发明不局限于所描述的实施方式。本发明的基本思路在于上述方案，对本领域普通技术人员而言，根据本发明的教导，设计出各种变形的模型、公式、参数并不需要花费创造性劳动。在不脱离本发明的原理和精神的情况下对实施方式进行的变化、修改、替换和变形仍落入本发明的保护范围内。

Claims (7)

1.一种基于非凸稀疏约束的地震资料噪声压制方法，其特征在于包括以下步骤：

1)获得原始地震资料s，并对原始地震资料s进行GBW小波变换，确定其标架算子K和伴随算子K^*；

2)基于确定的标架算子K和伴随算子K^*，构造基于GBW和非凸稀疏约束的随机噪声压制模型；

3)求解步骤2)中的随机噪声压制模型，获得去噪后的结果

4)利用小波变换的反变换公式，将去噪后的结果

转换到时域内，得到原始地震资料去噪后的时域信号

2.如权利要求1所述的一种基于非凸稀疏约束的地震资料噪声压制方法，其特征在于：所述步骤1)中，所述标架算子K和伴随算子K^*的计算公式为：

x＝<s,ψ_m,n>＝K^*s

s＝Kx

其中，ψ_m,n表示小波族；<·,·>表示二者的内积；x表示有效信号在小波域的系数。

3.如权利要求1所述的一种基于非凸稀疏约束的地震资料噪声压制方法，其特征在于：所述步骤2)中，所述随机噪声压制模型为：

其中，F(x)为优化问题的目标函数；K为GBW小波变换的标架算子；x为有效信号在小波域的系数；B为给定的矩阵；λ为规则化参数；s为原始地震信号；v为中间变量；ψ_B(x)为非凸的稀疏约束。

4.如权利要求1所述的一种基于非凸稀疏约束的地震资料噪声压制方法，其特征在于：所述步骤3)中，求解步骤2)中的随机噪声压制模型，获得去噪后的结果

的方法，包括以下步骤：

3.1)确定输入变量，包括原始地震信号y和标架算子K；

3.2)确定迭代参数，迭代参数包括正则化参数γ，且0<γ<1；中间变量μ，且

最大迭代次数N；正则化参数的范围[λ_min,λ_max]：其中，λ_min为正则化参数的最小值，λ_max为正则化参数的最大值；

3.3)对有效信号在小波域的系数x和中间变量v进行初始化，并根据确定的迭代参数对二者进行更新；

3.4)计算相邻两次迭代过程中有效信号在小波域的系数的差值，判断其差值是否小于预设阈值，如果小于预设阈值或达到迭代次数则迭代结束，否则进入步骤3.5)；

3.5)对正则化参数γ进行更新，并返回步骤3.3)，直到得到有效信号在小波域的最优系数

5.如权利要求4所述的一种基于非凸稀疏约束的地震资料噪声压制方法，其特征在于：所述步骤3.3)中，更新公式为：

x^k+1＝x^k-μK^*(K(x^k+γ(v^k-x^k))-s)

v^k+1＝v^k-μγK^*K(v^k-x^k)

x^k+1＝soft(x^k+1,λ^k+1)

v^k+1＝soft(v^k+1,λ^k+1)

其中，x^k和x^k+1分别为第k次和第k+1次迭代过程中有效信号在小波域的系数；v^k和v^k+1分别为第k次和第k+1次迭代过程的中间变量；μ为中间变量；γ为正则化参数；K和K^*分别为标架算子和伴随算子；soft(·,·)是软阈值策略，定义为：

6.如权利要求4所述的一种基于非凸稀疏约束的地震资料噪声压制方法，其特征在于：所述步骤3.4)中，差值计算公式为：

式中，x^k和x^k+1分别为第k次和第k+1次迭代过程中有效信号在小波域的系数；error为两次迭代过程中的系数差值，tolerance为预设阈值。

7.如权利要求4所述的一种基于非凸稀疏约束的地震资料噪声压制方法，其特征在于：所述步骤3.5)中，正则化参数λ的选择方法为：

其中，λ_max和λ_min分别是最大和最小正则化参数；N是最大迭代次数；λ^k是第k次迭代的正则化参数。

Images