Document Type : Research Paper
Authors
 Saman Gholtashi ^{1}
 Mohammad Amir Nazari Siahsar ^{2}
 Amin RoshandelKahoo ^{} ^{1}
 Hosein Marvi ^{2}
 Alireza Ahmadifard ^{2}
^{1} Department of Mining, Petroleum & Geophysics Engineering, University of Shahrood, Shahrood, Iran
^{2} Department of Electrical & Robotic Engineering, University of Shahrood, Shahrood, Iran
Abstract
Seismic waves are nonstationary due to its propagation through the earth. Timefrequency transforms are suitable tools for analyzing nonstationary seismic signals. Spectral decomposition can reveal the nonstationary characteristics which cannot be easily observed in the time or frequency representation alone. Various types of spectral decomposition methods have been introduced by some researchers. Conventional spectral decompositions have some restrictions such as Heisenberg uncertainty principle and crossterms which limit their applications in signal analysis. In this paper, synchrosqueezingbased transforms were used to overcome the mentioned restrictions; also, as an application of this new high resolution timefrequency analysis method, it was applied to random noise removal and the detection of lowfrequency shadows in seismic data. The efficiency of this method is evaluated by applying it to both synthetic and real seismic data. The results show that the mentioned transform is a proper tool for seismic data processing and interpretation.
Keywords
1. Introduction
Nowadays, timeseries analyses have vast applications in seismic data processing and interpretation (Castagna et al., 2003; Leite et al., 2008; Martelet et al., 2001; Sinha et al., 2005). By considering the low pass filtering behavior of the earth, the frequency content of seismic waves changes due to propagating through the earth (Roshandel Kahoo and Nejati Kalateh, 2011). Because of the nonstationary behavior of seismic signals, it is necessary to derive their time and frequency information simultaneously for many applications such as deconvolution, noise reduction, hydrocarbon reservoirs detection, and seismic attributes calculation. Although conventional methods in time and frequency domains have a wide application in signal processing, they cannot represent time and frequency information simultaneously. New eras of signal processing were created by introducing timefrequency representation (TFR) of signals, which increased the efficiency of signal processing. Hitherto, many approaches such as continuous wavelet transform (CWT) (Mallat, 1999), shorttime Fourier transform (STFT) (Gabor, 1946), WignerVille distribution (WVD) (Wigner, 1932), and S transform (Stockwell et al., 1996) have been proposed for timefrequency (TF) analyses. The timefrequency analysis, as a powerful and popular tool for the analysis of seismic data has widely been used in seismic data processing and interpretation. Chakraborty and Okaya (1995) suggested improved processing and interpretation algorithms for seismic data using various methods of spectral decomposition. Partyka et al. (1999) applied STFT to 3D seismic data to quantify thinbed interference and detect subtle discontinuities. The Stransform introduced by Stockwell et al. (1996) as a timefrequency analysis technique, which combines both CWT and STFT features, has widely been applied to seismic data processing (Askari and Siahkoohi, 2008; Ditommaso et al., 2010; Gao et al., 2003). Leite et al. (2008) used the wavelet transform to illuminate the coherent noise in reflection seismic data. Lin et al. (2013) and Deng et al. (2015) applied various versions of timefrequency peak filtering to seismic random noise attenuation.
Each of the aforementioned TFR’s has exclusive limitations, which affect their applications. They arise from various reasons such as Heisenberg Uncertainty Principle, quadratic superposition principle, and so on. However, these transforms have some advantages, which cannot be ignored. Looking for possible mathematical tools for preserving the existing advantages of the conventional methods and eliminating their disadvantages is the motivation of various studies. One of the improved TFR is the datadriven timefrequency analysis of multivariate signals, which is achieved through the empirical mode decomposition (EMD) algorithm (Huang et al., 1998), its multivariate extensions, the ensemble EMD (Wu and Huang, 2009), and its complete ensemble EMD (Torres et al., 2011).
Recently Daubechies et al.(2011) proposed synchrosqueezing transform by relying on wavelet (SSWT) in the context of speech signal processing by getting the philosophy of the EMD approach, but with a strong mathematical based theory. In comparison to other common TFR methods, SSWT, due to its high resolution in time and frequency, became a useful tool for analyzing nonstationary signals (Li and Liang, 2012; Marsset et al., 2014; Wang et al., 2014; Yangkang et al., 2014). Herrera et al. (2014) used the SSWT for analyzing the seismic signal and produced promising results on synthetic and field data examples. Herrera et al. (2015) used the SSWT to separate the body wave in the earthquake seismic waves. Oberlin et al. (2014) adapted the formulation of synchrosqueezing to the STFT and introduced the Fourier based synchrosqueezing (FSST).
In this paper, at first, synchrosqueezingbased method is briefly introduced. Then, its efficiency is evaluated in the detection of lowfrequency shadows. Finally, a new algorithm is presented for the random noise attenuation of seismic data based on the sparsity property of synchrosqueezingbased transform (SSBT).
2. Synchrosqueezingbased transform (SSBT)
SSWT as a sparse representation was introduced by Daubechies et al. (2011). This transform is a TFR based on wavelet transform which similar to EMD can decompose a signal into intrinsic mode functions (Herrera et al., 2014). Unlike EMD, SSWT has a firm theoretical foundation (Thakur et al., 2013; Wu et al., 2011). This transform is an adaptive and invertible tool which enhances the resolution of TFR.
The SSWT initially introduced in speech signal context (Daubechies et al., 2011; Daubechies and Maes, 1996). This transform decompose into intrinsic modes as given in Equation 1:
(1) 
where, is instantaneous amplitude; represents noise or error; is the number of the modes and stands for instantaneous phase. The instantaneous frequency of each mode is calculated by the deviation of instantaneous phase (Boashash, 2003) as reads:
(2) 
Transforms like STFT and continuous wavelet transform (CWT) spread signal energy around the instantaneous frequency of the original signal (Daubechies and Maes, 1996). The next part explains how synchrosqueezing wavelet transform is extracted from wavelet transform. The wellknown continuous wavelet transform of signal is expressed by Equation 3 (Mallat, 2008):
(3) 
where, is the mother wavelet and denotes the complex conjugate of ; is a scale parameter (Mallat, 2008). Wavelet transform spreads the energy of the signal around the scale axis, and it causes the representation to be blurred. If the energy smearing around the time axis is negligible (Daubechies and Maes, 1996), then the instantaneous frequency , around which the energy must be concentrated , can be derived by Equation 4 (Boashash, 2003):
(4) 
The synchrosqueezed transform can be computed only in the centers ( ) of frequency bins , with as given by Equation 5 (Daubechies et al., 2011).
(5) 
where, is the discrete variable of scale and . The SSBT is invertible and the original signal can be obtained by Equation 6.
(6) 
where, is a constant which depends on the mother wavelet in the continuous wavelet transform (Chen et al., 2014; Wang et al., 2014).
Oberlin et al. (2014) adapted the same formulation and theoretical foundation of SSWT to the STFT and introduced the Fourierbased synchrosqueezing (FSST) which has a more concentrated representation than STFT. In the next section, the performance of SSWT and FSST for a complex signal is shown.
3. Synthetic example
For the sake of clarity, a comparison between the timefrequency representation of SSWT, FSST with wavelet transform, and STFT for a complex signal depicted in Figure 1 with a sampling rate of 2 ms and 500 samples is shown in Figure 2. This signal consists of three sinusoidal components with a frequency of 10, 50, and 70 Hz, one Morlet wavelet with a central frequency of 25 Hz (.25 s), and three Ricker components with a central frequency of 30, 70, and 90 Hz respectively at 0.59, 0.77, and 0.94 s.
Figure 1
Signal and its Fourier spectrum; a) the original signal; b) the Fourier spectrum.
Figure 2
Results of timefrequency representation of the synthetic signal of Figure 1 obtained from different methods; (a) SSWT, (b) FSST, (c) Wavelet transform, and (d) STFT.
As can be seen in Figure 2, the resolution of STFT and wavelet transforms increases by synchrosqueezing. In the case of SSWT, Morlet wavelet is chosen as the mother wavelet and the number of voices per octave is 32; the Gaussian window is used for the FSST. The runtimes of SSWT, FSST, STFT, and CWT for this example are 0.9 s, 0.4 s, 0.2 s, and 0.1 s respectively. The performance of the reconstruction is shown in Figure 3. It is obvious that there are no major differences between the original and the reconstructed signals. The reconstructed error is indicated with the red line.
Figure 3
Reconstruction and error; a) the reconstructed signal with synchrosqueezing wavelet; b) the reconstructed signal with synchrosqueezing STFT; the red line indicates the reconstruction error.
Figure 4
(a) 2D real seismic data from an Iran hydrocarbon field and (b) its amplitude spectrum.
4. Lowfrequency shadow
Using spectral decomposition to detect the lowfrequency shadows associated with hydrocarbons is a common application of TFR (Chabyshova and Goloshubin, 2014; He et al., 2008). In fact, these shadows are often related to an additional energy occurring at low frequencies rather than the preferential attenuation of higher frequencies. One possible explanation is that these are locally converted shear waves which have traveled mostly as Pwaves and thus arrive slightly after the true primary event (Castagna et al., 2003). A lowfrequency shadow is a zone in seismic data characterized by anomalously low frequencies, which occurs beneath the gas reservoirs. One may identify the seismic lowfrequency shadows by comparing different common frequency slices. The existence of lowfrequency anomaly in common frequency slices rather than medium to high frequency slices is the indicator of lowfrequency shadow.
Here, a real example from one of Iran hydrocarbon fields is selected to indicate the performance of the SSTB and its ability to detect lowfrequency shadow. Since the output of SSWT and FSST is approximately similar and FSST had never been used for lowshadow frequency, FSST is used herein. Figure 4 shows the 2D real seismic section and its average amplitude spectrum. By considering the amplitude spectrum of data (Figure 4b), the 15 Hz and 55 Hz frequency slices were chosen as the low and high frequencies respectively. The 3D cubes of timefrequency transform of 2D seismic section obtained from three TFR’s (FSST, CWT, and STFT transforms) are shown in Figure 5. The common frequency slices (15 and 55 Hz) are illustrated in 3D and 2D views respectively in Figures 6 and 7. The positions of three lowfrequency shadows are indicated by three yellow rectangles in the 2D view. As can be seen, FSST has a much better resolution than the conventional spectral decompositions (STFT and CWT) and can potentially be used to detect hydrocarbons from a gas reservoir directly using lowfrequency shadows.
Figure 5
3D cubes of timefrequency transform of real 2D seismic data obtained from (a) STFT, (b) FSST, and (c) CWT methods; the positions of common frequency slices (30 and 60 Hz) on the 3D cubes are indicated by dashed lines.
Figure 6
3D view of common frequency slices (a) from STFT at 15 Hz, (b) from STFT at 55 Hz, (c) from CWT at 15 Hz, (d) from CWT at 55 Hz, (e) from FSST at 15 Hz, and (f) from FSST at 55 Hz.
5. Seismic random noise reduction
The sparse representation of signals is one of the useful properties of SSBT. In this paper, a new novel technique for random noise reduction is presented based on the mentioned property of the SSBT. Considering that seismic data are inherently lowrank (Ma, 2013) and this property is preserved after transforming seismic data to a new timefrequency domain, by using sparse TFR, the seismic data and their noise are represented as sparse as possible.
Consequently, it can be concluded that, if the proper sparse TFR is chosen, the decomposition of lowrank and sparse components can be useful for noise suppression. Nowadays, many methods are proposed for seismic random noise attenuation based on the extraction or estimation of the lowrank component from a noisy data (Cheng et al., 2015; Sacchi, 2009).
Figure 7
2D view of common frequency slices (a) from STFT at 15 Hz, (b) from STFT at 55 Hz, (c) from CWT at 15 Hz, (d) from CWT at 55 Hz, (e) from FSST at 15 Hz, and (f) from FSST at 55 Hz; the positions of three lowfrequency shadows are indicated by three yellow rectangles.
In these methods, classical principal component analysis (PCA) is the most widely used statistical tool for rank reduction. However, the validity of PCA has not been acceptable in grossly corrupted data. In fact, the principal components of data with a very low signaltonoise ratio will be changed. Therefore, many techniques are introduced to increase the robustness of PCA. Recently, robust principal component analysis (RPCA) (Candes et al., 2011) has been used to exactly decompose a signal into its lowrank and sparse components (Candes et al., 2011). In this method, nuclear norm ( ) and norm are utilized for separating the lowrank and sparse components of a signal respectively. The RPCA was herein utilized for decomposing the lowrank clean data from noisy observed data.
The proposed method first transformed the seismic trace into a new sparse subspace using SSWT. Then, the magnitude of sparse TFR matrix was decomposed into two parts, namely (a) a lowrank and (b) a sparse component, using RPCA. Now, let denote the magnitude of SSWT representation of the noisy observed signal by and represent its lowrank and sparse components by and respectively. Candes et al. (2011) showed that RPCA can solve the below constraint equation with unique answers:
(7) 
where, and are the maximum rank of lowrank component and the maximum number of nonzero elements in the sparse component respectively. One can recover the denoised seismic signal by minimizing the mixed norms objective function by considering the intrinsically lowrank property of the seismic data and the sparsity feature of the random noise in the sparse TFR domain. The mixed norms objective function is stated as follows:
(8) 
where, is the tradeoff parameter for balancing the sparsity condition and lowrank constraint. The recovered lowrank component is the SSWT magnitude of denoised data. Then, the clean data can be recovered by transforming back the denoised magnitude spectrum into the time domain.
5.1. Synthetic data
Herein, the 2D synthetic tx data set which has 30 traces with a total time of 0.36 s and a sampling interval of 2 ms is presented. The synthetic data includes two horizontal events, a slop event, and a curved event. The data are composed of Ricker wavelet with a central frequency of 30 Hz and contaminated with a random noise with a SNR equal to 1 dB. The original clean data and its noisy version are shown in Figure 8.
Both the proposed and classic singular spectrum analysis (SSA) methods (Sacchi, 2009) were applied to the noisy synthetic section. Figures 9a and 9b depict the results of denoising by the two mentioned methods. The parameter is set to 0.5 in the proposed method. Because of the existence of nonlinear event in the data, the rank parameter for SSA technique is selected to be 5. In the SSA method, there is a tradeoff between the noise reduction and event reconstruction. Therefore, if the higher values are chosen for rank in the SSA technique, the noise reduction will be decreased and the reconstruction of the event will be proper, and vice versa. In Figures 9c and 9d, the difference between the original noisy signal and the denoised version of the data are presented for comparing the performance of the two methods. It is clear that the efficiency of the method proposed in this work is much higher than the SSA algorithm. The standard SNR values calculated for the SSA and the proposed method are 0.8dB and +6.2dB respectively.
Figure 8
A synthetic seismic section; (a) clean data and (b) noisecontaminated data.
Figure 9
Denoising results for the synthetic data by (a) the proposed method and (b) the classical SSA; difference sections between the noisy and filtered data for (c) the proposed method and (d) the classical SSA.
5.2. Real data
Figure 10a shows a 44fold real CDP gather with 1015 time samples per trace with a time sampling interval of 0.003 s. Both methods (the proposed and the classical SSA) were applied to the real CDP gather. The rank parameter for the SSA method is set 15 and the parameter is chosen to be 5 for the proposed method. The obtained results by the two mentioned methods are presented in Figure 10b and 10c. For evaluating the efficiency of the two methods, the difference between the original input and the filtered gather was calculated (Figure 10d and 10e). As indicated, the performance of the proposed method in denoising the data, in contrast to the classical SSA method, is acceptable.
Figure 10
(a) A real CDP gather; denoising results for real CDP gather by (b) the proposed method and (c) the classical SSA; the difference sections between the noisy and filtered data for (d) the proposed method and (e) the classical SSA.
6. Conclusions
SSBT can be used to accurately map tx seismic signals into their TFR by relying on the frequency reassignment of CWT and STFT decompositions. It has a wellgrounded mathematical theory which facilitates the interpretations. Similar to other transform methods, it is a reversible function, which therefore allows for signal reconstruction, possibly after the removal of specific components such as noise. It can be seen that the SSBT provides a sparser image, which displays higher timefrequency resolution than the conventional methods in contrast to STFT and S transform. It was shown that FSST is a useful TFR to detect lowfrequency shadow anomalies with better resolution compared to STFT.
Then, a lowrank estimationbased scheme (RPCA) was introduced for seismic data denoising and the efficiency of the described method on synthetic and real data was shown. Because of the intrinsic low rank property of seismic data, this method can be applied to very different kinds of waves. The method was herein applied to synthetic data, to which random noise with an SNR of 1 dB was added. Then, both methods were applied to a real prestack CDP gather. The proposed method can improve signaltonoise ratio of the seismic data appropriately and provide higher signal fidelity compared with the classical SSA method.
Nomenclature
a 
: Scale 
A_{k}(t) 
: Instantaneous amplitude 
f_{k}(t) 
: Instantaneous frequency 
L 
: Lowrank component of matrix X 
S 
: Sparse component of matrix X 
t 
: Time 
T_{x}(w_{k},t) 
: Synchrosqueezing transform 
w_{x}(a,t) 
: Instantaneous frequency in wavelet domain 
W_{x}(a,t) 
: Coefficient of continuous wavelet transform of signal x(t) 
x(t) 
: Seismic trace in time domain 
x_{k}(t) 
: Intrinsic mode 
X 
: Magnitude of synchrosqueezing transform 
: Mother wavelet 

: Instantaneous phase 

: Tradeoff parameter 