- Email: [email protected]
Desert seismic noise suppression based on an improved low-rank matrix approximation method
Journal Pre-proof Desert seismic noise suppression based on an improved low-rank matrix approximation method
Juan Li, Wei Fan, Yue Li, Baojun Yang, Changgang Lu PII:
S0926-9851(19)30480-X
DOI:
https://doi.org/10.1016/j.jappgeo.2019.103926
Reference:
APPGEO 103926
To appear in:
Journal of Applied Geophysics
Received date:
25 May 2019
Revised date:
24 December 2019
Accepted date:
25 December 2019
Please cite this article as: J. Li, W. Fan, Y. Li, et al., Desert seismic noise suppression based on an improved low-rank matrix approximation method, Journal of Applied Geophysics(2019), https://doi.org/10.1016/j.jappgeo.2019.103926
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2019 Published by Elsevier.
Journal Pre-proof
Desert seismic noise suppression based on an improved low-rank matrix approximation method Juan Li,1 Wei Fan,1 Yue Li,1* Baojun Yang2 and Changgang Lu3 1
Department of Information, College of Communication Engineering, Jilin University, Changchun 130012, China. E-mail: [email protected] Department of Geophysics, Jilin University, Changchun 130026, China
3
College of Automotive Engineering, Jilin University, Changchun, China
of
2
ro
Abstract
-p
The suppression of random noise is a crucial step before seismic data analysis. Random noise in
re
desert areas has the characteristics of low frequency and non-stationary, and there is serious spectrum
lP
aliasing between random noise and effective signals, which makes it difficult to suppress such noise. In recent years, some methods based on signal rank minimization have achieved remarkable results
na
in seismic random noise suppression. Since the implementation of low rank matrix approximation is
Jo ur
an iterative process, noise estimation is an indispensable step before each iteration, but also an important step. The noise estimation method previously used is to calculate the residuals of the original noisy patch data and the corresponding iterative denoising version, which is intuitively considered as the filtered noise. This method may be very inaccurate in the case of high noise levels or complex seismic records. In this paper, a noise estimation method based on geometric texture is introduced to estimate the noise level by selecting weak textured patches in all seismic texture patches. At the same time, we reduce the loss of effective signals by truncating the singular values in each iteration. Experiments on both synthetic and field seismic data show that this method has better effect on suppressing random noise in desert areas.
1
Journal Pre-proof Key Words: Desert seismic noise; low rank matrix approximation; noise estimation; noise suppression; truncated singular value;
1 INTRODUCTION In the process of collecting seismic signals, it is inevitable that it will be affected by various types
of
of noise, which will have strong interference on the subsequent analysis of seismic data. If the
ro
random noise is suppressed, the subsequent analysis of seismic data will be more accurate (Zhang and Klemperer, 2005, 2010; Souza et al., 2016; Yuan et al., 2018; Li et al., 2017; Chen et al., 2016;
-p
Wu, 2017). Therefore, noise removal plays an important role in increasing the quality of seismic
re
exploration data. The frequency of random noise in desert areas is relatively low, and the frequency
lP
bands of effective signals and random noise are seriously overlapped, which makes some methods
na
suitable for suppressing white noise and separating signals and noise in frequency domain ineffective
Jo ur
on suppressing the random noise in desert areas.
Previously, abundant methods have been proposed and applied to suppress random noise in real seismic data. Sparse transform-based methods, such as Fourier transform (Sacchi et al., 1998; Zhai, 2014), wavelet transform (Gaci, 2014; Chen and Song, 2018; Shan et al., 2009), curvelet transform (Tang and Ma, 2010; Herrmann and Hennenfent, 2008; Gorszczyk et al., 2014; Liu et al., 2018; Lari and Gholami, 2014), shearlet transform ( Tang et al., 2018), seislet transform (Bai and Wu, 2019; Dalai et al., 2019), Radon transform (Xue et at., 2017) transform the seismic data into the sparse domain to obtain the corresponding sparse coefficients, then, the thresholding process is performed on the sparse coefficients, and the processed sparse coefficients are transformed back into the space-time domain to achieve the purpose of suppressing noise. Mode decomposition-based methods, 2
Journal Pre-proof
such as empirical mode decomposition (EMD), variational mode decomposition (VMD) and their improvements (Bekara and van der Baan, 2009; Kopsinis and McLaughlin, 2009; Liu et al., 2014; Chen, 2016; Yu and Ma, 2018; Liu et al., 2017; Zhou and Zhu, 2019), mainly decompose noisy seismic data into many different components, then, the components containing effective signals are selected and the components of random noise are discarded. Prediction-based methods , such as f-x
of
deconvolution and t-x predictive filtering (Gulunay, 1986; Abma and Claerbout, 1995), mainly utilize the predictability of effective signals to construct predictive filters to suppress noise and
ro
enhance signals. Some methods that use non-local self-similarity of seismic records are also used to
-p
process seismic records(Shao et al., 2019; Wang et al., 2019; Amani et al. 2017) . Some methods can
re
recover the signal effectively by use the low rank feature of seismic signal in time-frequency domain,
lP
for example, Rasoul Anvari et al. (2019) consider the ideal seismic signal in the transform domain as
na
both sparse and low rank. They first compute the short-time Fourier transform of the noisy seismic signal, then use a nonconvex penalty function to estimate the sparse low-rank matrix, and finally use
Jo ur
the estimated sparse low-rank matrix to synthesize the seismic signal. Rasoul Anvari et al. (2017) first use the synchrosqueezed wavelet transform to transform the seismic data trace by trace into the sparse domain to obtain the sparse time-frequency representation, and then use the Optshrink algorithm to decompose the sparse time-frequency representation into semilow-rank and sparse components, finally back-transform the semilow-rank component to the time domain to recover the denoised seismic trace using inverse synchrosqueezed wavelet transform. The above existing methods have played a certain role in dealing with large sets of oil and gas layers, random noise reduction in medium and shallow exploration, and improving the signal-to-noise ratio of seismic records, and obtained good results in seismic exploration engineering. 3
Journal Pre-proof
In recent years, a large number of scholars have conducted in-depth research on the low-rank method of signals. It can be considered that seismic signals are of low rank after some methods of signal rearrangement steps. For example, Rasoul Anvari et al. (2019) model the 3-D seismic data as a tensor data of size n1 n2 n3 , which can be expressed as a combination of a low rank component and a sparse component, then compute the Fourier transform of the tensor along the third dimension
of
of the tensor, and shrink the singular values in the Fourier domain to extract the low rank component. They repeat the above steps until the Frobenius norm of the error matrix reaches the desired value,
ro
the matrix constructed by stacking non-local similar patches from the noisy seismic signals is of low
-p
rank in the low-dimensional subspace of the given high-dimensional space (Li et al., 2017). EMD
re
can empirically decompose a multi-dip seismic image that is not of low rank into multiple single-dip
lP
seismic images that are low-rank individually (Chen et al., 2017). The Hankel matrix constructed
na
from seismic signals will have a lower rank (Chen and Sacchi, 2015). The denoising mechanism of these methods is also different. Li et al. (2017) takes advantage of the prior knowledge that the
Jo ur
singular values have clear physical meanings and should be treated differently, and assigns different weights to singular values of different amplitudes. Chen et al. (2017) constructs a hankel matrix for each empirical component after EMD decomposition, and then apply SVD to each hankel matrix, then retain the largest singular value of the same number as the rank of modal component and set others to zero. Chen and Sacchi (2015) decomposes the Hankel matrix constructed by seismic signals into the product of two lower dimensional factor matrices, then the two matrices are solved by weighted cost function.
Noise estimation is a very important step of the low rank matrix approximation denoising framework, the estimated noise standard deviation will be used to calculate the singular value of the 4
Journal Pre-proof
pure signal and the corresponding weight thresholds, which are decisive parameters for denoising, thus the accuracy of the noise estimation is directly related to the quality of the denoising results. The noise estimation method in (Gu et al., 2014) is to calculate the difference between the stacked similar patch matrices of noisy seismic records and the corresponding matrices of current iterative denoising results, which is intuitively assumed to be the noise. The denoising results of this method are greatly
of
influenced by the initial value of noise level, and the results of current iteration will affect the results of next iteration. If the initial value of noise level is not set accurately, or it is not only noise that is
ro
filtered out in a certain iteration, then the denoising results will not be very satisfactory. Especially in
-p
the field seismic data processing, because of the complexity of the field seismic data, it would be
re
inaccurate to use the original noise estimation method.
lP
To make up for this deficiency, we introduce the noise estimation method based on geometric
na
texture (Liu et al., 2012) into the denoising framework of low rank matrix approximation in this
Jo ur
paper. The null hypothesis test is used to select weak textured patches from all seismic texture patches and then principal component analysis (PCA) is performed on these weak textured patches to estimate the noise level. With this method, the denoising results no longer depend on the initial value of the noise level, and the noise estimations before each iteration are no longer dependent on the results of the previous iteration, which are independent of each other, so as to get the estimates closer to the real value. Meanwhile, taking into account the physical meaning of the singular value, the singular value of the effective signal is related to its energy. Generally, the singular value of the signal with higher energy is larger than that of the signal with lower energy. In the original framework, the smaller singular value will become zero after the weight threshold processing, which will result in the loss of low energy signal. So before each iteration, we divide the singular value 5
Journal Pre-proof
matrix into two matrices with the same dimension as the original one, one retaining larger singular value and the other retaining smaller singular value. Two seismic signals can be obtained from these two matrices, one representing high energy signals and the other representing low energy signals. We add this difference between the two seismic signals to the iterative regularization term to reduce the loss of effective signals. Tests with synthetic and field seismic data prove that the method proposed
Jo ur
na
lP
re
-p
ro
of
in this paper has better effect on suppressing seismic random noise in desert areas.
Fig. 1: Denoising records and corresponding filtered records after the 2nd, 4th and 6th iteration of WNNM. (a) Noisy record. (b),(c) the denoising records and corresponding filtered records after the 2nd iteration. (d),(e) the denoising records and corresponding filtered records after the 4th iteration. (f),(g) the denoising records and corresponding filtered records after the 6th iteration. 2 Low Rank Matrix Approximation for Seismic Noise Suppression in Desert Areas
6
Journal Pre-proof
2.1 Theory Noisy seismic data in desert areas can be expressed as follows:
Y X E,
(1)
where Y is the noisy seismic datas, X is the potential clean signals, and E is the random noise with a standard deviation of . Our purpose is to recover X from Y , In the case of complex
of
seismic signals, X itself is not low rank. Using the non-local self-similarity of seismic signals, we
ro
select a main patch PYj in a search window of a certain size, and then select patches similar to the
-p
main patch to stack them into a matrix M Yj . At the same time, the corresponding pure signal matrix
re
M X j and random noise matrix M E j can be found to represent (1) as follows：
lP
M Yj M X j M E j .
(2)
na
After this data rearrangement, it can be considered that the matrix M X j is of low rank, so the
Jo ur
method of low rank matrix approximation can be used to recover M X j from M Yj to achieve the purpose of denoising.
The most primitive method for solving (2) is defined as
ˆ X arg min rank (M X ) ,s.t. M Y M X M E . M j j j j j MX j
(3)
But the objective function in (3) is non-convex, and solving it is an NP-hard problem. However, it can be converted to nuclear norm minimization (NNM) by convex relaxation as follows (Gu et al., 2014; Cai et al., 2010):
7
Journal Pre-proof ˆ X arg min 1 || M Y M X || 2F || M X ||* , M j j j j 2 MXj
(4)
where is a fixed positive constant, ||||F is Frobenius norm and ||M X j||* is the nuclear norm of matrix M X j , || ||* i |λi ()| , i () represents the ith singular value of the matrix. The optimal solution of (4) can be obtained by singular value decomposition (Cai et al., 2010), i.e.,
M Yj UVT，ˆ ii S (ii ) max( ii ,0) ,
of
(5)
(6)
ro
ˆ X Uˆ VT , M j
-p
where ii represents the diagonal element of the matrix , ii i (M Yj ) , It can be observed that
re
NNM performs the same threshold for each singular value, which does not take into account the
lP
practical physical meaning of the singular value. The effective information of the seismic signal is mainly represented by a large singular value, so the larger singular value should be assigned a
na
smaller threshold than the smaller singular value. Gu et al. (2017) proposed weighted nuclear norm
which is defined as
Jo ur
minimization (WNNM), assigning different weight thresholds according to different singular values,
ˆ X arg min 1 || M Y M X || 2F || M X || w,* , M j j j j 2 MXj
(7)
where || ||w,* i |wi λi ()| , w w1, w2 , , wn , wi is the non-negative weight value corresponding to the singular value i () , The relationship between wi and i (M X j ) is as follows:
wi
c m , i (M X j )
(8)
8
Journal Pre-proof where c 0 is a constant and m is the number of patches constituting the matrix M Yj . To avoid the denominator being zero, we set 1016 . It can be seen from (8) that the larger the singular value is, the smaller the corresponding wi is. Although (7) is non-convex when the weight vector
w is in a non-descending order, Gu et al. (2017) prove that there exists a local minimum point which is its optimal solution. The solution of (7) is similar to that of (4). By replacing in (5) with
of
wi , we can get the following expression (9)
ro
ˆ ii S (ii ) max( ii wi ,0) .
-p
The other steps are the same as those of (4).
re
There is a problem in obtaining wi , that is, the singular value of the pure seismic signal cannot be
lP
known in advance. But it can be estimated by the singular value of noisy seismic signal and the
na
variance of noise in noisy seismic signal, which is defined as
i (M X j ) max( i2 (M Yj ) m 2 ,0) ,
Jo ur
(10)
where i (M Y j ) can be obtained by calculation, so the question now is how to get . Our work is also carried out around . From (8), (9), (10), we can see the accuracy of will directly affect the accuracy of i (M X j ) , and then it will affect the value of wi , and ultimately affect the effect of noise suppression.
Algorithm 1 Seismic Denoising Using the proposed method Input: Noisy seismic signal Y , parameter , , L , 0 ˆ 0 Y, X ˆ 0high Y, X ˆ low Y 1: Initialize Y0 Y, X
2: for k 1 : L do 3:
9
k ˆ k (Y X ˆ k ) (X ˆ khigh X ˆ low ) iterative regularization Y k X
k
Jo ur
na
lP
re
-p
ro
of
Journal Pre-proof
2.2 Existing noise estimation method in WNNM
In the WNNM framework, the noise standard deviation is an important parameter. Since WNNM is an iterative process, two noise values can be obtained after each iteration: filtered noise and residual noise. The residual noise standard deviation can be estimated from the original noise standard deviation and the standard deviation of the filtered noise. The estimation process is as 10
Journal Pre-proof
follows. First of all, the method cannot obtain the original noise standard deviation , so must be set artificially. It is assumed that the noise standard deviation filtered by WNNM at the k-th iteration can be expressed as kflt , the standard deviation of residual noise in the record to be processed in the k 1 (k+1)-th iteration can be represented by res . Suppose that the result of PYj after the k-th iterative
of
denoising is Pˆ Yk j 1 , because there is not only noise but also signal in Pˆ Yk j 1 , the standard deviation
ro
k 1 cannot be calculated directly in Pˆ Yk j 1 . However, when there is only noise in the record (i.e. res
-p
PYj Pˆ Ykj1 ) filtered by the k-th iteration, the standard deviation kflt can be directly calculated as
re
follows
(11)
lP
kflt || PYj Pˆ Ykj1 || 2 ,
formula
Jo ur
k 1 res 2 ( kflt ) 2 ,
na
then the standard deviation of the residual noise can be estimated by and kflt by the following
(12)
k 1 where res denotes the estimated standard deviation of the residual noise present at (k+1)-th
iteration. The constant 0 is a scaling factor, heuristically introduced to control the re-estimation
of the standard deviation.
There are two main problems in this noise estimation method: First, the original noise standard deviation needs to be set artificially, and its authenticity cannot be guaranteed. Second, if there is signal in the records filtered by k-th iteration, the calculation of the standard deviation of the filtered noise will be inaccurate, which will affect the estimation of the standard deviation of the residual 11
Journal Pre-proof
noise. In order to observe the filtered noise and residual noise obtained by WNNM in each iteration more intuitively, we use it to process the noisy records in Fig. 1(a), and obtain the denoising records and corresponding filtered records after the 2nd, 4th and 6th iteration respectively, as shown in Fig. 1(b)-(g). It can be seen from the experimental results that there are not only filtered noise but also obvious signals in the filtered records of these iterations.
of
The original noise standard deviation and the standard deviation of the filtered noise obtained
ro
at the kth iteration determine the standard deviation of the residual noise in the record to be
-p
processed in the (k+1)-th iteration. That is to say, if the is set inaccurately or there is not only
re
noise in the filtered records as shown in Fig.1, it will affect the denoising results of the (k+1)-th iteration.Therefore, we introduce a noise estimation method based on geometric texture. The noise
lP
standard deviation is directly estimated according to the geometric structure of different texture
na
patches of the noisy seismic signal. It is not necessary to set the initial value of the noise level
Jo ur
manually, and the effect of denoising in the previous iteration does not need to be considered.
3 Noise estimation based on geometric texture Geometric texture-based noise estimation was first proposed and applied to images by Liu et al. (2012). This estimated method has a good practical application effect in images, and the estimated noise standard deviation is very close to the actual value. So we apply this estimation method to noisy seismic signals. In view of the idea proposed by Liu et al. (2012), the noise level can be estimated from weak textured patches in seismic textured patches, which is defined as
ˆ n2 min (cov(W)) ,
(13)
12
Journal Pre-proof where min () is the minimum eigenvalue of the matrix, W represents the selected weak textured patches, and cov() denotes the covariance matrix. Therefore, the key issue in applying geometric texture-based noise estimation method to seismic signals is how to select weak textured patches from seismic texture patches. First, the gradient covariance matrix Cp j of the seismic texture patch can be written as
Dvp j ,
(14)
of
Cp j GpTj Gp j，Gp j Dhp j
ro
where p j denotes the texture patch matrix in matrix M Yj , D h and Dv represent the horizontal
-p
and vertical derivative operators, respectively. The larger eigenvalues of Cp j often reflect more
re
texture information of p j . Theoretically, the noisy flat seismic patch PM j in the reconstructed
lP
seismic signal M Yj can be decomposed into a perfectly flat seismic patch Pf j and a noise patch
PE j , as shown in follow,
na
PM j Pf j PE j .
(15)
Jo ur
It is easy to know that the gradients of the perfectly flat seismic patch are zero, so the expected value of the gradient covariance matrix of noisy flat seismic patches can be calculated as
E (CPM ) E (CPf ) E (CP ) E (CPE ) j
j
Ej
PET DTh Dh PE j E ( Tj T P D D P E j v h E j
j
PET j DTh Dv PE j ） PETj DTv Dv PE j
.
(16)
E (PET DTh Dh PE ) 0 j j 0 E (PET j DTv Dv PE j ) The statistical properties of two diagonal elements are the same. Liu et al. (2012) proved that the distribution of two diagonal elements PETj DT D PE j (where h or v ) can be approximated by gamma distribution. The shape parameter and the scale parameter of the gamma distribution 13
Journal Pre-proof
are as follows:
N 2 2 ， n tr (DT D ) , 2 N
(17)
where N denotes the number of elements in matrix PE j , n is the noise standard deviation to be estimated, and tr (DT D ) denotes the trace of the matrix DT D . In order to select the weak textured patches in the seismic texture patches, we refer to (Liu et al.
of
2012), also using the null hypothesis test method. The null hypothesis is “ the given patch is a flat
ro
patch with random noise”. When the maximum eigenvalue of the gradient covariance matrix of the
re
F 1 (v, , ) ,
lP
valid..The threshold is calculated by
-p
given seismic texture patch is less than a certain threshold, the hypothesis is valid, otherwise it is not
(18)
na
where F 1 denotes the inverse gamma cumulative distribution function and v denotes the
Jo ur
significance level, which must be set manually, generally set to 0.95~0.99. It can be seen from (17) and (18) that the noise standard deviation is also an indispensable parameter for calculating the threshold .
Liu et al (2012) proposes an iterative method to estimate the final noise level. We combine it with seismic signals and summarize it as follows: Firstly, the initial noise standard deviation n( 0) is estimated from the covariance matrix generated by all texture patches in M Yj by (13), and the initial threshold （0） is calculated by (18) according to n( 0) .Weak textured patches W (1) can be selected by null hypothesis test using the threshold （0）, and then n(1) can be estimated by (13). Repeat the above process until the estimated noise level n tends to be stable. We believe that the 14
Journal Pre-proof final n is the closest to the true noise level. To verify the effectiveness of this noise estimation, we add the actual desert noise with different noise levels to the synthetic seismic record in Fig. 2(a), and then use this estimation method to estimate the standard deviation of the noise we added. The true values and estimates of noise standard deviation are listed in Table 1. It can be seen that the noise standard deviation estimated by
of
this noise estimation method is very close to the true values.
1
2
3
4
Estimates
1.0005
1.9975
2.9920
3.9873
5
-p
Truth
ro
Table 1: Comparison of the true values and estimates of noise standard deviation. 15
20
9.9628
14.9432
19.9239
re
4.9830
10
lP
4 Reducing signal loss by truncating singular values
na
Low rank matrix approximation is an iterative denoising process, so it will inevitably lead to the loss of some seismic signals during the iterative process. We consider the physical meaning of the
Jo ur
singular value of the seismic signal. Generally, the singular values of the high-energy effective signals are relatively large, but there are some low-energy signals with small singular values. The singular values of this part may be changed to zero or very close to zero after the threshold processing by (9) in the iterative process, which will result in the loss of this part of the signal. In order to reduce this loss, in each iteration, we use threshold to separate the larger and smaller singular values of the signal as follows: ˆ ˆ high ˆ low ,
(19)
where ˆ high and ˆ low are matrices that retain only larger and smaller singular values in ˆ , respectively. The optimal solution of (7) can be written as 15
Journal Pre-proof ˆ X Uˆ VT U(ˆ high ˆ low )VT M ˆ high ˆ low M Xj M Xj . j
(20)
ˆ high and X ˆ low by aggregating all patch Then, we can obtain two seismic signal records X ˆ high ˆ low ˆ high and X ˆ low can reflect the seismic signal matrices M and M X j . The difference between X Xj
ˆ high and X ˆ low but with different energy levels. We add the difference between contained in both X high and low energy seismic signals to the iteration regularization term to reduce some unnecessary signal loss in the iteration denoising process.
Jo ur
na
lP
re
-p
ro
of
For a clearer description, the improved low rank matrix approximation algorithm is summarized in algorithm 1.
Fig. 2: Denoising results of synthetic seismic record. (a) Ideal signal. (b)Noisy signal. (c) Denoising result of Wavelet Transform. (d) Denoising result of f-x deconvolution. (e) Denoising result of WNNM. (f) Denoising result of the proposed method.
16
ro
of
Journal Pre-proof
re
-p
Fig. 3: Difference between before and after denoising. (a) Added noise. (b) Difference between before and after denoising of Wavelet Transform. (c) Difference between before and after denoising of f-x deconvolution. (d) Difference between before and after denoising of WNNM. (e) Difference between before and after denoising of the proposed method.
lP
5 Application to noise suppression for seismic records in desert areas
na
5.1 Experiments on synthetic seismic data
To verify the effectiveness of the proposed method, we first test it on the synthetic seismic data. The
Jo ur
synthetic seismic model has a total of seven effective events, which are generated by Ricker wavelets with dominant frequencies 30 Hz and 25 Hz, as shown in Fig. 2(a). Then we add real random noise in desert areas to the synthetic seismic data, so that the signal-to-noise ratio of the synthetic record is -5.1402 dB, as shown in Fig. 2(b). When processing this record, we set the size of the search window to 40 40 , and the size of the seismic texture patch to 1313 , this setting is only suitable for processing records at this noise level. When processing synthetic records with different noise levels, the sizes of search window and seismic patch are set differently. When the noise level is low, they are set smaller, and when the noise level is high, they are set larger. At the same time, we set parameters
0.1 , 0.4 , and 0.5 to process synthetic records at all signal-to-noise ratios. 17
Journal Pre-proof
In contrast experiments, Wavelet Transform, f-x deconvolution and WNNM are used to process the same synthetic records to compare with the method proposed in this paper. We adjust the parameters of each method to make the denoising result the best. The denoising results of these three methods are shown in Fig. 2(c), (d), (e). The denoising results of the proposed method are shown in Fig. 2(f). At the same time, to observe their retention of the signal, we obtained the difference
of
between before and after denoising, as shown in Fig. 3. To observe and compare the denoising results of these four methods more comprehensively, we achieve a more detailed knowledge on
ro
signal-preserving performance by analyzing a single trace. For this purpose, we extract the third trace
-p
of synthetic records for comparison, as shown in Fig. 4. Fig. 4(a)-(d) are a single trace denoising
re
result of wavelet transform, f-x deconvolution, WNNM and the proposed method in this paper,
lP
where the black line represents the ideal signal, the blue line represents the noisy signal, and the red
na
line represents the denoised signal. Then we use FK spectrum to analyze the denoising results of these four methods in the frequency domain. Fig. 5(a)-(f) represents the FK spectrum corresponding
Jo ur
to the denoising result in Fig 2(a)-(f), and Fig. 6(a)-(e) represents the FK spectrum corresponding to the difference between before and after denoising in Fig. 3(a)-(e).
Fig. 4: Comparison on the 3th trace of denoised records. (a) Wavelet transform. (b) f-x deconvolution. (c) WNNM. (d) The proposed method.
18
ro
of
Journal Pre-proof
Jo ur
na
lP
re
-p
Fig. 5: The FK spectrum of denoising results. (a) Ideal signal. (b) noisy signal. (c) Denoising result of Wavelet Transform. (d) Denoising result of f-x deconvolution. (e) Denoising result of WNNM. (f) Denoising result of the proposed method.
Fig. 6: The FK spectrum of the differences between before and after denoising. (a) Added noise. (b) Wavelet Transform. (c) f-x deconvolution. (d) WNNM. (e) The proposed method.
Table 2: SNRs of different denoising results. SNR(dB) the original data
Wavelet transform
f-x deconvolution
19
WNNM
The proposed method
Journal Pre-proof
11.1813
12.4269
17.0585
17.3263
7.0604
10.5032
12.1593
16.9041
17.2641
5.0581
9.5991
11.8136
16.4879
17.1718
3.0859
8.4850
11.3862
15.4712
17.0639
1.0935
7.1954
9.9314
11.2022
15.6333
-1.0050
5.5561
8.9776
-3.0694
3.8081
7.8492
-5.1402
1.9518
-7.0256
0.1969
-9.0900
-1.7741
14.0166
8.7879
12.8456
6.5472
10.7549
5.2267
3.3723
7.0238
3.6585
1.2603
4.5704
-p
ro
10.5486
re
of
9.0389
Jo ur
na
lP
7.0543
After the above all-round and multi-angle analysis, it can be seen that the noise suppression of wavelet transform is not complete and the effective signals are not recovered well; The noise suppression of f-x deconvolution is not very thorough, the noise in the same frequency band as the 20
Journal Pre-proof
effective signals is preserved, and some signals are lost; Although the advantage in noise suppression of it is very obvious, WNNM also has a large amount of signal loss; The proposed method suppresses the noise thoroughly, the recovery of the effective signal is the best and the processing result is the closest to the ideal signal.
Next, we use these four methods to process synthetic records of random noise with different levels
of
in desert areas. The signal-to-noise ratio (SNR) and mean square error (MSE) of denoising results are
ro
calculated to quantitatively analyze the effectiveness of the proposed method. SNR and MSE are
-p
calculated as follows:
i 1
N
| Xˆ (i) X(i) | i 1
1 N
N
(Xˆ (i) X(i))2 , i 1
(21)
(22)
na
MSE
), 2
lP
SNR(dB) 10 log10（
| X(i) |2
re
N
MSE
Jo ur
ˆ is the denoised signal. where X is the ideal signal and X
the original SNR
Wavelet transform
Table 3: MSEs of different denoising results
f-x deconvolution
WNNM
The proposed method
9.0389
0.0191
0.0144
0.0049
0.0047
7.0604
0.0224
0.0153
0.0051
0.0047
5.0581
0.0276
0.0166
0.0056
0.0048
3.0859
0.0356
0.0183
0.0071
0.0049
21
Journal Pre-proof
0.0479
0.0255
0.0191
0.0069
-1.0050
0.0699
0.0318
0.0222
0.0100
-3.0694
0.1046
0.0412
0.0332
0.0131
-5.1402
0.1604
0.0495
0.0557
0.0211
-7.0256
0.2402
0.0754
0.1156
0.0499
-9.0900
0.3781
0.1082
of
1.0935
0.0877
re
-p
ro
0.1880
The resulting SNR and MSE for synthetic seismic data are given in Table 2 and 3. The higher the
lP
SNR and the lower the MSE indicates that the denoising result is better. It can be seen from the two
na
tables that the SNR of the proposed method is always the highest and MSE is always the smallest,
Jo ur
which shows that the results of the proposed method are the best.
5.2 Experiments on field seismic data In order to verify the practicability of the proposed method, we use it to process field seismic data in desert areas. We have considered a common shot point seismic gather consisting of 127 traces and 2002 samples per trace for testing, as shown in Fig. 7(a). It can be seen that the reflection events on both sides of the record are almost obscured by random noise and are difficult to distinguish. Similarly, Wavelet transform, f-x deconvolution, WNNM and the proposed method are used to process the record. The results of these four methods are shown in fig. 7(b)-(e). When processing this record, we set the size of search window and seismic texture patch in the proposed method to be 22
Journal Pre-proof
30 30 and 7 7 , respectively, and adjust the results of other methods to the best. Similarly, we also draw the difference between before and after denoising of these methods to analyze their signal retention, as shown in fig. 8. From fig. 7(b)-(e) and Fig. 8, it can be seen that the effect on low-frequency noise suppression in Wavelet transform is not satisfactory, and the effective signals have been distorted. Although f-x deconvolution has a certain effect on the low-frequency
of
noise suppression as a whole, it does not restore the effective signals on both sides of the record.
ro
WNNM is not very good at suppressing low-frequency noise in the deep part of the record. The
-p
proposed method in this paper can suppress the low-frequency noise perfectly, and the effective
Jo ur
na
lP
re
signal on both sides of the record can still be recovered well.
Fig. 7: Results obtained with field seismic data. (a) Field seismic data (common shot point seismic data). (b) Denoising result of Wavelet Transform. (c) Denoising result of f-x deconvolution. (d) Denoising result of WNNM. (e) Denoising result of the proposed method.
23
Journal Pre-proof Fig. 8: Difference between before and after denoising. (a) Difference between before and after denoising of Wavelet Transform. (b) Difference between before and after denoising of f-x deconvolution. (c) Difference between before and after denoising of WNNM. (d) Difference between
ro
of
before and after denoising of the proposed method.
-p
Fig. 9: Partial enlarged views. (a) Original enlarged record. (b) Filtered enlarged record after using Wavelet transform. (c) Filtered enlarged record after using f-x deconvolution. (d) Filtered enlarged
re
record after using WNNM. (e) Filtered enlarged record after using the proposed method.
lP
In order to evaluate the quality of the denoising results of these four methods more clearly, we
na
enlarged the record between trace 77-110 and samples 1000-1900 as shown in Fig. 9. Fig. 9(a)-(e) are partial enlarged views corresponding to Fig. 7(a)-(e) respectively. By comparing the partial
Jo ur
enlargement views, it is more obvious that the proposed method in this paper can suppress the low-frequency noise in desert areas better than the other three methods, and can recover the effective signals well.
Table 3: The computation time of different methods in Fig.2(b) Wavelet transform Time(s)
0.1415
f-x deconvolution 0.0195
WNNM 163.5011
The proposed method 419.2504
6 Dicussion
Our algorithm mainly makes use of the non-local self-similarity of the seismic signal in the 24
Journal Pre-proof
space-time domain. The frequency of the signal does not affect this characteristic. In other words, as long as the stacked similar block matrices satisfy the low-rank property, the value of the dominant frequency of the signal will not have a great impact on the performance of our algorithm. Although our algorithm has a great advantage in suppressing random noise compared with several other methods, but for other types of noise, such as coherent noise, because the noise has a certain
of
correlation, it will cause not only signals are similar, but also noise is similar in the stacked similar block matrix, which will reduce the effectiveness of our method. And our algorithm has a cost in
ro
computation time because it needs to estimate the noise level of each stacked similar block matrix.
-p
We calculated the time spent in processing the record in Fig. 2(b) by the methods mentioned in this
re
paper, which are listed in Table 3. It can be seen that the computation time of our algorithm is much
lP
higher than other methods. In addition, under some geological conditions, the amplitude of the
na
seismic reflection signal will change with the distance between the shot point and the detector, but the presence of noise will make this change blurred and difficult to observe. In this case, our
Jo ur
algorithm can suppress the noise while retaining this change in amplitude, which is more conducive to AVO analysis. So our algorithm can be used as a pre-processing step before AVO analysis, providing a clearer and more effective seismic record for AVO analysis.
7 Conclusion
The method of low rank matrix approximation mainly considers the non-local similarity of seismic signals in space-time domain, and does not need to consider the frequency of the signal, which makes it have a certain effect on the suppression of seismic low-frequency noise in desert areas. Considering the direct relationship between the weight thresholds and the noise standard deviation, 25
Journal Pre-proof
this paper uses the noise estimation method based on geometric texture to obtain the noise standard deviation which is close to the real one, so as to obtain more accurate weight thresholds. Considering that the different singular values can also represent the effective signals of different energies, the method of truncating singular values is used to reduce the loss of effective signals in the iterative process. The comparison with Wavelet transform, f-x deconvolution and WNNM shows that the
of
proposed method has the best effect on the suppression of random noise in desert areas and the
ro
recovery of effective signals.
-p
Although the proposed method in this paper has a good effect on suppressing low-frequency noise
re
in desert areas, it is necessary to manually set the most suitable size of the search window and seismic texture patch when dealing with different seismic records, which cannot be adaptive. This
ACKNOWLEDGMENTS
na
lP
adds a lot of work, and we want to improve it in the future.
Jo ur
This work was supported by the key project of National Natural Science Foundation of China (Grant No. 41730422), The key Science Foundation of the Department of Science and Technology of Jilin Province (Grant No. 20180201081SF and 20190303082SF).
Reference Abma, R., & Claerbout, J., 1995, Lateral prediction for noise attenuation by t-x and f-x techniques: 26
Journal Pre-proof
Geophysics, 60, no. 6, 1887–1896.
https://doi.org/10.1190/1.1443920
Amani, S., Gholami, A. & Javaheri Niestanak, A., 2017, Seismic random noise attenuation via 3D block matching: Journal of Applied Geophysics, 136, 353-363.
of
https://doi.org/10.1016/j.jappgeo.2016.11.014
ro
Anvari, R. , Kahoo, A. R. , Mohammadi, M. , Khan, N. A. , & Chen, Y., 2019, Seismic random noise
-p
attenuation using sparse low-rank estimation of the signal in the time–frequency domain: IEEE
re
Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 12, no. 5,
lP
1612-1618.
na
https://doi.org/10.1109/JSTARS.2019.2906360
Anvari, R. , Siahsar, M. A. N. , Gholtashi, S. , Kahoo, A. R. , & Mohammadi, M., 2017, Seismic
Jo ur
random noise attenuation using synchrosqueezed wavelet transform and low-rank signal matrix approximation: IEEE Transactions on Geoscience and Remote Sensing, 55, no. 11, 6574-6581.
https://doi.org/10.1109/TGRS.2017.2730228
Anvari, R. , Mohammadi, M. , & Kahoo, A. R., 2018, Enhancing 3-d seismic data using the t-svd and optimal shrinkage of singular value: IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 12, no. 1, 382-388.
https://doi.org/10.1109/JSTARS.2018.2883404
Bai, M., & Wu, J., 2019, Iterative sparse deconvolution using seislet-domain constraint: Journal of 27
Journal Pre-proof
Seismic Exploration, 28, no. 1, 73-88.
Bekara, M., & van der Baan, M., 2009, Random and coherent noise attenuation by empirical mode decomposition: Geophysics, 74, no. 5, V89–V98.
https://doi.org/10.1190/1.3157244
of
Cai, J.F., Candès, E., & Shen, Z.W., 2010, A singular value thresholding algorithm for matrix
ro
completion: SIAM Journal on Optimization, 20, no.4, 1956–1982.
-p
https://doi.org/10.1137/080738970
re
Chen, Y.K., 2016, Dip-separated structural filtering using seislet transform and adaptive empirical
na
https://doi.org/10.1093/gji/ggw165
lP
mode decomposition based dip filter: Geophysical Journal International, 206, no. 1, 457-469.
Jo ur
Chen, K., & Sacchi, M.D., 2015, Robust reduced-rank filtering for erratic seismic noise attenuation: Geophysics, 80, no. 1, V1-V11.
https://doi.org/10.1190/GEO2014-0116.1
Chen, Y.K., Zhang, D., Jin, Z.Y., Chen, X.H., Zu, S.H., Huang, W.L., & Gan, S.W., 2016, Simultaneous denoising and reconstruction of 5-D seismic data via damped rank-reduction method: Geophysical Journal International, 206, no. 3, 1695-1717.
https://doi.org/10.1093/gji/ggw230
Chen, Y.K., Zhou, Y.T., Chen, W., Zu, S.H., Huang, W.L., & Zhang, D., 2017, Empirical low-rank approximation for seismic noise attenuation: IEEE Transactions On Geoscience And Remote 28
Journal Pre-proof
Sensing, 55, no. 8, 4696–4711.
https://doi.org/10.1109/TGRS.2017.2698342
Chen, W., & Song, H., 2018, Automatic noise attenuation based on clustering and empirical wavelet transform: Journal of Applied Geophysics. 159, 649-665.
of
https://doi.org/10.1016/j.jappgeo.2018.09.025
ro
Dalai, B., Kumar, P., Yuan, X.H., 2019, De-noising receiver function data using the Seislet
-p
Transform: Geophysical Journal International, 217, no. 3, 2047-2055.
re
https://doi.org/10.1093/gji/ggz135
lP
Gaci, S., 2014, The use of wavelet-based denoising techniques to enhance the first-arrival picking on
na
seismic traces: IEEE Transactions On Geoscience And Remote Sensing, 52, no. 8, 4558–4563.
Jo ur
https://doi.org/10.1109/TGRS.2013.2282422
Gu, S.H., Zhang, L., Zuo, W.M., & Feng, X.C., 2014, Weighted Nuclear Norm Minimization with Application to Image Denoising: IEEE Conference on Computer Vision and Pattern Recognition(CVPR), 2862-2869.
https://doi.org/10.1109/CVPR.2014.366
Gorszczyk, A., Adamczyk, A.,Michal, M., 2014. Application of curvelet denoising to 2D and 3D seismic data — practical considerations: Journal of Applied Geophysics. 105, 78–94.
https://doi.org/10.1016/j.jappgeo.2014.03.009
29
Journal Pre-proof
Gulunay, N., 1986, FXDECON and complex Wiener prediction filter: 56th Annual International Meeting, SEG, Expanded Abstracts, 279–281.
Herrmann, F.J., & Hennenfent, G., 2008, Non-parametric seismic data recovery with curvelet frames: Geophysical Journal International, 173, no. 1, 233-248.
of
https://doi.org/10.1111/j.1365-246X.2007.03698.x
Kopsinis, Y., & McLaughlin, S., 2009, Development of EMD-based denoising methods inspired by
-p
ro
wavelet thresholding: IEEE Transactions on Signal Processing, 57, no. 4, 1351–1362.
re
https://doi.org/10.1109/TSP.2009.2013885
lP
Li, J., Wang, D.X., Ji, S., Li, Y., & Qian, Z.H., 2017, Seismic noise suppression using weighted
na
nuclear norm minimization method: Journal of Applied Geophysics, 146, 214–220.
Jo ur
https://doi.org/10.1016/j.jappgeo.2017.09.013
Li, G.H., Li, Y., & Yang, B.J., 2017, Seismic Exploration Random Noise on Land: Modeling and Application to Noise Suppression: IEEE Transactions On Geoscience And Remote Sensing, 55, no. 8, 4668-4681.
https://doi.org/10.1109/TGRS.2017.2697444
Liu, W., Cao, S.Y., & Wang, Z.M., 2017, Application of variational mode decomposition to seismic random noise reduction: Journal of Geophysics and Engineering, 14, no. 4, 888–899.
https://doi.org/10.1088/1742-2140/aa6b28
Liu, X.H., Tanaka, M., Okutomi, M., 2012, Noise level estimation using weak textured patches of a 30
Journal Pre-proof
single noisy image: 19th IEEE International Conference on Image Processing(ICIP), 665–668.
https://doi.org/10.1109/ICIP.2012.6466947
Liu, Y.P., Li, Y., Lin, H.B., & Ma, H.T., 2014, An Amplitude-Preserved Time–Frequency Peak Filtering Based on Empirical Mode Decomposition for Seismic Random Noise Reduction: IEEE
of
Geoscience and Remote Sensing Letters, 11, no. 5, 896-900.
ro
https://doi.org/10.1109/LGRS.2013.2281202
re
Journal of Applied Geophysics, 151, 246–262.
-p
Liu, Z., Y. Chen, and J. Ma, 2018, Ground roll attenuation by synchrosqueezed curvelet transform:
lP
https://doi.org/10.1016/j.jappgeo.2018.02.016
na
Lari, H.H., & Gholami, A., 2014. Curvelet-tv regularized bregman iteration for seismic random noise
Jo ur
attenuation: Journal of Applied Geophysics. 109, 233–241.
https://doi.org/10.1016/j.jappgeo.2014.08.005
Sacchi, M.D., Ulrych, T.J., & Walker, C.J., 1998, Interpolation and extrapolation using a high-resolution discrete Fourier transform: IEEE Transactions on Signal Processing, 46, no. 1, 31-38.
https://doi.org/10.1.1.154.4747
Souza, R., Lumley, D., & Shragge, J., 2016, Estimation of reservoir fluid saturation from 4D seismic data: Effects of noise on seismic amplitude and impedance attributes: Journal of Geophysics and Engineering, 14, no. 1, 51–68. 31
Journal Pre-proof
https://doi.org/10.1088/1742-2132/14/1/51
Shao, O.Y., Wang, L.L., Hu, X.Y., & Long, Z.D, 2019, Seismic random noise attenuation using nonlocal means via smooth patch ordering: Journal of Applied Geophysics, 167, 1-10.
https://doi.org/10.1016/j.jappgeo.2019.05.005
ro
denoising. Journal of Applied Geophysics. 69, 103–115.
of
Shan, H., Ma, J. & Yang, H., 2009. Comparisons of wavelets, contourlets and curvelets in seismic
-p
https://doi.org/10.1016/j.jappgeo.2009.08.002
re
Tang, N., Zhao, X., Li, Y, & Zhu, D, 2018, Adaptive threshold shearlet transform for surface
lP
microseismic data denoising: Journal of Applied Geophysics, 153, 64-74.
na
https://doi.org/10.1016/j.jappgeo.2018.03.019
Jo ur
Tang, G., & Ma, J.W., 2010, Application of total-variation-based curvelet shrinkage for three-dimensional seismic data denoising: IEEE Geoscience and Remote Sensing Letters, 8, no. 1, 103–107.
https://doi.org/10.1109/LGRS.2010.2052345
Wu, X.M., 2017, Structure-, stratigraphy- and fault-guided regularization in geophysical inversion: Geophysical Journal International, 210, no. 1, 184-195.
https://doi.org/10.1093/gji/ggx150
Wang, H, Cao, S.Y., Jiang, K.K., Wang, H, & Zhang, Q.C., 2019, Seismic data denoising for complex structure using BM3D and local similarity: Journal of Applied Geophysics, 170. 32
Journal Pre-proof
https://doi.org/10.1016/j.jappgeo.2019.04.018
Xue, Y., Man, M., Zu, S., Chang, F., Chen, Y., 2017. Amplitude-preserving iterative deblending of simultaneous source seismic data using high-order radon transform. Journal of Applied Geophysics. 139, 79–90.
of
https://doi.org/10.1016/j.jappgeo.2017.02.010
Yu, S.W., & Ma, J.W., 2018, Complex variational mode decomposition for sloppreserving denoising:
re
https://doi.org/10.1109/TGRS.2017.2751642
-p
ro
IEEE Transactions On Geoscience And Remote Sensing, 56, no. 1, 586–597.
lP
Yuan, S.Y., Liu, J.W., Wang, S.X., Wang, T.Y., & Shi, P.D., 2018, Seismic waveform classification
Letters, 15, no. 2, 272–276.
na
and first-break picking using convolution neural networks: IEEE Geoscience and Remote Sensing
Jo ur
https://doi.org/10.1109/LGRS.2017.2785834
Zhai, M.Y., 2014. Seismic data denoising based on the fractional Fourier transformation: Journal of Applied Geophysics. 109 , 62–70.
https://doi.org/10.1016/j.jappgeo.2014.07.012
Zhang, Z.J., & Klemperer, S. L., 2005, West-east variation in crustal thickness in northern Lhasa block central Tibet from deep seismic sounding data: Journal of Geophysical Research, 110.
https://doi.org/10.1029/2004JB003139
Zhang, Z.J., & Klemperer, S. L., 2010, Crustal structure of the Tethyan Himalaya, southern Tibet: 33
Journal Pre-proof
newconstraints from old wide-angle seismic data: Geophysical Journal International, 181, no. 3,
1247–1260.
https://doi.org/10.1111/j.1365-246x.2010.04578.x
Zhou, Y.T., & Zhu, Z.L., 2019, A hybrid method for noise suppression using variational mode
of
decomposition and singular spectrum analysis: Journal of Applied Geophysics, 161, 105-115.
Jo ur
na
lP
re
-p
ro
https://doi.org/10.1016/j.jappgeo.2018.10.025
34
Journal Pre-proof CRediT author statement
Juan Li: Validation, Conceptualization, Supervision. Wei Fan: Methodology, Software, Writing-Original Draft. Yue Li: Data Curation, Supervision, Project administration. Baojun Yang: Data Curation, Resources.
Jo ur
na
lP
re
-p
ro
of
Changgang Lu: Conceptualization, Writing-Review & Editing.
35
Juan Li, Wei Fan, Yue Li, Baojun Yang, Changgang Lu PII:
S0926-9851(19)30480-X
DOI:
https://doi.org/10.1016/j.jappgeo.2019.103926
Reference:
APPGEO 103926
To appear in:
Journal of Applied Geophysics
Received date:
25 May 2019
Revised date:
24 December 2019
Accepted date:
25 December 2019
Please cite this article as: J. Li, W. Fan, Y. Li, et al., Desert seismic noise suppression based on an improved low-rank matrix approximation method, Journal of Applied Geophysics(2019), https://doi.org/10.1016/j.jappgeo.2019.103926
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2019 Published by Elsevier.
Journal Pre-proof
Desert seismic noise suppression based on an improved low-rank matrix approximation method Juan Li,1 Wei Fan,1 Yue Li,1* Baojun Yang2 and Changgang Lu3 1
Department of Information, College of Communication Engineering, Jilin University, Changchun 130012, China. E-mail: [email protected] Department of Geophysics, Jilin University, Changchun 130026, China
3
College of Automotive Engineering, Jilin University, Changchun, China
of
2
ro
Abstract
-p
The suppression of random noise is a crucial step before seismic data analysis. Random noise in
re
desert areas has the characteristics of low frequency and non-stationary, and there is serious spectrum
lP
aliasing between random noise and effective signals, which makes it difficult to suppress such noise. In recent years, some methods based on signal rank minimization have achieved remarkable results
na
in seismic random noise suppression. Since the implementation of low rank matrix approximation is
Jo ur
an iterative process, noise estimation is an indispensable step before each iteration, but also an important step. The noise estimation method previously used is to calculate the residuals of the original noisy patch data and the corresponding iterative denoising version, which is intuitively considered as the filtered noise. This method may be very inaccurate in the case of high noise levels or complex seismic records. In this paper, a noise estimation method based on geometric texture is introduced to estimate the noise level by selecting weak textured patches in all seismic texture patches. At the same time, we reduce the loss of effective signals by truncating the singular values in each iteration. Experiments on both synthetic and field seismic data show that this method has better effect on suppressing random noise in desert areas.
1
Journal Pre-proof Key Words: Desert seismic noise; low rank matrix approximation; noise estimation; noise suppression; truncated singular value;
1 INTRODUCTION In the process of collecting seismic signals, it is inevitable that it will be affected by various types
of
of noise, which will have strong interference on the subsequent analysis of seismic data. If the
ro
random noise is suppressed, the subsequent analysis of seismic data will be more accurate (Zhang and Klemperer, 2005, 2010; Souza et al., 2016; Yuan et al., 2018; Li et al., 2017; Chen et al., 2016;
-p
Wu, 2017). Therefore, noise removal plays an important role in increasing the quality of seismic
re
exploration data. The frequency of random noise in desert areas is relatively low, and the frequency
lP
bands of effective signals and random noise are seriously overlapped, which makes some methods
na
suitable for suppressing white noise and separating signals and noise in frequency domain ineffective
Jo ur
on suppressing the random noise in desert areas.
Previously, abundant methods have been proposed and applied to suppress random noise in real seismic data. Sparse transform-based methods, such as Fourier transform (Sacchi et al., 1998; Zhai, 2014), wavelet transform (Gaci, 2014; Chen and Song, 2018; Shan et al., 2009), curvelet transform (Tang and Ma, 2010; Herrmann and Hennenfent, 2008; Gorszczyk et al., 2014; Liu et al., 2018; Lari and Gholami, 2014), shearlet transform ( Tang et al., 2018), seislet transform (Bai and Wu, 2019; Dalai et al., 2019), Radon transform (Xue et at., 2017) transform the seismic data into the sparse domain to obtain the corresponding sparse coefficients, then, the thresholding process is performed on the sparse coefficients, and the processed sparse coefficients are transformed back into the space-time domain to achieve the purpose of suppressing noise. Mode decomposition-based methods, 2
Journal Pre-proof
such as empirical mode decomposition (EMD), variational mode decomposition (VMD) and their improvements (Bekara and van der Baan, 2009; Kopsinis and McLaughlin, 2009; Liu et al., 2014; Chen, 2016; Yu and Ma, 2018; Liu et al., 2017; Zhou and Zhu, 2019), mainly decompose noisy seismic data into many different components, then, the components containing effective signals are selected and the components of random noise are discarded. Prediction-based methods , such as f-x
of
deconvolution and t-x predictive filtering (Gulunay, 1986; Abma and Claerbout, 1995), mainly utilize the predictability of effective signals to construct predictive filters to suppress noise and
ro
enhance signals. Some methods that use non-local self-similarity of seismic records are also used to
-p
process seismic records(Shao et al., 2019; Wang et al., 2019; Amani et al. 2017) . Some methods can
re
recover the signal effectively by use the low rank feature of seismic signal in time-frequency domain,
lP
for example, Rasoul Anvari et al. (2019) consider the ideal seismic signal in the transform domain as
na
both sparse and low rank. They first compute the short-time Fourier transform of the noisy seismic signal, then use a nonconvex penalty function to estimate the sparse low-rank matrix, and finally use
Jo ur
the estimated sparse low-rank matrix to synthesize the seismic signal. Rasoul Anvari et al. (2017) first use the synchrosqueezed wavelet transform to transform the seismic data trace by trace into the sparse domain to obtain the sparse time-frequency representation, and then use the Optshrink algorithm to decompose the sparse time-frequency representation into semilow-rank and sparse components, finally back-transform the semilow-rank component to the time domain to recover the denoised seismic trace using inverse synchrosqueezed wavelet transform. The above existing methods have played a certain role in dealing with large sets of oil and gas layers, random noise reduction in medium and shallow exploration, and improving the signal-to-noise ratio of seismic records, and obtained good results in seismic exploration engineering. 3
Journal Pre-proof
In recent years, a large number of scholars have conducted in-depth research on the low-rank method of signals. It can be considered that seismic signals are of low rank after some methods of signal rearrangement steps. For example, Rasoul Anvari et al. (2019) model the 3-D seismic data as a tensor data of size n1 n2 n3 , which can be expressed as a combination of a low rank component and a sparse component, then compute the Fourier transform of the tensor along the third dimension
of
of the tensor, and shrink the singular values in the Fourier domain to extract the low rank component. They repeat the above steps until the Frobenius norm of the error matrix reaches the desired value,
ro
the matrix constructed by stacking non-local similar patches from the noisy seismic signals is of low
-p
rank in the low-dimensional subspace of the given high-dimensional space (Li et al., 2017). EMD
re
can empirically decompose a multi-dip seismic image that is not of low rank into multiple single-dip
lP
seismic images that are low-rank individually (Chen et al., 2017). The Hankel matrix constructed
na
from seismic signals will have a lower rank (Chen and Sacchi, 2015). The denoising mechanism of these methods is also different. Li et al. (2017) takes advantage of the prior knowledge that the
Jo ur
singular values have clear physical meanings and should be treated differently, and assigns different weights to singular values of different amplitudes. Chen et al. (2017) constructs a hankel matrix for each empirical component after EMD decomposition, and then apply SVD to each hankel matrix, then retain the largest singular value of the same number as the rank of modal component and set others to zero. Chen and Sacchi (2015) decomposes the Hankel matrix constructed by seismic signals into the product of two lower dimensional factor matrices, then the two matrices are solved by weighted cost function.
Noise estimation is a very important step of the low rank matrix approximation denoising framework, the estimated noise standard deviation will be used to calculate the singular value of the 4
Journal Pre-proof
pure signal and the corresponding weight thresholds, which are decisive parameters for denoising, thus the accuracy of the noise estimation is directly related to the quality of the denoising results. The noise estimation method in (Gu et al., 2014) is to calculate the difference between the stacked similar patch matrices of noisy seismic records and the corresponding matrices of current iterative denoising results, which is intuitively assumed to be the noise. The denoising results of this method are greatly
of
influenced by the initial value of noise level, and the results of current iteration will affect the results of next iteration. If the initial value of noise level is not set accurately, or it is not only noise that is
ro
filtered out in a certain iteration, then the denoising results will not be very satisfactory. Especially in
-p
the field seismic data processing, because of the complexity of the field seismic data, it would be
re
inaccurate to use the original noise estimation method.
lP
To make up for this deficiency, we introduce the noise estimation method based on geometric
na
texture (Liu et al., 2012) into the denoising framework of low rank matrix approximation in this
Jo ur
paper. The null hypothesis test is used to select weak textured patches from all seismic texture patches and then principal component analysis (PCA) is performed on these weak textured patches to estimate the noise level. With this method, the denoising results no longer depend on the initial value of the noise level, and the noise estimations before each iteration are no longer dependent on the results of the previous iteration, which are independent of each other, so as to get the estimates closer to the real value. Meanwhile, taking into account the physical meaning of the singular value, the singular value of the effective signal is related to its energy. Generally, the singular value of the signal with higher energy is larger than that of the signal with lower energy. In the original framework, the smaller singular value will become zero after the weight threshold processing, which will result in the loss of low energy signal. So before each iteration, we divide the singular value 5
Journal Pre-proof
matrix into two matrices with the same dimension as the original one, one retaining larger singular value and the other retaining smaller singular value. Two seismic signals can be obtained from these two matrices, one representing high energy signals and the other representing low energy signals. We add this difference between the two seismic signals to the iterative regularization term to reduce the loss of effective signals. Tests with synthetic and field seismic data prove that the method proposed
Jo ur
na
lP
re
-p
ro
of
in this paper has better effect on suppressing seismic random noise in desert areas.
Fig. 1: Denoising records and corresponding filtered records after the 2nd, 4th and 6th iteration of WNNM. (a) Noisy record. (b),(c) the denoising records and corresponding filtered records after the 2nd iteration. (d),(e) the denoising records and corresponding filtered records after the 4th iteration. (f),(g) the denoising records and corresponding filtered records after the 6th iteration. 2 Low Rank Matrix Approximation for Seismic Noise Suppression in Desert Areas
6
Journal Pre-proof
2.1 Theory Noisy seismic data in desert areas can be expressed as follows:
Y X E,
(1)
where Y is the noisy seismic datas, X is the potential clean signals, and E is the random noise with a standard deviation of . Our purpose is to recover X from Y , In the case of complex
of
seismic signals, X itself is not low rank. Using the non-local self-similarity of seismic signals, we
ro
select a main patch PYj in a search window of a certain size, and then select patches similar to the
-p
main patch to stack them into a matrix M Yj . At the same time, the corresponding pure signal matrix
re
M X j and random noise matrix M E j can be found to represent (1) as follows：
lP
M Yj M X j M E j .
(2)
na
After this data rearrangement, it can be considered that the matrix M X j is of low rank, so the
Jo ur
method of low rank matrix approximation can be used to recover M X j from M Yj to achieve the purpose of denoising.
The most primitive method for solving (2) is defined as
ˆ X arg min rank (M X ) ,s.t. M Y M X M E . M j j j j j MX j
(3)
But the objective function in (3) is non-convex, and solving it is an NP-hard problem. However, it can be converted to nuclear norm minimization (NNM) by convex relaxation as follows (Gu et al., 2014; Cai et al., 2010):
7
Journal Pre-proof ˆ X arg min 1 || M Y M X || 2F || M X ||* , M j j j j 2 MXj
(4)
where is a fixed positive constant, ||||F is Frobenius norm and ||M X j||* is the nuclear norm of matrix M X j , || ||* i |λi ()| , i () represents the ith singular value of the matrix. The optimal solution of (4) can be obtained by singular value decomposition (Cai et al., 2010), i.e.,
M Yj UVT，ˆ ii S (ii ) max( ii ,0) ,
of
(5)
(6)
ro
ˆ X Uˆ VT , M j
-p
where ii represents the diagonal element of the matrix , ii i (M Yj ) , It can be observed that
re
NNM performs the same threshold for each singular value, which does not take into account the
lP
practical physical meaning of the singular value. The effective information of the seismic signal is mainly represented by a large singular value, so the larger singular value should be assigned a
na
smaller threshold than the smaller singular value. Gu et al. (2017) proposed weighted nuclear norm
which is defined as
Jo ur
minimization (WNNM), assigning different weight thresholds according to different singular values,
ˆ X arg min 1 || M Y M X || 2F || M X || w,* , M j j j j 2 MXj
(7)
where || ||w,* i |wi λi ()| , w w1, w2 , , wn , wi is the non-negative weight value corresponding to the singular value i () , The relationship between wi and i (M X j ) is as follows:
wi
c m , i (M X j )
(8)
8
Journal Pre-proof where c 0 is a constant and m is the number of patches constituting the matrix M Yj . To avoid the denominator being zero, we set 1016 . It can be seen from (8) that the larger the singular value is, the smaller the corresponding wi is. Although (7) is non-convex when the weight vector
w is in a non-descending order, Gu et al. (2017) prove that there exists a local minimum point which is its optimal solution. The solution of (7) is similar to that of (4). By replacing in (5) with
of
wi , we can get the following expression (9)
ro
ˆ ii S (ii ) max( ii wi ,0) .
-p
The other steps are the same as those of (4).
re
There is a problem in obtaining wi , that is, the singular value of the pure seismic signal cannot be
lP
known in advance. But it can be estimated by the singular value of noisy seismic signal and the
na
variance of noise in noisy seismic signal, which is defined as
i (M X j ) max( i2 (M Yj ) m 2 ,0) ,
Jo ur
(10)
where i (M Y j ) can be obtained by calculation, so the question now is how to get . Our work is also carried out around . From (8), (9), (10), we can see the accuracy of will directly affect the accuracy of i (M X j ) , and then it will affect the value of wi , and ultimately affect the effect of noise suppression.
Algorithm 1 Seismic Denoising Using the proposed method Input: Noisy seismic signal Y , parameter , , L , 0 ˆ 0 Y, X ˆ 0high Y, X ˆ low Y 1: Initialize Y0 Y, X
2: for k 1 : L do 3:
9
k ˆ k (Y X ˆ k ) (X ˆ khigh X ˆ low ) iterative regularization Y k X
k
Jo ur
na
lP
re
-p
ro
of
Journal Pre-proof
2.2 Existing noise estimation method in WNNM
In the WNNM framework, the noise standard deviation is an important parameter. Since WNNM is an iterative process, two noise values can be obtained after each iteration: filtered noise and residual noise. The residual noise standard deviation can be estimated from the original noise standard deviation and the standard deviation of the filtered noise. The estimation process is as 10
Journal Pre-proof
follows. First of all, the method cannot obtain the original noise standard deviation , so must be set artificially. It is assumed that the noise standard deviation filtered by WNNM at the k-th iteration can be expressed as kflt , the standard deviation of residual noise in the record to be processed in the k 1 (k+1)-th iteration can be represented by res . Suppose that the result of PYj after the k-th iterative
of
denoising is Pˆ Yk j 1 , because there is not only noise but also signal in Pˆ Yk j 1 , the standard deviation
ro
k 1 cannot be calculated directly in Pˆ Yk j 1 . However, when there is only noise in the record (i.e. res
-p
PYj Pˆ Ykj1 ) filtered by the k-th iteration, the standard deviation kflt can be directly calculated as
re
follows
(11)
lP
kflt || PYj Pˆ Ykj1 || 2 ,
formula
Jo ur
k 1 res 2 ( kflt ) 2 ,
na
then the standard deviation of the residual noise can be estimated by and kflt by the following
(12)
k 1 where res denotes the estimated standard deviation of the residual noise present at (k+1)-th
iteration. The constant 0 is a scaling factor, heuristically introduced to control the re-estimation
of the standard deviation.
There are two main problems in this noise estimation method: First, the original noise standard deviation needs to be set artificially, and its authenticity cannot be guaranteed. Second, if there is signal in the records filtered by k-th iteration, the calculation of the standard deviation of the filtered noise will be inaccurate, which will affect the estimation of the standard deviation of the residual 11
Journal Pre-proof
noise. In order to observe the filtered noise and residual noise obtained by WNNM in each iteration more intuitively, we use it to process the noisy records in Fig. 1(a), and obtain the denoising records and corresponding filtered records after the 2nd, 4th and 6th iteration respectively, as shown in Fig. 1(b)-(g). It can be seen from the experimental results that there are not only filtered noise but also obvious signals in the filtered records of these iterations.
of
The original noise standard deviation and the standard deviation of the filtered noise obtained
ro
at the kth iteration determine the standard deviation of the residual noise in the record to be
-p
processed in the (k+1)-th iteration. That is to say, if the is set inaccurately or there is not only
re
noise in the filtered records as shown in Fig.1, it will affect the denoising results of the (k+1)-th iteration.Therefore, we introduce a noise estimation method based on geometric texture. The noise
lP
standard deviation is directly estimated according to the geometric structure of different texture
na
patches of the noisy seismic signal. It is not necessary to set the initial value of the noise level
Jo ur
manually, and the effect of denoising in the previous iteration does not need to be considered.
3 Noise estimation based on geometric texture Geometric texture-based noise estimation was first proposed and applied to images by Liu et al. (2012). This estimated method has a good practical application effect in images, and the estimated noise standard deviation is very close to the actual value. So we apply this estimation method to noisy seismic signals. In view of the idea proposed by Liu et al. (2012), the noise level can be estimated from weak textured patches in seismic textured patches, which is defined as
ˆ n2 min (cov(W)) ,
(13)
12
Journal Pre-proof where min () is the minimum eigenvalue of the matrix, W represents the selected weak textured patches, and cov() denotes the covariance matrix. Therefore, the key issue in applying geometric texture-based noise estimation method to seismic signals is how to select weak textured patches from seismic texture patches. First, the gradient covariance matrix Cp j of the seismic texture patch can be written as
Dvp j ,
(14)
of
Cp j GpTj Gp j，Gp j Dhp j
ro
where p j denotes the texture patch matrix in matrix M Yj , D h and Dv represent the horizontal
-p
and vertical derivative operators, respectively. The larger eigenvalues of Cp j often reflect more
re
texture information of p j . Theoretically, the noisy flat seismic patch PM j in the reconstructed
lP
seismic signal M Yj can be decomposed into a perfectly flat seismic patch Pf j and a noise patch
PE j , as shown in follow,
na
PM j Pf j PE j .
(15)
Jo ur
It is easy to know that the gradients of the perfectly flat seismic patch are zero, so the expected value of the gradient covariance matrix of noisy flat seismic patches can be calculated as
E (CPM ) E (CPf ) E (CP ) E (CPE ) j
j
Ej
PET DTh Dh PE j E ( Tj T P D D P E j v h E j
j
PET j DTh Dv PE j ） PETj DTv Dv PE j
.
(16)
E (PET DTh Dh PE ) 0 j j 0 E (PET j DTv Dv PE j ) The statistical properties of two diagonal elements are the same. Liu et al. (2012) proved that the distribution of two diagonal elements PETj DT D PE j (where h or v ) can be approximated by gamma distribution. The shape parameter and the scale parameter of the gamma distribution 13
Journal Pre-proof
are as follows:
N 2 2 ， n tr (DT D ) , 2 N
(17)
where N denotes the number of elements in matrix PE j , n is the noise standard deviation to be estimated, and tr (DT D ) denotes the trace of the matrix DT D . In order to select the weak textured patches in the seismic texture patches, we refer to (Liu et al.
of
2012), also using the null hypothesis test method. The null hypothesis is “ the given patch is a flat
ro
patch with random noise”. When the maximum eigenvalue of the gradient covariance matrix of the
re
F 1 (v, , ) ,
lP
valid..The threshold is calculated by
-p
given seismic texture patch is less than a certain threshold, the hypothesis is valid, otherwise it is not
(18)
na
where F 1 denotes the inverse gamma cumulative distribution function and v denotes the
Jo ur
significance level, which must be set manually, generally set to 0.95~0.99. It can be seen from (17) and (18) that the noise standard deviation is also an indispensable parameter for calculating the threshold .
Liu et al (2012) proposes an iterative method to estimate the final noise level. We combine it with seismic signals and summarize it as follows: Firstly, the initial noise standard deviation n( 0) is estimated from the covariance matrix generated by all texture patches in M Yj by (13), and the initial threshold （0） is calculated by (18) according to n( 0) .Weak textured patches W (1) can be selected by null hypothesis test using the threshold （0）, and then n(1) can be estimated by (13). Repeat the above process until the estimated noise level n tends to be stable. We believe that the 14
Journal Pre-proof final n is the closest to the true noise level. To verify the effectiveness of this noise estimation, we add the actual desert noise with different noise levels to the synthetic seismic record in Fig. 2(a), and then use this estimation method to estimate the standard deviation of the noise we added. The true values and estimates of noise standard deviation are listed in Table 1. It can be seen that the noise standard deviation estimated by
of
this noise estimation method is very close to the true values.
1
2
3
4
Estimates
1.0005
1.9975
2.9920
3.9873
5
-p
Truth
ro
Table 1: Comparison of the true values and estimates of noise standard deviation. 15
20
9.9628
14.9432
19.9239
re
4.9830
10
lP
4 Reducing signal loss by truncating singular values
na
Low rank matrix approximation is an iterative denoising process, so it will inevitably lead to the loss of some seismic signals during the iterative process. We consider the physical meaning of the
Jo ur
singular value of the seismic signal. Generally, the singular values of the high-energy effective signals are relatively large, but there are some low-energy signals with small singular values. The singular values of this part may be changed to zero or very close to zero after the threshold processing by (9) in the iterative process, which will result in the loss of this part of the signal. In order to reduce this loss, in each iteration, we use threshold to separate the larger and smaller singular values of the signal as follows: ˆ ˆ high ˆ low ,
(19)
where ˆ high and ˆ low are matrices that retain only larger and smaller singular values in ˆ , respectively. The optimal solution of (7) can be written as 15
Journal Pre-proof ˆ X Uˆ VT U(ˆ high ˆ low )VT M ˆ high ˆ low M Xj M Xj . j
(20)
ˆ high and X ˆ low by aggregating all patch Then, we can obtain two seismic signal records X ˆ high ˆ low ˆ high and X ˆ low can reflect the seismic signal matrices M and M X j . The difference between X Xj
ˆ high and X ˆ low but with different energy levels. We add the difference between contained in both X high and low energy seismic signals to the iteration regularization term to reduce some unnecessary signal loss in the iteration denoising process.
Jo ur
na
lP
re
-p
ro
of
For a clearer description, the improved low rank matrix approximation algorithm is summarized in algorithm 1.
Fig. 2: Denoising results of synthetic seismic record. (a) Ideal signal. (b)Noisy signal. (c) Denoising result of Wavelet Transform. (d) Denoising result of f-x deconvolution. (e) Denoising result of WNNM. (f) Denoising result of the proposed method.
16
ro
of
Journal Pre-proof
re
-p
Fig. 3: Difference between before and after denoising. (a) Added noise. (b) Difference between before and after denoising of Wavelet Transform. (c) Difference between before and after denoising of f-x deconvolution. (d) Difference between before and after denoising of WNNM. (e) Difference between before and after denoising of the proposed method.
lP
5 Application to noise suppression for seismic records in desert areas
na
5.1 Experiments on synthetic seismic data
To verify the effectiveness of the proposed method, we first test it on the synthetic seismic data. The
Jo ur
synthetic seismic model has a total of seven effective events, which are generated by Ricker wavelets with dominant frequencies 30 Hz and 25 Hz, as shown in Fig. 2(a). Then we add real random noise in desert areas to the synthetic seismic data, so that the signal-to-noise ratio of the synthetic record is -5.1402 dB, as shown in Fig. 2(b). When processing this record, we set the size of the search window to 40 40 , and the size of the seismic texture patch to 1313 , this setting is only suitable for processing records at this noise level. When processing synthetic records with different noise levels, the sizes of search window and seismic patch are set differently. When the noise level is low, they are set smaller, and when the noise level is high, they are set larger. At the same time, we set parameters
0.1 , 0.4 , and 0.5 to process synthetic records at all signal-to-noise ratios. 17
Journal Pre-proof
In contrast experiments, Wavelet Transform, f-x deconvolution and WNNM are used to process the same synthetic records to compare with the method proposed in this paper. We adjust the parameters of each method to make the denoising result the best. The denoising results of these three methods are shown in Fig. 2(c), (d), (e). The denoising results of the proposed method are shown in Fig. 2(f). At the same time, to observe their retention of the signal, we obtained the difference
of
between before and after denoising, as shown in Fig. 3. To observe and compare the denoising results of these four methods more comprehensively, we achieve a more detailed knowledge on
ro
signal-preserving performance by analyzing a single trace. For this purpose, we extract the third trace
-p
of synthetic records for comparison, as shown in Fig. 4. Fig. 4(a)-(d) are a single trace denoising
re
result of wavelet transform, f-x deconvolution, WNNM and the proposed method in this paper,
lP
where the black line represents the ideal signal, the blue line represents the noisy signal, and the red
na
line represents the denoised signal. Then we use FK spectrum to analyze the denoising results of these four methods in the frequency domain. Fig. 5(a)-(f) represents the FK spectrum corresponding
Jo ur
to the denoising result in Fig 2(a)-(f), and Fig. 6(a)-(e) represents the FK spectrum corresponding to the difference between before and after denoising in Fig. 3(a)-(e).
Fig. 4: Comparison on the 3th trace of denoised records. (a) Wavelet transform. (b) f-x deconvolution. (c) WNNM. (d) The proposed method.
18
ro
of
Journal Pre-proof
Jo ur
na
lP
re
-p
Fig. 5: The FK spectrum of denoising results. (a) Ideal signal. (b) noisy signal. (c) Denoising result of Wavelet Transform. (d) Denoising result of f-x deconvolution. (e) Denoising result of WNNM. (f) Denoising result of the proposed method.
Fig. 6: The FK spectrum of the differences between before and after denoising. (a) Added noise. (b) Wavelet Transform. (c) f-x deconvolution. (d) WNNM. (e) The proposed method.
Table 2: SNRs of different denoising results. SNR(dB) the original data
Wavelet transform
f-x deconvolution
19
WNNM
The proposed method
Journal Pre-proof
11.1813
12.4269
17.0585
17.3263
7.0604
10.5032
12.1593
16.9041
17.2641
5.0581
9.5991
11.8136
16.4879
17.1718
3.0859
8.4850
11.3862
15.4712
17.0639
1.0935
7.1954
9.9314
11.2022
15.6333
-1.0050
5.5561
8.9776
-3.0694
3.8081
7.8492
-5.1402
1.9518
-7.0256
0.1969
-9.0900
-1.7741
14.0166
8.7879
12.8456
6.5472
10.7549
5.2267
3.3723
7.0238
3.6585
1.2603
4.5704
-p
ro
10.5486
re
of
9.0389
Jo ur
na
lP
7.0543
After the above all-round and multi-angle analysis, it can be seen that the noise suppression of wavelet transform is not complete and the effective signals are not recovered well; The noise suppression of f-x deconvolution is not very thorough, the noise in the same frequency band as the 20
Journal Pre-proof
effective signals is preserved, and some signals are lost; Although the advantage in noise suppression of it is very obvious, WNNM also has a large amount of signal loss; The proposed method suppresses the noise thoroughly, the recovery of the effective signal is the best and the processing result is the closest to the ideal signal.
Next, we use these four methods to process synthetic records of random noise with different levels
of
in desert areas. The signal-to-noise ratio (SNR) and mean square error (MSE) of denoising results are
ro
calculated to quantitatively analyze the effectiveness of the proposed method. SNR and MSE are
-p
calculated as follows:
i 1
N
| Xˆ (i) X(i) | i 1
1 N
N
(Xˆ (i) X(i))2 , i 1
(21)
(22)
na
MSE
), 2
lP
SNR(dB) 10 log10（
| X(i) |2
re
N
MSE
Jo ur
ˆ is the denoised signal. where X is the ideal signal and X
the original SNR
Wavelet transform
Table 3: MSEs of different denoising results
f-x deconvolution
WNNM
The proposed method
9.0389
0.0191
0.0144
0.0049
0.0047
7.0604
0.0224
0.0153
0.0051
0.0047
5.0581
0.0276
0.0166
0.0056
0.0048
3.0859
0.0356
0.0183
0.0071
0.0049
21
Journal Pre-proof
0.0479
0.0255
0.0191
0.0069
-1.0050
0.0699
0.0318
0.0222
0.0100
-3.0694
0.1046
0.0412
0.0332
0.0131
-5.1402
0.1604
0.0495
0.0557
0.0211
-7.0256
0.2402
0.0754
0.1156
0.0499
-9.0900
0.3781
0.1082
of
1.0935
0.0877
re
-p
ro
0.1880
The resulting SNR and MSE for synthetic seismic data are given in Table 2 and 3. The higher the
lP
SNR and the lower the MSE indicates that the denoising result is better. It can be seen from the two
na
tables that the SNR of the proposed method is always the highest and MSE is always the smallest,
Jo ur
which shows that the results of the proposed method are the best.
5.2 Experiments on field seismic data In order to verify the practicability of the proposed method, we use it to process field seismic data in desert areas. We have considered a common shot point seismic gather consisting of 127 traces and 2002 samples per trace for testing, as shown in Fig. 7(a). It can be seen that the reflection events on both sides of the record are almost obscured by random noise and are difficult to distinguish. Similarly, Wavelet transform, f-x deconvolution, WNNM and the proposed method are used to process the record. The results of these four methods are shown in fig. 7(b)-(e). When processing this record, we set the size of search window and seismic texture patch in the proposed method to be 22
Journal Pre-proof
30 30 and 7 7 , respectively, and adjust the results of other methods to the best. Similarly, we also draw the difference between before and after denoising of these methods to analyze their signal retention, as shown in fig. 8. From fig. 7(b)-(e) and Fig. 8, it can be seen that the effect on low-frequency noise suppression in Wavelet transform is not satisfactory, and the effective signals have been distorted. Although f-x deconvolution has a certain effect on the low-frequency
of
noise suppression as a whole, it does not restore the effective signals on both sides of the record.
ro
WNNM is not very good at suppressing low-frequency noise in the deep part of the record. The
-p
proposed method in this paper can suppress the low-frequency noise perfectly, and the effective
Jo ur
na
lP
re
signal on both sides of the record can still be recovered well.
Fig. 7: Results obtained with field seismic data. (a) Field seismic data (common shot point seismic data). (b) Denoising result of Wavelet Transform. (c) Denoising result of f-x deconvolution. (d) Denoising result of WNNM. (e) Denoising result of the proposed method.
23
Journal Pre-proof Fig. 8: Difference between before and after denoising. (a) Difference between before and after denoising of Wavelet Transform. (b) Difference between before and after denoising of f-x deconvolution. (c) Difference between before and after denoising of WNNM. (d) Difference between
ro
of
before and after denoising of the proposed method.
-p
Fig. 9: Partial enlarged views. (a) Original enlarged record. (b) Filtered enlarged record after using Wavelet transform. (c) Filtered enlarged record after using f-x deconvolution. (d) Filtered enlarged
re
record after using WNNM. (e) Filtered enlarged record after using the proposed method.
lP
In order to evaluate the quality of the denoising results of these four methods more clearly, we
na
enlarged the record between trace 77-110 and samples 1000-1900 as shown in Fig. 9. Fig. 9(a)-(e) are partial enlarged views corresponding to Fig. 7(a)-(e) respectively. By comparing the partial
Jo ur
enlargement views, it is more obvious that the proposed method in this paper can suppress the low-frequency noise in desert areas better than the other three methods, and can recover the effective signals well.
Table 3: The computation time of different methods in Fig.2(b) Wavelet transform Time(s)
0.1415
f-x deconvolution 0.0195
WNNM 163.5011
The proposed method 419.2504
6 Dicussion
Our algorithm mainly makes use of the non-local self-similarity of the seismic signal in the 24
Journal Pre-proof
space-time domain. The frequency of the signal does not affect this characteristic. In other words, as long as the stacked similar block matrices satisfy the low-rank property, the value of the dominant frequency of the signal will not have a great impact on the performance of our algorithm. Although our algorithm has a great advantage in suppressing random noise compared with several other methods, but for other types of noise, such as coherent noise, because the noise has a certain
of
correlation, it will cause not only signals are similar, but also noise is similar in the stacked similar block matrix, which will reduce the effectiveness of our method. And our algorithm has a cost in
ro
computation time because it needs to estimate the noise level of each stacked similar block matrix.
-p
We calculated the time spent in processing the record in Fig. 2(b) by the methods mentioned in this
re
paper, which are listed in Table 3. It can be seen that the computation time of our algorithm is much
lP
higher than other methods. In addition, under some geological conditions, the amplitude of the
na
seismic reflection signal will change with the distance between the shot point and the detector, but the presence of noise will make this change blurred and difficult to observe. In this case, our
Jo ur
algorithm can suppress the noise while retaining this change in amplitude, which is more conducive to AVO analysis. So our algorithm can be used as a pre-processing step before AVO analysis, providing a clearer and more effective seismic record for AVO analysis.
7 Conclusion
The method of low rank matrix approximation mainly considers the non-local similarity of seismic signals in space-time domain, and does not need to consider the frequency of the signal, which makes it have a certain effect on the suppression of seismic low-frequency noise in desert areas. Considering the direct relationship between the weight thresholds and the noise standard deviation, 25
Journal Pre-proof
this paper uses the noise estimation method based on geometric texture to obtain the noise standard deviation which is close to the real one, so as to obtain more accurate weight thresholds. Considering that the different singular values can also represent the effective signals of different energies, the method of truncating singular values is used to reduce the loss of effective signals in the iterative process. The comparison with Wavelet transform, f-x deconvolution and WNNM shows that the
of
proposed method has the best effect on the suppression of random noise in desert areas and the
ro
recovery of effective signals.
-p
Although the proposed method in this paper has a good effect on suppressing low-frequency noise
re
in desert areas, it is necessary to manually set the most suitable size of the search window and seismic texture patch when dealing with different seismic records, which cannot be adaptive. This
ACKNOWLEDGMENTS
na
lP
adds a lot of work, and we want to improve it in the future.
Jo ur
This work was supported by the key project of National Natural Science Foundation of China (Grant No. 41730422), The key Science Foundation of the Department of Science and Technology of Jilin Province (Grant No. 20180201081SF and 20190303082SF).
Reference Abma, R., & Claerbout, J., 1995, Lateral prediction for noise attenuation by t-x and f-x techniques: 26
Journal Pre-proof
Geophysics, 60, no. 6, 1887–1896.
https://doi.org/10.1190/1.1443920
Amani, S., Gholami, A. & Javaheri Niestanak, A., 2017, Seismic random noise attenuation via 3D block matching: Journal of Applied Geophysics, 136, 353-363.
of
https://doi.org/10.1016/j.jappgeo.2016.11.014
ro
Anvari, R. , Kahoo, A. R. , Mohammadi, M. , Khan, N. A. , & Chen, Y., 2019, Seismic random noise
-p
attenuation using sparse low-rank estimation of the signal in the time–frequency domain: IEEE
re
Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 12, no. 5,
lP
1612-1618.
na
https://doi.org/10.1109/JSTARS.2019.2906360
Anvari, R. , Siahsar, M. A. N. , Gholtashi, S. , Kahoo, A. R. , & Mohammadi, M., 2017, Seismic
Jo ur
random noise attenuation using synchrosqueezed wavelet transform and low-rank signal matrix approximation: IEEE Transactions on Geoscience and Remote Sensing, 55, no. 11, 6574-6581.
https://doi.org/10.1109/TGRS.2017.2730228
Anvari, R. , Mohammadi, M. , & Kahoo, A. R., 2018, Enhancing 3-d seismic data using the t-svd and optimal shrinkage of singular value: IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 12, no. 1, 382-388.
https://doi.org/10.1109/JSTARS.2018.2883404
Bai, M., & Wu, J., 2019, Iterative sparse deconvolution using seislet-domain constraint: Journal of 27
Journal Pre-proof
Seismic Exploration, 28, no. 1, 73-88.
Bekara, M., & van der Baan, M., 2009, Random and coherent noise attenuation by empirical mode decomposition: Geophysics, 74, no. 5, V89–V98.
https://doi.org/10.1190/1.3157244
of
Cai, J.F., Candès, E., & Shen, Z.W., 2010, A singular value thresholding algorithm for matrix
ro
completion: SIAM Journal on Optimization, 20, no.4, 1956–1982.
-p
https://doi.org/10.1137/080738970
re
Chen, Y.K., 2016, Dip-separated structural filtering using seislet transform and adaptive empirical
na
https://doi.org/10.1093/gji/ggw165
lP
mode decomposition based dip filter: Geophysical Journal International, 206, no. 1, 457-469.
Jo ur
Chen, K., & Sacchi, M.D., 2015, Robust reduced-rank filtering for erratic seismic noise attenuation: Geophysics, 80, no. 1, V1-V11.
https://doi.org/10.1190/GEO2014-0116.1
Chen, Y.K., Zhang, D., Jin, Z.Y., Chen, X.H., Zu, S.H., Huang, W.L., & Gan, S.W., 2016, Simultaneous denoising and reconstruction of 5-D seismic data via damped rank-reduction method: Geophysical Journal International, 206, no. 3, 1695-1717.
https://doi.org/10.1093/gji/ggw230
Chen, Y.K., Zhou, Y.T., Chen, W., Zu, S.H., Huang, W.L., & Zhang, D., 2017, Empirical low-rank approximation for seismic noise attenuation: IEEE Transactions On Geoscience And Remote 28
Journal Pre-proof
Sensing, 55, no. 8, 4696–4711.
https://doi.org/10.1109/TGRS.2017.2698342
Chen, W., & Song, H., 2018, Automatic noise attenuation based on clustering and empirical wavelet transform: Journal of Applied Geophysics. 159, 649-665.
of
https://doi.org/10.1016/j.jappgeo.2018.09.025
ro
Dalai, B., Kumar, P., Yuan, X.H., 2019, De-noising receiver function data using the Seislet
-p
Transform: Geophysical Journal International, 217, no. 3, 2047-2055.
re
https://doi.org/10.1093/gji/ggz135
lP
Gaci, S., 2014, The use of wavelet-based denoising techniques to enhance the first-arrival picking on
na
seismic traces: IEEE Transactions On Geoscience And Remote Sensing, 52, no. 8, 4558–4563.
Jo ur
https://doi.org/10.1109/TGRS.2013.2282422
Gu, S.H., Zhang, L., Zuo, W.M., & Feng, X.C., 2014, Weighted Nuclear Norm Minimization with Application to Image Denoising: IEEE Conference on Computer Vision and Pattern Recognition(CVPR), 2862-2869.
https://doi.org/10.1109/CVPR.2014.366
Gorszczyk, A., Adamczyk, A.,Michal, M., 2014. Application of curvelet denoising to 2D and 3D seismic data — practical considerations: Journal of Applied Geophysics. 105, 78–94.
https://doi.org/10.1016/j.jappgeo.2014.03.009
29
Journal Pre-proof
Gulunay, N., 1986, FXDECON and complex Wiener prediction filter: 56th Annual International Meeting, SEG, Expanded Abstracts, 279–281.
Herrmann, F.J., & Hennenfent, G., 2008, Non-parametric seismic data recovery with curvelet frames: Geophysical Journal International, 173, no. 1, 233-248.
of
https://doi.org/10.1111/j.1365-246X.2007.03698.x
Kopsinis, Y., & McLaughlin, S., 2009, Development of EMD-based denoising methods inspired by
-p
ro
wavelet thresholding: IEEE Transactions on Signal Processing, 57, no. 4, 1351–1362.
re
https://doi.org/10.1109/TSP.2009.2013885
lP
Li, J., Wang, D.X., Ji, S., Li, Y., & Qian, Z.H., 2017, Seismic noise suppression using weighted
na
nuclear norm minimization method: Journal of Applied Geophysics, 146, 214–220.
Jo ur
https://doi.org/10.1016/j.jappgeo.2017.09.013
Li, G.H., Li, Y., & Yang, B.J., 2017, Seismic Exploration Random Noise on Land: Modeling and Application to Noise Suppression: IEEE Transactions On Geoscience And Remote Sensing, 55, no. 8, 4668-4681.
https://doi.org/10.1109/TGRS.2017.2697444
Liu, W., Cao, S.Y., & Wang, Z.M., 2017, Application of variational mode decomposition to seismic random noise reduction: Journal of Geophysics and Engineering, 14, no. 4, 888–899.
https://doi.org/10.1088/1742-2140/aa6b28
Liu, X.H., Tanaka, M., Okutomi, M., 2012, Noise level estimation using weak textured patches of a 30
Journal Pre-proof
single noisy image: 19th IEEE International Conference on Image Processing(ICIP), 665–668.
https://doi.org/10.1109/ICIP.2012.6466947
Liu, Y.P., Li, Y., Lin, H.B., & Ma, H.T., 2014, An Amplitude-Preserved Time–Frequency Peak Filtering Based on Empirical Mode Decomposition for Seismic Random Noise Reduction: IEEE
of
Geoscience and Remote Sensing Letters, 11, no. 5, 896-900.
ro
https://doi.org/10.1109/LGRS.2013.2281202
re
Journal of Applied Geophysics, 151, 246–262.
-p
Liu, Z., Y. Chen, and J. Ma, 2018, Ground roll attenuation by synchrosqueezed curvelet transform:
lP
https://doi.org/10.1016/j.jappgeo.2018.02.016
na
Lari, H.H., & Gholami, A., 2014. Curvelet-tv regularized bregman iteration for seismic random noise
Jo ur
attenuation: Journal of Applied Geophysics. 109, 233–241.
https://doi.org/10.1016/j.jappgeo.2014.08.005
Sacchi, M.D., Ulrych, T.J., & Walker, C.J., 1998, Interpolation and extrapolation using a high-resolution discrete Fourier transform: IEEE Transactions on Signal Processing, 46, no. 1, 31-38.
https://doi.org/10.1.1.154.4747
Souza, R., Lumley, D., & Shragge, J., 2016, Estimation of reservoir fluid saturation from 4D seismic data: Effects of noise on seismic amplitude and impedance attributes: Journal of Geophysics and Engineering, 14, no. 1, 51–68. 31
Journal Pre-proof
https://doi.org/10.1088/1742-2132/14/1/51
Shao, O.Y., Wang, L.L., Hu, X.Y., & Long, Z.D, 2019, Seismic random noise attenuation using nonlocal means via smooth patch ordering: Journal of Applied Geophysics, 167, 1-10.
https://doi.org/10.1016/j.jappgeo.2019.05.005
ro
denoising. Journal of Applied Geophysics. 69, 103–115.
of
Shan, H., Ma, J. & Yang, H., 2009. Comparisons of wavelets, contourlets and curvelets in seismic
-p
https://doi.org/10.1016/j.jappgeo.2009.08.002
re
Tang, N., Zhao, X., Li, Y, & Zhu, D, 2018, Adaptive threshold shearlet transform for surface
lP
microseismic data denoising: Journal of Applied Geophysics, 153, 64-74.
na
https://doi.org/10.1016/j.jappgeo.2018.03.019
Jo ur
Tang, G., & Ma, J.W., 2010, Application of total-variation-based curvelet shrinkage for three-dimensional seismic data denoising: IEEE Geoscience and Remote Sensing Letters, 8, no. 1, 103–107.
https://doi.org/10.1109/LGRS.2010.2052345
Wu, X.M., 2017, Structure-, stratigraphy- and fault-guided regularization in geophysical inversion: Geophysical Journal International, 210, no. 1, 184-195.
https://doi.org/10.1093/gji/ggx150
Wang, H, Cao, S.Y., Jiang, K.K., Wang, H, & Zhang, Q.C., 2019, Seismic data denoising for complex structure using BM3D and local similarity: Journal of Applied Geophysics, 170. 32
Journal Pre-proof
https://doi.org/10.1016/j.jappgeo.2019.04.018
Xue, Y., Man, M., Zu, S., Chang, F., Chen, Y., 2017. Amplitude-preserving iterative deblending of simultaneous source seismic data using high-order radon transform. Journal of Applied Geophysics. 139, 79–90.
of
https://doi.org/10.1016/j.jappgeo.2017.02.010
Yu, S.W., & Ma, J.W., 2018, Complex variational mode decomposition for sloppreserving denoising:
re
https://doi.org/10.1109/TGRS.2017.2751642
-p
ro
IEEE Transactions On Geoscience And Remote Sensing, 56, no. 1, 586–597.
lP
Yuan, S.Y., Liu, J.W., Wang, S.X., Wang, T.Y., & Shi, P.D., 2018, Seismic waveform classification
Letters, 15, no. 2, 272–276.
na
and first-break picking using convolution neural networks: IEEE Geoscience and Remote Sensing
Jo ur
https://doi.org/10.1109/LGRS.2017.2785834
Zhai, M.Y., 2014. Seismic data denoising based on the fractional Fourier transformation: Journal of Applied Geophysics. 109 , 62–70.
https://doi.org/10.1016/j.jappgeo.2014.07.012
Zhang, Z.J., & Klemperer, S. L., 2005, West-east variation in crustal thickness in northern Lhasa block central Tibet from deep seismic sounding data: Journal of Geophysical Research, 110.
https://doi.org/10.1029/2004JB003139
Zhang, Z.J., & Klemperer, S. L., 2010, Crustal structure of the Tethyan Himalaya, southern Tibet: 33
Journal Pre-proof
newconstraints from old wide-angle seismic data: Geophysical Journal International, 181, no. 3,
1247–1260.
https://doi.org/10.1111/j.1365-246x.2010.04578.x
Zhou, Y.T., & Zhu, Z.L., 2019, A hybrid method for noise suppression using variational mode
of
decomposition and singular spectrum analysis: Journal of Applied Geophysics, 161, 105-115.
Jo ur
na
lP
re
-p
ro
https://doi.org/10.1016/j.jappgeo.2018.10.025
34
Journal Pre-proof CRediT author statement
Juan Li: Validation, Conceptualization, Supervision. Wei Fan: Methodology, Software, Writing-Original Draft. Yue Li: Data Curation, Supervision, Project administration. Baojun Yang: Data Curation, Resources.
Jo ur
na
lP
re
-p
ro
of
Changgang Lu: Conceptualization, Writing-Review & Editing.
35