Earthquake Reaearch in China  2019, Vol. 33 Issue (2): 220-234     DOI: 10.19743/j.cnki.0891-4176.201902011
Weak Seismic Signal Extraction Based on the Curvelet Transform
TAN Junqing1, YANG Runhai2, WANG Bin2, XIANG Ya3
1. Yunnan University, Kunming 650091, China;
2. Yunnan Earthquake Agency, Kunming 650224, China;
3. Key Laboratory of Earthquake Geodesy, Institute of Seismology, China Earthquake Administration, Wuhan 430071, China
Abstract: Seismic signal denoising is a key step in seismic data processing. Airgun signals are easy to be interfered with by noise when it travels a long distance due to the weak energy of active source signal of the airgun. Aiming to solve this problem, and considering that the conventional Curvelet transform threshold processing method does not use the seismic spectrum information, we independently process the Curvelet scale layer corresponding to valid data based on the characteristics of the Curvelet transform of multi-scale, multi-direction and capable of expressing the sparse seismic signals in order to fully excavate the information features. Combined with the Curvelet adaptive threshold denoising the algorithm, we apply the Curvelet transform to denoising seismic signals while retaining the weak information in the signal as much as possible. The simulation experiments show that the improved threshold denoising method based on Curvelet transform is superior to the frequency domain filtering, wavelet denoising and traditional Curvelet denoising method in detailed information extraction and signal denoising of low SNR signals. The calculation accuracy of the relative wave velocity variation of underground medium is improved.
Key words: Seismic signal denoising     Airgun active source signal     Curvelet transform     The velocity of the underground medium

INTRODUCTION

Seismic waves are one of the few signals that can penetrate the whole earth. They are effective tools for studying the structure of the earth. People's understanding of the structural composition, state and evolution of the Earth is almost based on the analysis of seismic waves. The seismic waves are "illuminating the earth". A bright light (Chen Yong et al., 2005; Chen Yong et al., 2017), which carries abundant underground media structure information and physical information, is an important basis for us to study the internal structure of the earth (Wang Bin et al., 2012).

Earthquakes are caused by the instability of active faults under crustal stress. In recent years, the use of seismology to study the physical properties of underground media has received more and more attention. Although the stress state can also be obtained by studying the attenuation and anisotropy of seismic waves, the accurate measurement of the seismic wave velocity change is a more effective method (Chen Yong et al., 2007a, 2007b; Xu Hui et al., 2015; Zhang Yunpeng et al., 2017; Li Yonggang et al., 2003; Luan Yi et al., 2016; Luo Guichun et al., 2008a). The determination of seismic wave velocity changes requires the use of repetitive seismic sources that produce similar waveforms. The difference between the seismic waveforms of repetitive earthquakes with the same source location and propagation path recorded by the same station reflects the changes in the physical properties of the propagation path medium (Wang Baoshan et al., 2011; Wang Wenzhi et al., 2015; Xiang Ya et al., 2017; Chen Jianxiong et al., 2011; Wei Yunyun et al., 2016). The limitations of natural earthquakes in time and space distribution and in number limit their application in local crustal exploration. However, the time and location of an artificial source are known, and the seismic waves generated by it can be used to achieve high-precision measurement of subsurface medium changes (Chen Yong et al., 2007a, 2007b; Chen Yong et al., 2008; Wang Baoshan et al., 2016; Hu Jiupeng et al., 2017; Luo Guichun et al., 2008b; Tang Jie et al., 2009; Wang Weitao et al., 2018).

In the past, the most used artificial sources to generate seismic waves were blasting sources. With an increasing emphasis put on environmental protection, the use of blasting sources is increasingly restricted (Wang Weitao et al., 2017). The airgun source is a new type of source suitable for studying underground media structures by instantaneously releasing high-pressure air under water to excite seismic waves (Su Jinbo et al., 2017; Chen Yong et al., 2007a, 2007b; Xu Yihe et al., 2016). At present, Binchuan Yunnan, Hutubi Xinjiang, Qilianshan Gansu airgun active source signal transmitting stations are established to conduct experiments and research on underground structures (Li Xiaobin et al., 2017; Wang Bin et al., 2016; Zhang Yuansheng et al., 2016a; Zhang Yuansheng et al., 2016b; Su Jinbo et al., 2016). In active source airgun signal recording, with the increase of the distance of the receiving station, the scattering and attenuation of the signal gradually submerges the seismic signal in the ambient noise, so we need to extract the effective seismic signal in low SNR seismic records (Liu Bideng et al., 2011; Xu Yihe et al., 2016; Jiang Shengmiao et al., 2017; You Xiuzhen et al., 2018; Yao Jiaqi et al., 2017).

Common methods for improving the signal strength in seismic records are: frequency domain filtering, median filtering, superposition techniques, etc. (Yang Wei et al., 2016; Xia Ji, 2017; Li Songlin, 2018). Later, seismologists proposed the use of the deconvolution method to obtain the approximate Green's function between the airgun source and the receiving station to eliminate the influence of incomplete repetition of the source (Zhai Qiushi et al., 2016; Liu Zifeng et al., 2015). In recent years, two-dimensional wavelet transform has been widely used in seismic data processing and image processing. When processing image data, it can only express the cross-edge characteristics of the wavelet coefficients, but cannot express the along-the-edge characteristics, so the ability to process the image edge in detail is general. Donoho D. L., Johnstione I. M. et al. proposed the first generation Curvelet transform based on the Rideglet transform in 1999. Due to the complexity of the algorithm and the redundancy of data, they proposed the second-generation Curvelet transform, making analysis directly from the frequency domain. This method is simple and easy to understand, and compared to wavelet transform, its strongest feature is that it is a highly anisotropic, bandpassing, multiresolutional, local function analysis method (Candès E. J. et al., 2004; Candès E. et al., 2006; Zhao Qingmin et al., 2013; Qiu Aizhong, 2011; Lin Chun et al., 2009; Zhuang Zhemin et al., 2014; Jin Changkun et al., 2015; Wang Chenyi et al., 2008).

In view of the fact that the general denoising method is not ideal for the processing of low SNR data, Curvelet transform is a signal analysis method with high directionality and time-frequency resolution, which is suitable for the processing and analysis of seismic signals. This paper is based on the Curvelet adaptive threshold denoising (Liu Huipu et al., 2014; Yang Guoliang et al., 2011; Shen Yang, 2011; Liu Wei et al., 2013), improving the threshold determination and processing of Curvelet coefficients, through the active source of airguns. The processing and analysis of analog signals and actual data demonstrate the effectiveness and feasibility of the improved data processing method based on Curvelet in low SNR signal processing.

1 PRINCIPLE OF CURVELET TRANSFORM

Curvelet transform, like wavelet transform and ridgelet transform, belongs to the category of coefficient theory. The sparse representation of signals is realized by the inner product of basis function and signal (Wang Chenyi et al., 2008; Candès E. J. et al., 2004). In two-dimensional continuous space R2, r, θ is the polar coordinate in frequency domain, ω the frequency domain variable and x the spatial position variable. The "radius window" W(r) in the frequency domain, and "angle window" V(t) of the smooth, non-negative, real value in the frequency domain is defined to satisfy

 $\left\{\begin{array}{c}{\sum_{j=-\infty}^{\infty} W^{2}\left(2^{j} r\right)=1, r \in\left(\frac{3}{4}, \frac{3}{2}\right)} \\ {\sum_{l=-\infty}^{\infty} V^{2}(t-l)=1, t \in\left(-\frac{1}{2}, \frac{1}{2}\right)}\end{array}\right.$ (1)

where, j is the radius, l the angle parameter variable.

For all the radi of j≥1, the "frequency domain window" in the Fourier frequency domain is defined as

 ${U_j}(r, \theta) = {2^{ - \frac{{3j}}{4}}}W\left({{2^{ - jr}}} \right)V\left({\frac{{{2^{\left\lfloor {\frac{j}{2}} \right\rfloor \theta }}}}{{2{\rm{ \mathtt{ π} }}}}} \right)$ (2)

where, ${\left\lfloor {\frac{j}{2}} \right\rfloor }$ represents the integer part of ${\frac{j}{2}}$, the support area of Uj in the polar coordinate of the frequency domain is a "wedge window" in the polar coordinate of frequency domain. As shown in Fig. 1, at the jth level and kth angle in the frequency domain, the Curvelet transform coefficients at position k=(k1, K2)∈Z2 are defined as:

 Fig. 1 (a) Continuous Curvelet frequency domain block diagram (Emmanuel Candes et al.) and (b) discrete Curvelet frequency domain block diagram
 $C(i, j, k)\frac{1}{{{{(2\pi)}^2}}}\int {\hat f} (\omega){U_j}\left({{R_\theta }\omega } \right){e^{i < x_k^{(j, l)}}}, \omega > {\rm{d}}\omega$ (3)

In a continuous Curvelet transform, the frequency domain window function is able to filter out concentric circles corresponding to different frequencies (Fig. 1(a)). In the discrete Curvelet transform, the concentric square region is used instead (Fig. 1(b)).

The local frequency domain window in Cartesian coordinates is defined as (Candès E. J. et al., 2004; Liu Huipu et al., 2014; Wang Chenyi et al., 2008; Tong Zhongfei et al., 2008):

 ${{\tilde U}_{j, l}}(\omega) = {{\tilde W}_j}(\omega){V_j}\left({{S_{\theta l}}\omega } \right)$ (4)
 $\left\{ {\begin{array}{*{20}{l}} {{W_j}(\omega) = \sqrt {\left({\Phi _{j + 1}^2} \right.} - \Phi _j^2(\omega)), j \ge 0}\\ {{V_j}(\omega) = V\left({{2^{j/2}}{\omega _1}/{\omega _2}} \right)} \end{array}} \right.$ (5)

Where, Φ is defined as the inner product of two one-dimensional low-pass windows.

The discrete Curvelet transform is defined as:

 $\begin{array}{l} c(j, l, k) = \int {\hat f} (\omega){{\tilde U}_j}\left({S_{\theta l}^T\omega } \right){e^{i\left({{S^{ - T}}{\theta _l}b, \omega } \right)}}{\rm{d}}\omega \\ = \int {\hat f} \left({{S_{{\theta _l}}}\omega } \right){{\tilde U}_j}(\omega){e^{i(b, \omega)}}{\rm{d}}\omega \end{array}$ (6)

There are two kinds of discrete Curvelet transform algorithms. This paper uses Wrapping-based fast Curvelet transform algorithm, which has fast calculation speed and high algorithm efficiency (Zhao Qingmin et al., 2013; Shen Yang, 2011; Candès E. J. et al., 2004).

2 PRINCIPLE OF CURVELET DENOISING 2.1 Curvelet Transform Denoising Principle

The Curvelet transform can represent a signal as a curvelet coefficient of different scales, different directions, and different positions, and is a linear transformation. For seismic signals, there are:

 $\mathrm{d}(x, t)=s(x, t)+n(x, t)$ (7)

Where d is seismic data with random noise, s is the signal, n is random noise, and x and t are spatial and temporal coordinates.

Curvelet transform coefficients refer to the sum of signal curve transform coefficients and noise curve transform coefficients, namely,

 ${C_{\rm{d}}}\left({j, l, {k_x}, {k_t}} \right) = {C_s}\left({j, l, {k_x}, {k_t}} \right) + {C_n}\left({j, l, {k_x}, {k_t}} \right)$ (8)

The Curvelet transform coefficients of d, s, n are Cd, Cs and Cn, which represent scale factor, angle factor and translation factor respectively (Bao Qianzong et al., 2010).

The seismic signals have strong correlations at different angles in the space-time domain. The Curvelet transform coefficients are distributed in the finite region of the Curvelet domain, and the transform coefficients have strong amplitudes. The random noise has no correlation in the space-time domain, and the Curvelet coefficients are distributed throughout the coefficient domain. In the domain, the transform coefficients have weak amplitudes. Therefore, the interference signal can be removed by using the difference between the amplitude of the signal and the noise (Shaswary E., 2010; Liu Huipu et al., 2014).

2.2 Target Scale Layer Corresponding to Curvelet Coefficient Processing

The seismic data contains many small fluctuations and detailed information. The related information mainly exists in several scale layers in the Curvelet domain. When using the general Curvelet denoising process to process seismic signals, the data signal-to-noise ratio can be improved, but it is also very much likely to remove the effective signal as noise interference.

The frequency range of the Curvelet decomposition is related to the sampling frequency. If the signal of the sampling rate Fs is decomposed by layer N, the frequency bands are approximated to $\frac{\mathrm{Fs} / 2}{2^{N}}$. For the seismic signal, when the sampling frequency is 100Hz, according to the Nyquist theorem, the frequency band corresponding to the layer Fine at layer N is 25Hz to 50Hz; for the layer N-1, the corresponding frequency band is 12.5Hz to 25Hz. The layer N-2 corresponds to 6.25Hz to 12.5Hz, and so on. In this paper, we refer to the scale layer corresponding to the main frequency band of the seismic signal as the target scale layer. By combining the Curvelet adaptive threshold denoising algorithm, the target scale layer of the seismic signal's Curvelet domain is processed separately, that is, the coefficients in the target scale layer are taken as absolute values and arranged in order of magnitude, and parts of the coefficients are taken in appropriate proportion to the direction of large amplitude, and the remaining data is removed. This way, most of the noise of the target scale layer is removed while preserving the effective data information as much as possible and reducing the loss of weak information of the seismic signal.

2.3 Curvelet Threshold Estimation

The definition of multi-scale unified threshold estimation proposed by Donoho D.L. is as follows (Shen Yang, 2011; Liu Huipu et al., 2014):

 $T=\sigma_{n} \sqrt{2 \operatorname{Ln} N}$ (9)

Among them, N is the Fine scale layer coefficient, the noise variance is estimated by the robust median proposed by Donoho D. L., namely:

 $\sigma_{n}=\frac{\operatorname{Median}\left(\left|C_{j_{0}, l_{0}}\right|\right)}{0.6745}$ (10)

where Cj0, l0 is the coefficient of the scale layer Fine.

The simulation experiment found that the standard deviation of noise estimated by formula (10) in the Curvelet domain has a certain deviation from the theoretical value. Through multiple simulation experiments, the optimized robust median estimation formula is as follows:

 $\hat{\sigma}_{n}=\frac{\operatorname{Median}\left(\left|C_{j, i}\right|\right)}{0.5843}$ (11)

Fig. 2 is the Gaussian white noise with a mean of 0 added to the analog signal. The Curvelet transform is performed. The standard deviation of the noise signal can be calculated by formula (10) and (11) respectively. The estimation accuracy is higher by the optimized definition formula in the Curvelet domain.

 Fig. 2 The comparison of robust median optimization algorithm with original algorithm

The coefficients of the second-generation Curvelet transform are distributed in different scales and directions, so using the same threshold in all scales will inevitably have greater impact on the noise reduction results. Multi-scale and multi-directional properties of the Curvelet transform can be exploited by using different thresholds at different scales and directions (Liu Huipu et al., 2014; Abrahim B. A. et al., 2011).

In the determination of thresholds of different scales and directions in the Curvelet domain, we need to use the sub-band adaptive threshold calculation method based on the Bayes criterion to minimize the Bayes risk and take different values between different Curvelet coefficient layers. So, the threshold calculation based on the minimization conditions of Bayes risk is too complicated, and its approximate definition formula is adopted (Yang Guoliang et al., 2011; Liu Wei et al., 2013) as,

 $T_{j, l}=\frac{\hat{\sigma}_{n}^{2}(j, l)}{\hat{\sigma}_{g}(j, l)}$ (12)

Where Tj, l is the threshold in the l direction at the j scale, ${\hat \sigma _n^2(j, l)}$ is the estimated value of the noise variance in the l direction at the j scale; ${{{\hat \sigma }_g}}$ is the estimated value of the ideal signal standard deviation in the l direction at the scale.

The optimized robust median method is used to estimate the noise standard deviation of the noised signal:

 $\hat{\sigma}_{n}=\frac{\operatorname{Median}\left(\left|C_{j, l}\right|\right)}{0.5843}$ (13)

Where Cj, l is the coefficient in the l direction at j scale.

The BayesShrink threshold is determined by the noise variance $\sigma _v^2(j, l)$ and the standard deviation ${\sigma _g}(j, l)$ of the ideal signal, and because the noisy signal $\sigma _g^2(j, l)$, the ideal noise-free signal $\sigma _v^2(j, l)$, and the noise signal $\sigma _v^2(j, l)$ are independent of each other, the relationship between their variance of the sub-band after the Curvelet transform is as follows:

 $\sigma_{f}^{2}(j, l)=\sigma_{g}^{2}(j, l)+\sigma_{n}^{2}(j, l)$ (14)

Further calculation of the sub-band variance can be obtained as

 $\sigma _f^2(j, l) = \frac{1}{{mn}}\sum\limits_{k = 1}^m {\sum\limits_{i = 1}^n {C_{j, l}^2} }$ (15)

Among them, mn is the size of each sub-band.

From formulae (13), (14), (15), the estimated variance of each sub-band can be calculated as

 $\hat{\sigma}_{g}(j, l)=\sqrt{\max \left(\hat{\sigma}_{f}^{2}(j, l)-\hat{\sigma}_{n}^{2}(j, l), 0\right)}$ (16)

Finally, the noise factor of the Curvelet can be removed by combining appropriate threshold processing methods. The compromise threshold method combines the advantages of soft and hard threshold calculations, and has a small amount of calculation. It has achieved good results in seismic denoising. This paper takes this method as an example to carry out experimental calculations (Shen Yang, 2011; Bao Qianzong et al., 2010; Yang Kai et al., 2013; Qin Yi et al., 2011).

3 EXPERIMENTAL SIMULATION AND RESULT ANALYSIS 3.1 Analog Signal

The simulation data (Fig. 3) simulates the excitation process of the airgun active source signal for three years. There are a total (365×3) of trace data, with a sampling rate of 100 and a sampling duration of 14s for each trace signal, including the information of monthly and annual variation of the velocity of the underground medium, as shown in Fig. 3(b). White Gaussian noise is added to the signal. This paper takes the -5dB SNR data as an example for the illustration. The noised data is shown in Figs. 3(c) and (d).

 Fig. 3 The simulated signal of airgun active source repeat excitation
3.2 Calculation and Evaluation Criteria

The improved data processing flow is shown in Fig. 4.

 Fig. 4 The flowchart of Curvelet denoising

In threshold processing, the Coarse scale layer corresponds to the large-scale contour information of the image, and all of them are retained, while the Detail and Finest scale layer data outside the target scale layer are processed and calculated by the Curvelet adaptive compromise threshold method.

The evaluation of the data results is an important step in judging the denoising effect. This paper comprehensively considers various factors to define the evaluation criteria as follows:

 $Eva = \frac{{r \cdot {\mathop{\rm SNR}\nolimits} }}{{|\mathit{Res}|}}$ (17)

Where, r is the correlation coefficient between the deceleration wave velocity curve and the theoretical signal wave velocity curve, 0 < r≤1.

Signal to noise ratio:

 $\mathrm{SNR}=\frac{P_{s}}{P_{n}}$ (18)

In equation (18), Ps is the two-dimensional denoising signal power, and Pn is the two-dimensional standard signal power.

Residual:

 $Res = \sum\limits_{i = 1}^n {\left| {{v_i} - {{\hat v}_i}} \right|}$ (19)

Equation (19) calculates the residual of the relative curve of the wave velocity, the closer the residual is to 0, the closer the calculated wave velocity curve is to the theoretical value.

3.3 Frequency Domain Filtering and Wavelet Denoising Calculation Results

Fig. 5 and Fig. 6 show the processing results of the noisy data by frequency domain filtering and wavelet hard threshold denoising. The signal-to-noise ratio of the data after frequency domain filtering is 5.68dB and the Eva value is 4.28. The correlation coefficient between the calculated and the theoretical wave velocity change rate is 0.9395; the signal-to-noise ratio of the data after wavelet processing is 11.33dB, the Eva value is 7.17, and the correlation coefficient of the wave velocity curve is 0.9359. Obviously, the correlation coefficient and Eva value of the rate of change of the wave velocity calculated by the two methods are relatively low. Observing Fig. 5, it can be seen that there is still more noise remaining after frequency domain filtering, and less noise residual after wavelet processing. In the wave velocity change rate curve (Fig. 6), the wave velocity change rate curve calculated after frequency domain filtering is more chaotic, while the monthly change information by wavelet hard threshold denoising is not effectively represented due to the over-processing of data.

 Fig. 5 The results of frequency domain filtering and wavelet denoising

 Fig. 6 The comparison of relative wave velocity curves
3.4 Curvelet Adaptive Denoising and Improved Algorithm

Fig. 7 shows the data processing results using the Curvelet adaptive threshold denoising and improved Curvelet methods.

 Fig. 7 Data processing results using Curvelet adaptive threshold denoising and the improved denoising methods

After Curvelet adaptive threshold denoising, the signal-to-noise ratio of the data is 19.41dB, the correlation coefficient of the calculated relative variation curve of wave velocity is 0.9277 and the Eva value is 12.88. The improved method calculates the signal-to-noise ratio of 17.53dB, the correlation coefficient is 0.9898 and the Eva value is 33.73. The signal-to-noise ratio of the data denoised by the two Curvelet methods is obviously improved, while the Eva value is greatly increased by the improved method. The morphology of the data is close to the theoretical value in the single-trace waveform and spectrogram, but in the relative variation of wave velocity (Fig. 8), the monthly variation information in the wave velocity curve is lost by Curvelet adaptive threshold denoising, which can only express the annual variation, while the improved Curvelet threshold processing method can recover the information of annual and monthly variations, and the correlation between the wave velocity variation curve and the theoretical one is high.

 Fig. 8 Relative curve of wave velocity calculated by Curvelet denoising data
3.5 Comparison of different SNR data

The Eva values of the data with different SNR are compared by different methods (Fig. 9).From Fig. 9, it can be found that at different SNR levels, the restoration effect of the Curvelet algorithm is better than that of wavelet denoising and frequency domain filtering, and the effect of wavelet denoising and frequency domain filtering is comparable. The improved method has a higher Eva value when SNR is lower than 5dB. The experimental results show that the improved Curvelet denoising algorithm can better restore the monthly variation of analog signals with SNR of above -8~5dB. While the unimproved Curvelet threshold denoising method can only restore the monthly variation information for signals with SNR above 3dB. Considering that the object of this method is the signals with SNR of less than 0dB, the improved method is an optimal method for extracting weak signals.

 Fig. 9 The comparison of Eva by different methods
4 ACTUAL DATA PROCESSING

The following is a description of the low signal-to-noise ratio airgun signal recorded by Station 53036 about 67km from the airgun source in Yunnan Province. The data recording time is from January 4, 2016 to February 22, 2016. The airgun source was fired 10-15 times per day, the excitation interval is about half an hour, the signal sampling rate is 100, and the first 60s of the recorded signal is calculated and analyzed. This signal has been subjected to the frequency domain filtering of 2-7Hz before applying the method. The processing results are shown in Fig. 10.

 Fig. 10 The comparison of the signals before and after data processing using the improved Curvelet threshold method

It can be seen from Fig. 10 that a large amount of noise remains in the data processed by the frequency domain filtering (Fig. 10(a), (c)); The result of further processing of this signal with the method proposed in this paper is shown in Fig. 10(b), (d), the signal-to-noise ratio of the signal is obviously improved. The relative curve of wave velocity for each of the signals in Fig. 10(a)(b) are calculated, the results are shown in Fig. 11.

 Fig. 11 The comparison of relative curve of wave velocity before and after the Curvelet processing

It can be seen from Fig. 11(a) that the relative change curve of the wave velocity calculated with the frequency domain filtered data is more discrete and irregular. The curve of the wave velocity calculated with the data processed by the modified Curvelet algorithm (Fig. 11(b)) is easy to identify and analyze, and also proves the ability and advantages of the improved method of this paper in actual data processing.

5 CONCLUSION

Focused on the difficulty in removing noise from low signal-to-noise ratio seismic signals, this paper studies the improved Curvelet threshold processing method for low signal-to-noise ratio signals of active sources, and compares it with other methods through experiments.

(1) Generally speaking, two-dimensional denoising is better than one-dimensional denoising, because one-dimensional denoising only considers single-trace signals, while two-dimensional denoising optimizes the processing from a global perspective, which overcomes the problem that one-dimensional denoising cannot distinguish signals and noises in the corresponding frequency band.

(2) When using a robust median to estimate noise variance, it is found that there is a certain deviation between the estimated noise variance and the real value directly in the Curvelet domain. After many experiments and analysis, this paper formulates an improved robust median formula for the Curvelet domain, which improves the accuracy of the Curvelet threshold estimation.

(3) The improved threshold denoising method based on the Curvelet transform proposed in this paper solves the over-processing problem of Curvelet denoising algorithm in low signal-to-noise ratio (SNR) signals by combining Curvelet adaptive threshold algorithm with special processing of target scale layer coefficients. The simulation experiment shows that this method has better ability of denoising and detail recovery for low SNR analog airgun signal. While removing noise interference, it enhances the detail recovery ability of low SNR (less than 0dB) signal.

(4) The application of the improved method to the actual data processing results show that the denoising data retains the effective signal as much as possible, while the noise has been well removed, the calculation accuracy of velocity change rate is improved significantly, and the denoising effect is obvious. Of course, this method still has many details to think about, which requires further study and improvement in the future.

ACKNOWLEDGEMENT

The author is grateful to the reviewers for their valuable comments on the content of this article.

REFERENCES