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
Abrahim B.A., Kadah Y. Speckle noise reduction method combining total variation and wavelet shrinkage for clinical ultrasound imaging. In: Proceedings of 2011 1st Middle East Conference on Biomedical Engineering. Sharjah, United Arab Emirates: IEEE, 2011: 80-83.
Bao Qianzong, Chen Wenchao, Gao Jinghuai. Seismic data random noise attenuation based on the second generation Curvelet transform[J]. Coal Geology & Exploration, 2010, 38(1): 66-70 (in Chinese with English abstract).
Candès E., Demanet L., Donoho D., Ying Lexing. Fast discrete Curvelet transforms[J]. Multiscale Modeling & Simulation, 2006, 5(3): 861-899.
Candès E.J., Donoho D.L. New tight frames of curvelets and optimal representations of objects with piecewise C2 singularities[J]. Communications on Pure and Applied Mathematics, 2004, 57(2): 219-266.
Chen Jianxiong, Wang Baoshan, Ge Hongkui, Lin Jianmin, Chen Yong. Using information of airgun to contract the low velocity zone in crustal in North China[J]. Earthquake Research in China, 2011, 27(1): 49-55 (in Chinese with English abstract).
Chen Yong, Zhu Rixiang. Proposed project of "underground bright lump"[J]. Advances in Earth Science, 2005, 20(5): 485-489 (in Chinese with English abstract).
Chen Yong, Wang Baoshan, Ge Honhkui, Xu Ping, Zhang Wei. Proposed of transmitted seismic stations[J]. Advances in Earth Science, 2007a, 22(5): 441-446 (in Chinese with English abstract).
Chen Yong, Zhang Xiankang, Qiu Xuelin, Ge Hongkui, Liu Baojin, Wang Baoshan. A new way to generate seismic waves for continental crustal exploration[J]. Chinese Science Bulletin, 2007b, 52(16): 2264-2268. DOI:10.1007/s11434-007-0247-4
Chen Yong, Liu Labo, Ge Hongkui, Liu Baojun, Qiu Xuelin. Using an airgun array in a land reservoir as the seismic source for seismotectonic studies in northern China:experiments and preliminary results[J]. Geophysical Prospecting, 2008, 56(4): 601-612. DOI:10.1111/j.1365-2478.2007.00679.x
Chen Yong, Wang Baoshan, Yao Huajian. Seismic airgun exploration of continental crust structures[J]. Science China Earth Sciences, 2017, 60(10): 1739-1751. DOI:10.1007/s11430-016-9096-6
Hu Jiupeng, Wang Baoshan, Chen Yong. The influence of water shape on land airgun triggered signals[J]. Journal of Seismological Research, 2017, 40(4): 543-549 (in Chinese with English abstract).
Jiang Shengmiao, Wang Baoshan, Zhang Yunpeng, Chen Yong. The influence of noise on the stacking effect of airgun signal and the automatic data screening method[J]. Journal of Seismological Research, 2017, 40(4): 534-542 (in Chinese with English abstract).
Jin Changkun, Jiang Yazhou, Zhang Jianzhong. An automatic picking method for slopes of locally coherent reflection events based on Curvelet transform[J]. Geophysical Prospecting for Petroleum, 2015, 54(4): 414-419 (in Chinese with English abstract).
Li Songlin. Using Ambient Noise Correlation Technique Based on the Time-frequency Phase-weighted Stacking to Study the Mantle Discontinuities' Structures in Northeast China[D]. Beijing: China University of Geosciences (Beijing), 2018 (in Chinese with English abstract).
Li Xiaobin, Chen Jia, Gao Qiong, Ye Beng, Long Zhiqiang, Zhang Yunpeng, Yang Jun. Estimating the receiving efficiency of stations by using the statistical characteristics of noise power spectral densities on active source signals[J]. Journal of Seismological Research, 2017, 40(4): 572-580 (in Chinese with English abstract).
Li Yonggang, Vidale J.E., Day S.M., Oglesby D.D., Cochran E. Postseismic fault healing on the rupture zone of the 1999 M 7.1 Hector mine, California, Earthquake[J]. Bulletin of the Seismological Society of America, 2003, 93(2): 854-869. DOI:10.1785/0120020131
Lin Chun, Wang Xuben. Application of second generation Curvelet transform in random noise decaying of seismic data[J]. Computer Engineering and Applications, 2009, 45(25): 222-224 (in Chinese with English abstract).
Liu Bideng, Li Xiaojun, Zhou Zhenghua, Liu Peixuan, Wang Yushi. Ground motion effect analysis of airgun source[J]. Acta Seismologica Sinica, 2011, 33(4): 539-544 (in Chinese).
Liu Wei, Cao Siyuan, Wang Zheng, Wang Yaliang, Dong Shuili. The self-adaptive random noise attenuation in curvelet domain based on Bayes estimation[J]. Geophysical Prospecting for Petroleum, 2013, 52(2): 115-120 (in Chinese).
Liu Huipu, Lv Jian, Pan Weijie, Liu Zhenghong. Research of an adaptive threshold value de-noising method based on the second Curvelet transform[J]. Microcomputer & Its Applications, 2014, 33(10): 76-79 (in Chinese with English abstract).
Liu Zifen, Su Youjin, Wang Baoshan, Wang Bin, Yang Jun, Li Xiaobin. Study on analysis method of travel time variations of seismic wave of active source in Binchuan[J]. Journal of Seismological Research, 2015, 38(4): 591-597 (in Chinese with English abstract).
Luan Yi, Yang Hongfeng, Wang Baoshan. Large volume air-gun waveform data processing (Ⅰ):Binchuan, Yunnan[J]. Earthquake Research in China, 2016, 32(2): 305-318 (in Chinese with English abstract).
Luo Guichun, Ge Hongkui, Wang Baoshan, Hu Ping, Chen Yong. Process in precise measurement of seismic velocity variation by correlated detection[J]. Progress in Geophysics, 2008a, 23(1): 56-62 (in Chinese with English abstract).
Luo Guichun, Ge Hongkui, Wang Baoshan, Hu Ping, Mu Hongwang, Chen Yong. Fired models of air-gun source and its application[J]. Earthquake Research in China, 2008b, 22(2): 112-120.
Qin Yi, Wang Jiaxu, Mao Yongfang. Signal denoising based on soft thresholding and reconstruction from dyadic wavelet transform modulus maxima[J]. Journal of Vibration, Measurement & Diagnosis, 2011, 31(5): 543-547 (in Chinese).
Qiu Aizhong. A novel denoising method with an improved wavelet thresholding[J]. Communications Technology, 2011, 44(8): 113-115 (in Chinese with English abstract).
Shaswary E. A new algorithm for time-delay estimation in ultrasonic echo signals and a comparative study of its performance evaluation in ultrasound elastography imaging[D]. Toronto, Canada: Ryerson University, 2010.
Shen Yang. The Research of De-noising Method Based on Curvelet Transform[D]. Harbin: Harbin Engineering University, 2011 (in Chinese with English abstract).
Su Jinbo, Wang Qiong, Wang Haitao, Shi Yongjun, Chen Xiangjun, Feng Lei, Chen Hao, Zhang Wenxiu. A research on large volum airgun source's signal reception in Hutubi Xinjiang using seismic network data[J]. Earthquake Research in China, 2016, 32(2): 202-208 (in Chinese with English abstract).
Su Jinbo, Wang Haitao, Wang Qiong, Zhang Wenlai, Ji Zhanbo, Wei Yunyun, Chen Hao, Zhang Wenxiu. A study on the repeatability of large volum arigun seismic source in Hutubi, Xinjiang using cluster analysis method[J]. Journal of Seismological Research, 2017, 40(4): 613-618 (in Chinese with English abstract).
Tang Jie, Wang Baoshan, Gehongkui, Chen Yong. Experiment and simulation of large capacity air-guns in deep structure exploration[J]. Earthquake Research in China, 2009, 23(4): 372-382.
Tong Zhongfei, Wang Deli, Liu Bing. Seismic data denoise based on Curvelet transform with the threshold method[J]. Journal of Jilin University (Earth Science Edition), 2008, 38(S1): 48-52 (in Chinese with English abstract).
Wang Baoshan, Wang Weitao, Ge Hongkui, Xu Ping, Wang Bin. Monitoring subsurface changes with active sources[J]. Advances in Earth Science, 2011, 26(3): 249-256 (in Chinese with English abstract).
Wang Baoshan, Ge Hongkui, Wang Bin, Wang Haitao, Zhang Yuansheng, Cai Huiteng, Chen Yong. Practices and advances in exploring the subsurface structure and its temporal evolution with repeatable artificial sources[J]. Earthquake Research in China, 2016, 32(2): 168-179 (in Chinese with English abstract).
Wang Bin, Yang Runhai, Wang Baoshan, Wang Weitao. The test study of the precise measurement of seismic wave travel time variety[J]. Journal of Yunnan University, 2012, 34(S2): 15-20 (in Chinese with English abstract).
Wang Bin, Li Xianbin, Liu Zifeng, Yang Jun, Yang Runhai, Wang Baoshan, Yang Wei. The source and observation system of Binchuan earthquake signal transmitting seismic station and its preliminary observation results[J]. Earthquake Research in China, 2016, 32(2): 193-201 (in Chinese with English abstract).
Wang Chenyi, Wang Jianjun. Information hiding method of second generation Curvelet transform[J]. Information and Electronic Engineering, 2008, 6(2): 105-110 (in Chinese with English abstract).
Wang Weitao, Wang Baoshan, Jiang Shengmiao, Hu Jiupeng, Zhang Yuansheng. A perspective review of seismological investigation on the crust at regional scale using the active airgun source[J]. Journal of Seismological Research, 2017, 40(4): 514-524 (in Chinese with English abstract).
Wang Weitao, Wang Baoshan, Jiang Shengmiao, Hu Jiupeng, Zhang Yuansheng. Regional scale seismological investigation on the continent crust using airgun sources——A perspective review[J]. Earthquake Research in China, 2018, 32(3): 315-329.
Wei Yunyun, Wang Haitao, Su Jinbo, Chen Xiangjun, Wang Qiong. The preliminary study on travel time abnormal variation of reflection wave phase of air-gun in Xinjiang before two earthquakes with MS5.0[J]. Earthquake Research in China, 2016, 32(2): 270-281 (in Chinese with English abstract).
Xia Ji. Study on the Characteristics of Large Volume Air-gun Source[D]. Harbin: Institute of Engineering Mechanics, China Earthquake Administration, 2017 (in Chinese with English abstract).
Xiang Ya, Wang Bin, Yang Runhai, Wang Weitao, Yang Haiyan. Comparison of propagation characteristics of signals triggered by air-gun active source and natural earthquakes[J]. Journal of Seismological Research, 2017, 40(4): 605-612 (in Chinese with English abstract).
Xu Hui, Liu Xuejun, Wang Bin, Wang Baoshan. Monitoring temporal variation of subsurface seismic velocity in Xiaojiang fault zone using direct wave cross correlation time delay technique with active source[J]. Journal of Seismological Research, 2015, 38(1): 7-15 (in Chinese with English abstract).
Xu Yihe, Wang Baoshan, Wang Weitao. Characteristics of air-gun signals excited in the Yangtze River from analysis of permanent stations' data[J]. Earthquake Research in China, 2016, 32(2): 282-294 (in Chinese with English abstract).
Yang Guoliang, Lei Songze. The image denoising method of soft and hard adaptive thresholding based on Curvelet transform and Bayesian estimation[J]. Journal of Xi'an Polytechnic University, 2011, 25(6): 857-861, 866 (in Chinese with English abstract).
Yang Kai, Liu Wei, Pan Yong. Random noise attenuation based on soft and hard threshold compromise in Curvelet domain[J]. Chinese Journal of Engineering Geophysics, 2013, 10(4): 437-441 (in Chinese with English abstract).
Yang Wei, Wang Baoshan, Liu Zhengyi, Yang Jun, Li Xiaobin, Chen Yong. Study on the source characteristic of downhole air-gun with different exciting environment[J]. Earthquake Research in China, 2016, 32(2): 231-240 (in Chinese with English abstract).
Yao Jiaqi, Cao Wenzhong, Ye Beng, Zhang Wei. Effects of nonlinear stacking and zero-phase processing on signal extraction from airgun source waveform[J]. Journal of Seismological Research, 2017, 40(4): 581-594 (in Chinese with English abstract).
You Xiuzhen, Li Jun, Lin Binhua, Huang Yandan, Wu Lihua, Guo Yang. Analysis of frequency domain water-level deconvolution method based on reservoir airgun source data[J]. Earthquake Research in China, 2018, 34(1): 14-24 (in Chinese with English abstract).
Zhai Qiushi, Yao Huajian, Wang Baoshan. Study on the deconvolution method and processing flow of air-gun source data[J]. Earthquake Research in China, 2016, 32(2): 295-304 (in Chinese with English abstract).
Zhang Yuansheng, Guo Xiao, Qin Manzhong, Liu Xuzhou, Wei Congxin, Shen Xuzhang, Yan Wenhua, Zou Rui. The construction of active source repeated monitoring in the Qilian Mountains of Gansu Province[J]. Earthquake Research in China, 2016a, 32(2): 209-215 (in Chinese with English abstract).
Zhang Yuansheng, Zou Rui, Yan Wenhua. The drift control technology of wire-suspension for air-gun source platform in the Qilian Mountain[J]. Earthquake Research in China, 2016b, 32(2): 438-442 (in Chinese with English abstract).
Zhang Yunpeng, Li Xiaobin, Wang Weitao, Wang Baoshan, Ye Beng, Yang Jun, Wang Bin. Data service system and data quality sssessment for mobile observation of transmitting seismic stations in Binchuan, Yunnan[J]. Journal of Seismological Research, 2017, 40(4): 525-533 (in Chinese with English abstract).
Zhao Qingmin, Ouyang Huan, Gu Daoping. Application study on Curvelet transform in face recognition[J]. Journal of Jiangxi Normal University (Natural Sciences Edition), 2013, 37(4): 397-400 (in Chinese with English abstract).
Zhuang Zhemin, Yao Weike, Yang Jinyao, Li Fenlan, Yuan Ye. Curvelet denoising algorithm for medical ultrasound image based on adaptive threshold[J]. Chinese Journal of Medical Instrumentation, 2014, 38(6): 398-401 (in Chinese with English abstract).