CN104730576A

CN104730576A - Curvelet transform-based denoising method of seismic signals

Info

Publication number: CN104730576A
Application number: CN201510175579.XA
Authority: CN
Inventors: 龙云; 王洪超; 陈祖斌; 林君; 杨泓渊
Original assignee: Jilin University
Current assignee: Jilin University
Priority date: 2015-04-14
Filing date: 2015-04-14
Publication date: 2015-06-24

Abstract

The present invention relates to a seismic signal denoising method based on Curvelet transformation. A wedge-shaped model with different wavelet lengths is established. There is only one thin layer in the model, and different wavelets have different abilities to distinguish thin layers and wedges. The time resolution is selected. The low seismic records are used as matching traces, and the seismic records with high time resolution are used as reference traces. Matching processing of matching trace seismic records and seismic records mixed with Gaussian white noise is carried out through the matched filtering method of Curvelet transform, and the resolution of the reference traces is compared with that of the reference traces. High seismic records for comparison. The actual comparative test results show that the FITSA algorithm processes the influences of the two aspects at the same time, the waveform consistency of the matching result is high, and the noise is well suppressed, which significantly improves the resolution and signal-to-noise ratio of seismic data processing. While improving the signal-to-noise ratio, it greatly improves the processing efficiency of seismic data.

Description

Seismic Signal Denoising Method Based on Curvelet Transform

技术领域technical field

本发明涉及一种地震数据的处理方法，尤其是基于Curvelet变换的地震信号去噪方法。The invention relates to a method for processing seismic data, in particular to a seismic signal denoising method based on Curvelet transformation.

技术背景technical background

地震勘探中，不同震源的地震数据受匹配滤波方法的限制，匹配滤波后的地震数据仍存在信噪比低、分辨率低等一系列差异；或新老地震资料进行匹配处理后出现分层不清晰，使得波阻抗反演的精度难以保证和地震属性参数的提取困难；或不同区块地震资料的拼接也存在同相轴不连续、信噪比低等，给地震数据匹配处理带来了困难，这一系列问题都严重制约着地震数据处理的分辨率、信噪比。In seismic exploration, the seismic data of different sources are limited by the matched filtering method, and there are still a series of differences in the seismic data after the matched filtering, such as low signal-to-noise ratio and low resolution; clear, which makes it difficult to guarantee the accuracy of wave impedance inversion and the extraction of seismic attribute parameters; or there are discontinuous events and low signal-to-noise ratio in the splicing of seismic data from different blocks, which brings difficulties to seismic data matching processing. These series of problems seriously restrict the resolution and signal-to-noise ratio of seismic data processing.

吉林大学2008年硕士论文公开了《基于Curvelet变换的地震数据去噪方法研究》，介绍了一种基于Curvelet变换的地震数据处理技术研究；将第二代Curvelet变换引入地震数据处理,通过选取恰当的阈值,对地震数据进行Curvelet变换。计算了模拟数据的Curvelet变换,对比了不同信噪比下的去噪效果。并且通过Curvelet变换对实际地震数据进行去噪处理，实验结果得到了比较好的效果。Jilin University’s 2008 master’s thesis published "Research on Seismic Data Denoising Method Based on Curvelet Transform", which introduced a seismic data processing technology research based on Curvelet transform; the second-generation Curvelet transform was introduced into seismic data processing. Threshold, Curvelet transformation is performed on the seismic data. The Curvelet transform of the simulated data is calculated, and the denoising effect under different signal-to-noise ratios is compared. In addition, the actual seismic data is denoised by Curvelet transformation, and the experimental results have obtained relatively good results.

CN103543469A公开了《一种基于小波变换的小尺度阈值去噪方法》，首先，小尺度地对地震数据进行扫描，得到一个小时窗内的相关系数值，其次，设置一个阈值对这个小尺度时窗内的数据进行判断，是以地震信号为主还是以噪声信号为主，然后采用合适的分频数及合适的小波阈值进行去噪，最后将去噪后的小波尺度进行小波重构，从而得到去噪后信噪比较高的地震道集。CN103543469A discloses "A Small-Scale Threshold Denoising Method Based on Wavelet Transform". First, the seismic data is scanned on a small scale to obtain the correlation coefficient value within an hour window. Secondly, a threshold is set for this small-scale time window. It is judged whether it is mainly seismic signal or noise signal based on the data inside, and then the denoising is carried out with the appropriate frequency division number and appropriate wavelet threshold, and finally the wavelet scale after denoising is reconstructed by wavelet, so as to obtain Seismic gathers with high signal-to-noise ratio after denoising.

CN103399348A公开了一种《基于Shearlet变换的地震信号去噪方法》，读入二维地震剖面数据，将二维地震数据S扩展为一长宽为奇数的方阵S1，构造频域方向滤波器组，将各个变换矩阵分别与S1信号向量相乘，再分别做二维傅里叶反变换,得到各个方向和尺度的Shearlet系数C_i,j，阈值处理，将经过阈值处理的Shearlet变换系数做Shearlet反变换获得去噪之后的地震数据。本发明将含有噪声的地震信号做拉普拉斯分解然后利用剪切波函数进行滤波处理，得到对应剪切波系数，通过阈值处理滤除噪声信号，再通过逆非下采样剪切波变换恢复去噪的信号，获得较好的去噪效果，具有较好的实用价值。CN103399348A discloses a "Seismic Signal Denoising Method Based on Shearlet Transform", which reads in two-dimensional seismic profile data, expands the two-dimensional seismic data S into a square matrix S1 whose length and width are odd numbers, and constructs a frequency domain direction filter bank , each transformation matrix is multiplied by the S1 signal vector, and then the two-dimensional inverse Fourier transform is performed respectively to obtain the Shearlet coefficients C _i,j of each direction and scale, threshold value processing, and the Shearlet transformation coefficients processed by the threshold value are used as Shearlet Inverse transform to obtain the denoised seismic data. In the present invention, the noise-containing seismic signal is decomposed by Laplace and then filtered by the shear wave function to obtain the corresponding shear wave coefficient, and the noise signal is filtered out by threshold value processing, and then restored by the inverse non-subsampled shear wave transform The denoised signal can obtain a better denoising effect and has better practical value.

上述现有技术相比虽然在一定程度上能够有效的压制噪声，但过程复杂，计算速度较慢，信噪比低，地震数据处理效果不是很好。Compared with the above prior art, although the noise can be effectively suppressed to a certain extent, the process is complex, the calculation speed is slow, the signal-to-noise ratio is low, and the seismic data processing effect is not very good.

发明内容Contents of the invention

本发明的目的就是针对上述现有技术的不足，提出了将采集到的地震数据进行匹配滤波与噪声压制，通过FITSA算法同时处理两方面的影响，在改善信噪比的同时极大地提高地震资料处理效率。The purpose of the present invention is to address the above-mentioned deficiencies in the prior art, and proposes to carry out matched filtering and noise suppression on the collected seismic data, and simultaneously process the influence of the two aspects through the FITSA algorithm, and greatly improve the seismic data while improving the signal-to-noise ratio. Processing efficiency.

本发明的目的是通过以下方式实现的：The purpose of the present invention is achieved in the following manner:

建立一个不同子波长度的楔形模型，模型中只存在一个薄层，不同子波分辨薄层和楔形的能力也不同，选用时间分辨率低的地震记录作为匹配道，时间分辨率高的地震记录作为参考道，通过Curvelet变换的匹配滤波方法对匹配道地震记录以及掺加高斯白噪声的地震记录进行匹配处理，并与参考道分辨率高的地震记录进行对比。Establish a wedge-shaped model with different wavelet lengths. There is only one thin layer in the model, and different wavelets have different abilities to distinguish between thin layers and wedges. Seismic records with low temporal resolution are selected as matching traces, and seismic records with high temporal resolution are selected as matching traces. As a reference trace, match the matched trace seismic records and the seismic records mixed with Gaussian white noise by the matched filtering method of Curvelet transform, and compare them with the high-resolution reference trace seismic records.

基于Curvelet变换的地震信号去噪方法，包括以下步骤：The seismic signal denoising method based on Curvelet transform comprises the following steps:

A、输入含噪声信号道集；A. Input the noise-containing signal gather;

B、利用Curvelet变换，由原始地震数据得到匹配算子A；B. Using Curvelet transform, the matching operator A is obtained from the original seismic data;

C、利用匹配算子A构建等式Ax＝b+n……(1)，C. Utilize the matching operator A to construct the equation Ax=b+n...(1),

在匹配滤波过程中：A为匹配算子，x为待匹配数据，b为匹配道数据，n为随机噪声；In the matched filtering process: A is the matching operator, x is the data to be matched, b is the matching channel data, and n is random noise;

在压制随机噪声过程中：A表示为产生随机噪声的算子矩阵，x是不含随机噪声的数据，b为伴随随机噪声的地震记录；以上两种情况的n都是未知的噪声数据；In the process of suppressing random noise: A represents the operator matrix that generates random noise, x is the data without random noise, and b is the seismic record with random noise; n in the above two cases is unknown noise data;

D、将等式(1)中数据转换到Curvelet域变为稀疏数据；D, the data in the equation (1) is converted to the Curvelet domain and becomes sparse data;

E、将等式转化为目标函数： $J = {| | Dm - b | |}_{2}^{2} + λ {| | m | |}_{1} \cdot \cdot \cdot (2),$ E. Transform the equation into an objective function: $J = {| | Dm - b | |}_{2}^{2} + λ {| | m | |}_{1} \cdot \cdot \cdot (2),$

式中：m是x在Curvelet域中稀疏系数，D＝AR，R是x的Curvelet逆变换，即x＝Rm，λ是正则化参数；In the formula: m is the sparse coefficient of x in the Curvelet domain, D=AR, R is the Curvelet inverse transformation of x, that is, x=Rm, and λ is the regularization parameter;

F、令目标函数j最小，利用FISTA算法求解x；F. Make the objective function j the smallest, and use the FISTA algorithm to solve x;

G、输出滤波后的地震数据。G. Outputting the filtered seismic data.

有益效果：建立一个不同子波长度的楔形模型，对无噪声和对含高斯白噪声的地震数据使用FISTA算法进行Curvelet域精细尺度下地震数据匹配,经实际对比试验，匹配后的结果波形一致性高，并且噪声得到了好的压制，明显的提高了地震数据处理的分辨率、信噪比。通过实验结果对比可以看出，匹配滤波效果好，可以有效地压制噪声，提高了地震数据信噪比、分辨率，增强了地震勘探效果。Beneficial effect: establish a wedge-shaped model with different wavelet lengths, use the FISTA algorithm to match the seismic data at the fine scale in the Curvelet domain for the seismic data without noise and Gaussian white noise, and after the actual comparison test, the waveform consistency of the matching results High, and the noise has been well suppressed, significantly improving the resolution and signal-to-noise ratio of seismic data processing. Through the comparison of experimental results, it can be seen that the matched filter has a good effect, can effectively suppress noise, improve the signal-to-noise ratio and resolution of seismic data, and enhance the seismic exploration effect.

附图说明：Description of drawings:

图1为基于Curvelet变换的地震信号去噪流程图。Figure 1 is a flow chart of seismic signal denoising based on Curvelet transform.

图2(a)为模拟数据1波形图，由最小相位子波合成的褶积模型构成。Fig. 2(a) is the waveform diagram of the simulation data 1, which is composed of the convolution model synthesized by the minimum phase wavelet.

图2(b)是对模拟数据1进行频谱分析得到的振幅谱。Figure 2(b) is the amplitude spectrum obtained by performing spectrum analysis on the simulation data 1.

图2(c)是对模拟数据1进行频谱分析得到的相位谱。Figure 2(c) is the phase spectrum obtained by performing spectrum analysis on the analog data 1.

图3(a)为模拟数据2波形图，由零相位子波合成的褶积模型构成。Figure 3(a) is the waveform diagram of the simulation data 2, which is composed of a convolution model synthesized by zero-phase wavelets.

图3(b)是对模拟数据2进行频谱分析得到的振幅谱。Figure 3(b) is the amplitude spectrum obtained by performing spectrum analysis on the analog data 2.

图3(c)是对模拟数据2进行频谱分析得到的相位谱。Figure 3(c) is the phase spectrum obtained by performing spectrum analysis on the analog data 2.

图4(a)为模拟数据1在Curvelet域经过FISTA计算匹配得到的地震波形图。Fig. 4(a) is the seismic waveform obtained by FISTA calculation and matching of simulated data 1 in the Curvelet domain.

图4(b)是对图4(a)波形图进行频谱分析得到振幅谱图。Fig. 4(b) is an amplitude spectrogram obtained by performing spectrum analysis on the waveform diagram in Fig. 4(a).

图4(c)是对图4(a)波形图进行频谱分析得到相位谱图。Fig. 4(c) is a phase spectrogram obtained by performing spectrum analysis on the waveform diagram in Fig. 4(a).

图5(a)为模拟数据3波形图，由模拟数据1加入一定量的高斯白噪声构成，信噪比为-7.45dB。Fig. 5(a) is a waveform diagram of analog data 3, which is formed by adding a certain amount of Gaussian white noise to analog data 1, and the signal-to-noise ratio is -7.45dB.

图5(b)为模拟数据3进行频谱分析得到的振幅谱图。Fig. 5(b) is the amplitude spectrogram obtained by the spectrum analysis of the simulated data 3.

图5(c)为模拟数据3进行频谱分析得到的相位谱图。Fig. 5(c) is the phase spectrogram obtained by performing spectrum analysis on the simulated data 3.

图6(a)为模拟数据4波形图，由模拟数据2加入一定量的高斯白噪声构成，信噪比为1.09dB。Fig. 6(a) is a waveform diagram of analog data 4, which is formed by adding a certain amount of Gaussian white noise to analog data 2, and the signal-to-noise ratio is 1.09dB.

图6(b)为模拟数据4进行频谱分析得到的振幅谱图。Fig. 6(b) is the amplitude spectrogram obtained by performing spectrum analysis on the simulated data 4.

图6(c)为模拟数据4进行频谱分析得到的相位谱图。Fig. 6(c) is the phase spectrogram obtained by performing spectrum analysis on the simulated data 4.

图7(a)是对加入噪声的模拟数据3进行Curvelet域匹配滤波的结果。Fig. 7(a) is the result of matching filtering in Curvelet domain on the analog data 3 with noise added.

图7(b)是对图7(a)波形图进行频谱分析得到振幅谱图。Fig. 7(b) is an amplitude spectrogram obtained by performing spectrum analysis on the waveform diagram in Fig. 7(a).

图7(c)是对图7(a)波形图进行频谱分析得到相位谱图。Fig. 7(c) is a phase spectrogram obtained by performing spectrum analysis on the waveform diagram in Fig. 7(a).

具体实施方式：Detailed ways:

下面结合附图和实施例对本发明作进一步详细的说明：Below in conjunction with accompanying drawing and embodiment the present invention will be described in further detail:

如图1是基于Curvelet变换的地震信号去噪流程图，Figure 1 is a flow chart of seismic signal denoising based on Curvelet transform.

A、输入含噪声信号道集；A. Input the noise-containing signal gather;

C、利用匹配算子A构建等式Ax＝b+n………(1)，C. Utilize the matching operator A to construct the equation Ax=b+n...(1),

E、将等式转化为目标函数： $J = {| | Dm - b | |}_{2}^{2} + λ {| | m | |}_{1} \cdot \cdot \cdot (2),$ E. Transform the equation into an objective function: $J = {| | Dm - b | |}_{2}^{2} + λ {| | m | |}_{1} &Center Dot; \cdot \cdot (2),$

G、输出滤波后的地震数据。G. Outputting the filtered seismic data.

通过Curvelet变换，由原始地震数据得到匹配算子A。Through Curvelet transformation, the matching operator A is obtained from the original seismic data.

再利用匹配算子A构建等式Ax＝b+n，在匹配滤波过程中：A匹配滤波算子矩阵，x为待匹配数据，b为匹配道数据；在压制随机噪声过程中：A表示为产生随机噪声的算子矩阵，x是不含随机噪声的数据，b为伴随随机噪声的地震记录；以上两种情况下。n都是未知的噪声数据；Then use the matching operator A to construct the equation Ax=b+n, in the matched filtering process: A matched filtering operator matrix, x is the data to be matched, b is the matching channel data; in the process of suppressing random noise: A is expressed as Operator matrix for generating random noise, x is the data without random noise, b is the seismic record with random noise; in the above two cases. n is unknown noise data;

对等式(1)中数据进行Curvelet变换，转换到Curvelet域变为稀疏数据；Carry out Curvelet transformation to the data in equation (1), convert to Curvelet domain and become sparse data;

离散Curvelet变换原理如下：The principle of discrete Curvelet transform is as follows:

定义一对窗函数：径向窗函数W(r)和角度窗函数V(t)，其中，r∈(1/2,2),t∈[-1,1]。它们都满足可允许条件：Define a pair of window functions: radial window function W(r) and angular window function V(t), where r∈(1/2,2), t∈[-1,1]. They all satisfy the allowable conditions:

${Σ Σ}_{j j = = - - \infty \infty}^{\infty \infty} {W W}^{22} ((22^{j j} r r)) = = 11,, r r &Element; &Element; ((33 / / 4,3 4,3 / / 22)) \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; ((33))$

${Σ Σ}_{j j = = - - \infty \infty}^{\infty \infty} {V V}^{22} ((t t - - l l)) = = 11,, t t &Element; &Element; ((- - 11 / / 2,1 2,1 / / 22)) \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; ((44))$

对于每一个j≥j₀，在频域中定义频窗U_j如下：For each j≥j ₀ , the frequency window U _j is defined in the frequency domain as follows:

${U u}_{j j} ((r r,, θ θ)) = = 22^{- - 33 j j / / 44} W W ((22^{j j} r r)) V V ((\frac{22^{[[j j / / 22]]} θ θ}{22 π π})) \cdot \cdot \cdot \cdot \cdot \cdot ((55))$

其中[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²。where [j/2] is the integer part of j/2. The support interval of U _j is restricted by the support intervals of W and V to obtain a wedge-shaped area, and the wedge-shaped area conforms to the characteristics of anisotropic scale. Assuming that the frequency expression of atom φ _j is φ _j (ω)=U _j (ω), if φ _j on scale j is known, then Curvelets on other 2 ^-j scales can be obtained by φ _j through rotation and translation . Suppose the rotation angle sequence θ _l = 2π·2 ^-[j/2 ]·l at equal intervals, l is a non-negative integer, 0≤θ≤2π; the translation parameter is: k=(k ₁ ,k ₂ )∈Z ² .

综合以上概念，定义尺度在2^-j，方向θ_l，平移参数(k₁,k₂)处的Curvelet为：Combining the above concepts, define the Curvelet at the scale 2 ^-j , the direction θ _l , and the translation parameters (k ₁ ,k ₂ ) as:

${φ φ}_{j j,, l l,, k k} ((x x)) = = {φ φ}_{j j} (({R R}_{{θ θ}_{l l}} ((x x - - {x x}_{k k}^{((j j,, k k))})))) \cdot \cdot \cdot &Center Dot; \cdot \cdot ((66))$

其中 $x_{k}^{(j, l)} = R_{θl}^{- 1} (k_{1} \cdot 2^{- j}, k_{2} \cdot 2^{- j / 2}),$ R_θ和θ_l旋转获得。Curvelet变换可以表示为：in $x_{k}^{(j, l)} = R_{θl}^{- 1} (k_{1} &Center Dot; 2^{- j}, k_{2} &Center Dot; 2^{- j / 2}),$ R _θ and θ _l are obtained by rotation. Curvelet transformation can be expressed as:

$c c ((j j,, l l,, k k)) = = < < f f,, {φ φ}_{j j,, l l,, k k} > > = = {&Integral; &Integral;}_{{R R}^{22}} f f ((x x)) \overset{&OverBar; &OverBar;}{{φ φ}_{j j,, l l,, k k} ((x x))} dx dx \cdot &Center Dot; \cdot \cdot \cdot &Center Dot; ((77))$

通过笛卡尔坐标系下的f[t₁,t₂],0≤t₁,t₂＜n为信号的输入，则Curvelet变换的离散形式为：With f[t ₁ ,t ₂ ], 0≤t ₁ , t ₂ <n in the Cartesian coordinate system as the signal input, the discrete form of the Curvelet transform is:

$c c [[j j,, l l,, k k]] = = \underset{00 \leq \leq {t t}_{11},, {t t}_{22} < < n no}{Σ Σ} f f [[{t t}_{11},, {t t}_{22}]] \overset{&OverBar; &OverBar;}{{φ φ}_{j j,, l l,, k k} [[{t t}_{11},, {t t}_{22}]]} \cdot \cdot \cdot \cdot \cdot \cdot ((88))$

其中φ_j,l,k[t₁,t₂]为离散Curvelet波形。Among them φ _{j, l, k} [t ₁ , t ₂ ] are discrete Curvelet waveforms.

通常情况下受勘探条件以及噪声的影响，会使得线性问题求解存在多解性，需要引入正则化方法对解进行约束。从信息论的角度来看，最优匹配处理时一个欠定问题，Chen(2001)等则在理论上证明了在一定条件下，用稀疏建模具有最佳的重构效果，因此，为降低多解性，采用稀疏约束反演法Under normal circumstances, due to the influence of exploration conditions and noise, there will be multiple solutions in the solution of linear problems, and regularization methods need to be introduced to constrain the solutions. From the perspective of information theory, the optimal matching process is an underdetermined problem. Chen (2001) et al. proved in theory that under certain conditions, sparse modeling has the best reconstruction effect. Therefore, in order to reduce multiple Solvability, using the sparse constraint inversion method

将数据转换到Curvelet域变为稀疏数据后，将方程式(1)转化为下式，只需将目标函数J最小，便可求得解：After transforming the data into the Curvelet domain and turning it into sparse data, the equation (1) is converted into the following formula, and the solution can be obtained only by minimizing the objective function J:

$J J = = {| | | | Dm Dm - - b b | | | |}_{22}^{22} + + λ λ {| | | | m m | | | |}_{11} \cdot \cdot \cdot &Center Dot; \cdot \cdot ((99))$

其中，m是公式(1)中x在Curvelet域中的稀疏系数，D＝AR，R是x的Curvelet逆变换，即，x＝Rm。λ是正则化参数，其控制L1范数的约束项与L2范数的误差项之间的权重。方程式(9)是curvelet域中待求解的稀疏系数。由于方程式(9)的D可以用AR表示，可以将方程式写成如下式子：Wherein, m is the sparse coefficient of x in the Curvelet domain in formula (1), D=AR, and R is the inverse Curvelet transform of x, ie, x=Rm. λ is a regularization parameter that controls the weight between the constraint term of the L1 norm and the error term of the L2 norm. Equation (9) is the sparse coefficients to be solved in the curvelet domain. Since D in equation (9) can be expressed by AR, the equation can be written as follows:

$J J = = {| | | | ARm ARm - - b b | | | |}_{22}^{22} + + λ λ {| | | | m m | | | |}_{11} \cdot &Center Dot; \cdot \cdot \cdot &Center Dot; ((1010))$

令目标函数j最小，利用FISTA算法对方程式(2)进行迭代运算求解，得到待匹配数据x。Make the objective function j the smallest, use the FISTA algorithm to iteratively solve the equation (2), and obtain the data x to be matched.

采用FISTA(Fast Iterative Shrinkage-Thresholding Algorithm)(Berkand Teboulle,2009)求解方程式(10)。FISTA作为一种有效的解决线性反演问题的算法，具有收敛速度快，迭代次数少的优点，大大提高计算效率。FISTA求解过程如下：Use FISTA (Fast Iterative Shrinkage-Thresholding Algorithm) (Berk and Teboulle, 2009) to solve equation (10). As an effective algorithm for solving linear inversion problems, FISTA has the advantages of fast convergence speed and fewer iterations, which greatly improves computational efficiency. The FISTA solution process is as follows:

第1步：b₁＝x₀和t₁＝0Step 1: b ₁ =x ₀ and t ₁ =0

第k步：step k:

$\begin{matrix} {x x}_{k k} = = soft soft (({x x}_{k k - - 11} + + \frac{11}{α α} {A A}^{H h} ((b b - - {Ax Ax}_{k k - - 11})),, \frac{11}{α α})),, \\ {t t}_{k k + + 11} = = \frac{11 + + \sqrt{11 + + 44 {t t}_{k k}^{22}}}{22},, \\ {b b}_{k k + + 11} = = {x x}_{k k} + + \frac{{t t}_{k k - - 11}}{{t t}_{k k + + 11}} (({x x}_{k k} - - {x x}_{k k - - 11})) . . \end{matrix} \cdot &Center Dot; \cdot \cdot \cdot &Center Dot; ((1111))$

其中，soft(x,t)定义为软阈值，当输入x为复数时，令x＝ze^iω。软阈值可表示为：Among them, soft(x,t) is defined as a soft threshold, and when the input x is a complex number, set x=ze ^iω . The soft threshold can be expressed as:

$soft soft ((x x,, t t)) = = \{\begin{matrix} ((z z - - t t)) {e e}^{iω iω},, & if if & z z > > t t \\ 00,, & if if & z z \leq \leq t t \end{matrix} \cdot \cdot \cdot \cdot \cdot &Center Dot; ((1212))$

FISTA算法具有实现简单，收敛速度快，并能将复数的运算也包括在内，扩大了使用范围。The FISTA algorithm has the advantages of simple implementation, fast convergence speed, and can also include complex number operations, which expands the scope of use.

不同种类地震子波合成的楔形模型代表不同震源或者年代的地震记录。子波波长较长，时间分辨率较低，反之亦然。不同子波分辨薄层和楔形的能力也不同，文中给出两种不同类型子波的褶积模型。为了更好地模拟真实情况以证明本方法的有效性，本文建立一个不同子波长度的楔形模型，且模型中存在一个薄层。本发明中选用时间分辨率低的地震记录作为匹配道，时间分辨率高的地震记录作为参考道。The wedge-shaped model synthesized by different types of seismic wavelets represents seismic records of different sources or ages. Wavelets with longer wavelengths have lower temporal resolution, and vice versa. Different wavelets have different ability to distinguish thin layer and wedge shape, and the convolution models of two different types of wavelets are given in this paper. In order to better simulate the real situation and prove the effectiveness of this method, a wedge model with different wavelet lengths is established in this paper, and there is a thin layer in the model. In the present invention, seismic records with low temporal resolution are selected as matching traces, and seismic records with high temporal resolution are selected as reference traces.

模型共有45道，采样时间2ms，零相位和最小相位子波长度均为140ms，模型厚度是4000ms，薄层厚度是120ms。对模型进行Curvelet域多尺度分解，而只分解了三个尺度的Curvelet变换。The model has a total of 45 channels, the sampling time is 2ms, the length of the zero phase and the minimum phase wavelet are both 140ms, the model thickness is 4000ms, and the thickness of the thin layer is 120ms. The multi-scale decomposition of the Curvelet domain is performed on the model, and only the three-scale Curvelet transformation is decomposed.

如图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域最优约束范数匹配具有很好的效果。Figure 2 shows the convolution model of the minimum phase wavelet synthesis, and Figure 2(b) and Figure 2(c) respectively show the amplitude spectrum and phase spectrum obtained from the spectrum analysis of the convolution model of the minimum phase seismic wavelet. Range 10-80Hz. Due to the influence of the wavelet type, the ability to distinguish wedges and thin layers in the convolution model is greatly reduced, which brings inconvenience to subsequent processing and interpretation. Therefore, the convolution model obtained by the minimum phase wavelet is used as the matching channel. Figure 3(a) is the waveform diagram of the simulation data 2, which is composed of a convolution model synthesized by zero-phase wavelets. The amplitude spectrogram 3(b) and the phase spectrogram 3(c) are obtained by analyzing the spectrum of the simulated data 2. It can be seen from the spectrum analysis that the main frequency energy of the seismic wave is also concentrated around 10-80Hz, and the phase spectrum 3(c) is a nonlinear spectral line, and there is a big difference in phase compared with the phase spectrum 2(c). The simulated data 2 has a higher time resolution, which facilitates the subsequent processing and interpretation of seismic data, so the simulated data 2 is used as the reference seismic data. Figure 4(a) is the seismic waveform obtained by FISTA calculation and matching in the Curvelet domain; the amplitude spectrum 4(b) and phase spectrum 4(c) are obtained by performing spectrum analysis on the matched filtering results. The waveform, amplitude spectrum, and phase spectrum of the reference trace seismic data 2 are compared and analyzed. After matching, the waveform, amplitude, and phase of the seismic data are highly consistent, indicating that the optimal constraint norm matching in the Curvelet domain has a good effect.

在实际地震勘探中，地震数据往往存在噪声。为更好的模拟实际情况，本发明对模拟数据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。对比振幅谱和相位谱,它们也具很较高的一致性。In actual seismic exploration, seismic data often have noise. In order to better simulate the actual situation, the present invention adds a certain amount of Gaussian white noise to the simulated data 1 and simulated data 2. Adding noise to analog data 1 is called analog data 3. Figure 5(a) is the waveform of analog data 3, and the signal-to-noise ratio is -7.45dB. Perform spectrum analysis on the simulated data 3 to obtain the amplitude spectrogram 5(b) and the phase spectrogram 5(c). Add Gaussian white noise to the simulation data 2 to obtain the simulation data 4, the signal-to-noise ratio is 1.09dB, as shown in Fig. 6(a). Spectrum analysis is performed on the simulated data with noise added to obtain the amplitude spectrogram 6(b) and the phase spectrogram 6(c). Due to the high signal-to-noise ratio and zero-phase wavelet type, it has a strong ability to distinguish wedges and thin layers, so the simulated data 4 is used as the reference trace seismic data. Fig. 7(a) is the result of Curvelet domain matching filtering on the noise-added analog data 3, and the amplitude spectrum 7(b) and phase spectrum 7(c) are obtained through spectrum analysis on the matched seismic data. Compared with the simulation data 4, the waveform has a high similarity, and the signal-to-noise ratio is increased to 26dB. Comparing the amplitude spectrum and phase spectrum, they also have a high consistency.

Curvelet变换都具有很高的稀疏性，对地震数据进行多尺度分解，为信号和噪声的分离奠定了基础，因而利用稀疏约束求解匹配问题更加准确。本发明方法取得了更好的效果，为匹配处理提供了新的思路,将最优范数和Curvelet分解结合在一起，为分辨薄层和新老地震数据提供了有价值的信息。FISTA算法具有较快的计算效率，且具收敛性能好等特点，易于计算机编程等特点,能够很好地应用于数据拟合的测量实践中。Curvelet域匹配滤波与噪声压制的结果具有同相轴连续性好、波形、振幅谱和相位谱一致性高的特点。The Curvelet transform has high sparsity, and multi-scale decomposition of seismic data lays the foundation for the separation of signal and noise, so it is more accurate to use sparse constraints to solve the matching problem. The method of the invention achieves better results, provides a new idea for matching processing, combines the optimal norm and Curvelet decomposition, and provides valuable information for distinguishing thin layers and new and old seismic data. The FISTA algorithm has the characteristics of fast calculation efficiency, good convergence performance, and easy computer programming. It can be well applied in the measurement practice of data fitting. The result of matched filtering and noise suppression in Curvelet domain has the characteristics of good event continuity and high consistency of waveform, amplitude spectrum and phase spectrum.

Claims

1. A seismic signal denoising method based on Curvelet transform, characterized in that a wedge-shaped model with different wavelet lengths is established, and there is only one thin layer in the model, and the ability of different wavelets to distinguish thin layers and wedges is also different. Seismic records with low temporal resolution are used as matching traces, and seismic records with high temporal resolution are used as reference traces. Matching trace seismic records and seismic records mixed with Gaussian white noise are processed by Curvelet transform matching filtering method, and compared with reference traces. High-resolution seismic records for comparison.

2. according to the seismic signal denoising method based on Curvelet transform according to claim 1, it is characterized in that, comprising the following steps:

A. Input the noise-containing signal gather;

B. Using Curvelet transform, the matching operator A is obtained from the original seismic data;

C. Using the matching operator A to construct the equation Ax=b+n...(1),

In the matched filtering process: A is the matching operator, x is the data to be matched, b is the matching channel data, and n is random noise;

In the process of suppressing random noise: A represents the operator matrix that generates random noise, x is the data without random noise, and b is the seismic record with random noise; n in the above two cases is unknown noise data;

D, the data in the equation (1) is converted to the Curvelet domain and becomes sparse data;

E. Transform the equation into an objective function:

J = {| | Dm - b | |}_{2}^{2} + λ {| | m | |}_{1} . . . (2),

In the formula: m is the sparse coefficient of x in the Curvelet domain, D=AR, R is the Curvelet inverse transformation of x, that is, x=Rm, and λ is the regularization parameter;

F. Make the objective function j the smallest, and use the FISTA algorithm to solve x;

G. Output the filtered seismic data.