Application of Local Wave Decomposition in Seismic Signal Processing

Local wave decomposition (LWD) method plays an important role in seismic signal processing for its superiority in significantly revealing the frequency content of a seismic signal changes with time variation. The LWD method is an effective way to decompose a seismic signal into several individual components. Each component represents a harmonic signal localized in time, with slowly varying amplitudes and frequencies, potentially highlighting different geologic and stratigraphic information. Empirical mode decomposition (EMD), the synchrosqueezing transform (SST), and variational mode decomposition (VMD) are three typical LWD methods. We mainly study the application of the LWD method especially EMD, SST, and VMD in seismic signal processing including seismic signal de‐noising, edge detection of seismic images, and recovery of the target reflection near coal seams.


Introduction
The main characteristics of a nonlinear and nonstationary signal are that the signal is changed with the time and its frequency is transient and only in the presence of a local time, that is, a local wave signal.The conventional Fourier transform is only suitable for the global wave, and it shows great limitation for the analysis of the local wave.The thereafter appeared timefrequency methods including short-time Fourier transform, wavelet transform, Wigner-Ville Since a seismic signal is a nonlinear and nonstationary signal, the LWD method is more suitable for the analysis of the seismic signal than the other traditional Fourier-and wavelet-based methods.Each obtained IMF with different narrow frequency band from LWD method represents a harmonic signal localized in time, with slowly varying amplitudes and frequencies, potentially highlighting different geologic and stratigraphic information and their pore fluids.
In this chapter, we mainly study the application of the LWD method especially three typical LWD methods, that is, EMD, SST, and VMD, in seismic signal processing including seismic signal de-noising, edge detection of seismic images, and recovery of the target reflection near coal seams.

The general theory of local wave decomposition
In this chapter, not all the LWD methods that might be useful in practice can be examined.Guided by the various previous literatures, we only briefly review a few concepts and tools on the most commonly used and three new typical LWD methods: EMD, SST, and VMD.

Empirical mode decomposition
The purpose of EMD is to obtain IMFs.Huang et al. [1] believe that any data consists of different simple intrinsic modes of oscillations, that is, IMFs.Each IMF involves only one mode of oscillation due to the forbiddance of no complex riding waves.With respect to the "local mean," the oscillation will also be symmetric.A meaningful instantaneous frequency of a multicomponent signal can be obtained by reducing the original signal to a collection of IMFs by employing EMD through a sifting process [19].Only the instantaneous frequency of an IMF has physical meaning [1,19].An IMF is defined as a signal whose number of extrema and zerocrossings must either equal or differ at most by one and the mean value of the upper and lower envelopes, respectively, defined by the local maxima and minima is zero [1].The IMF definition guarantees the narrow-band requirement for a stationary Gaussian process and the unwanted fluctuations induced by a symmetric waveform will not be emerged in the instantaneous frequency [1,20].The IMF is a linear or nonlinear signal that has constant amplitude and frequency as in a simple harmonic component or variable amplitude and frequency as functions of time.
EMD mainly includes the following steps: 1. Compute the average curve () of the original signal () for the upper envelope () and the lower envelope (), which are respectively fitted by all the maxima and minima of (): in which, () ≤ () ≤ ().
4. When there are more than one extremes (neither the constant nor the trend term) in (), the EMD process is conducted until the remaining portion of the resulting signal is a monotone or a value less than a predetermined given value.
After the EMD process, the original signal () can be expressed as the sum of the IMFs and the margin: where () is the ith IMF( = 1 ∼ )and () is the residue. ()and () are respectively the instantaneous amplitude and phase of the ith IMF ().
As shown above, the local mode of the signal is successively isolated from high frequency to low frequency based on the characteristic time scales by the sifting process in EMD.It is proved that EMD acts as an adaptive, multiband overlapping filter bank [21].

The synchrosqueezing transform
SST algorithm is another technique that aims to decompose a multicomponent signal into a series of IMFs like EMD.But unlike EMD, it introduces a precise mathematical definition for IMFs that can be viewed as a superposition of a reasonably small number of approximate harmonic components [4].SST method uses frequency reassignment to improve the readability of wavelet-based time-frequency analysis.
The continuous wavelet transform (CWT) of a signal () is [22]: where a is the scale factor, is the time shift factor, and * is the complex conjugate of the mother wavelet.The coefficients (a,b) which represent a concentrated time-frequency picture can be used to extract the instantaneous frequencies [22].Rewrite Eq. (3) using Plancherel's theorem which states that the energy in the time domain equals energy in the frequency domain: where is the angular frequency. ()and () are respectively the Fourier transform of () and ().= −1.
Considering a single harmonic signal with the form: Its Fourier pair can be expressed as Then Eq. ( 4) can be transformed into Since the wavelet *() is condensed around its central frequency 0 , (a,b) will be condensed around the horizontal line = 0 /, which represents the ratio of the central frequency of the wavelet to the central frequency of the signal.But in fact (a,b) always spreads out around this horizontal line, which makes a blurred projection in time-scale representation.This main smearing occurs in the scale dimension along the time shift factor .And little smearing in dimension along the scale axis can be found.It is proved that the instantaneous frequency (a,b) can be computed using the following equation when the smear in the b-dimension can be neglected: where for any point (a,b), (a,b) ≠ 0.
Then map the information from the time-scale plane to the time-frequency plane and convert every point (, ) to (, (, )), which is named as synchrosqueezing operation.Since and are discrete values, a scaling step = − 1 − can be used for any , where (, ) is computed.Then the SST (, ) is determined only at the centers with the frequency range Eq. (9) shows that the new time-frequency representation of the signal (, ) is synchrosqueezed along the frequency axis only [23].The coefficients of the CWT are reallocated by SST to get a concentrated image over the time-frequency plane.The instantaneous frequency is also extracted from the new time-frequency representation [24].
The individual components can be reconstructed by the discrete SST : where is a constant dependent on the selected wavelet.Re( • ) takes the real part of the representation.is the discretized version of ( , ), that is, ( , ). is the discrete time, = 0 + , = 0, 1, …, − 1. is the total number of samples in the discrete signal .
is the sampling rate.

Variational mode decomposition
VMD is a newly developed and theoretically well-founded EMD-like tool.It can nonrecursively decompose a multicomponent signal into a series of band-limited IMFs.Note here, IMF is defined more restrictively in VMD compared with the original IMF definition in EMD as the following [5]: where () resembles an IMF which is a narrow-band amplitude-modulated, frequencymodulated signal.The phase () is a nondecreasing function.The instantaneous frequency () = ′() and the envelope () are nonnegative.() and () vary much slower than () itself.
VMD is a constrained variational problem as shown by Eq. ( 12) [5,25]: Earthquakes -Tectonics, Hazard and Risk Mitigation where and are respectively the kth mode and their center frequency.For each mode, its bandwidth is estimated through the H 1 Gaussian smoothness of the demodulated signal.To address Eq. ( 12), a quadratic penalty and Lagrangian multipliers are used.The augmented Lagrangian is introduced as follows: ( ) , where represents the balancing parameter of the data-fidelity constraint.
The main steps of VMD include the following: 1. Modes update.Eq. ( 14) shows how the modes + 1 are updated: Eq. ( 14) clearly demonstrates that the modes + 1 acts as a Wiener filtering of the current residual with signal prior 1/( − ) 2 . The mode in time domain can be obtained by taking the real part of the inverse Fourier transform of this filtered analytic signal.

Center frequencies
+ 1 update.The new center frequencies + 1 are put at the center of gravity of the corresponding mode's power spectrum which are shown in Eq. ( 15): 3. Dual ascent update.The Lagrangian multiplier + 1 is used as a dual ascent to enforce exact signal reconstruction.For all ≥ 0, + 1 is updated as shown in Eq. ( 16) until convergence ∑ The detailed complete algorithm of VMD can be found in [5].

Synthetic data
In this section, synthetic data is first used to compare EMD with the SST and VMD methods.
Then, the equivalent band-limited filter behavior of EMD, SST, and VMD are respectively investigated.Finally, to evaluate the noise robustness of the three methods, the synthetic data with added noise is analyzed.
The synthetic signal is shown in  When EMD is applied to the synthetic signal, three IMFs and one residue are obtained (Figure 2).It is clear that IMF1 does not represent the background cosine wave or the Ricker wavelets singly.Mode mixing occurred seriously in the IMF1.Therefore, the IMF1 lacks physical meaning.Furthermore, affected by the IMF1, the following IMFs are all distorted.The SST result is shown in Figure 3.We can find that the reconstructed components of the SST better represent the background cosine wave and the Ricker wavelets, respectively.Mode mixing is better suppressed than EMD.The spectrums of the EMD output, the SST output, and the VDM output are shown in Figure 5.We can find that an overlap of half bandwidth occurred between the two adjacent IMFs for the IMFs in EMD (Figure 5a).But for the spectrum of the SST output and the VMD output, the band-pass filters characteristics are shown with the increasing predominant center frequencies (Figure 5b and c).The time-frequency spectrum comparison of the EMD-, SST-, and VMD-based method with short-time Fourier transform (STFT) and wavelet transform is shown in Figure 6.Here, the EMD-, SST-, and VMD-based method are respectively using EMD, SST, and VMD combined with Hilbert transform to provide the time-frequency distribution.The EMD-, SST-, and VMD-based method show higher temporal and spatial resolution than STFT and wavelet transform.
We also can find that SST-and VMD-based method give more precise time-frequency distribution of the component than EMD-based method.To evaluate the noise robustness of the three methods, we use a noisy signal (Figure 7b) by adding the Gaussian noise only distributed within the time 0.45-0.55s (Figure 7a) to the synthetic signal in Figure 1.

Seismic data
Here, one 2D prestack seismic section with random noise interference from Ordos Basin in China is used for analysis (Figure 12).The data is sampled at 1 ms.The SST method is mainly used for seismic data de-noising.Then, the de-noising component and the first component are used to generate the de-noising seismic section (Figure 14).Compared with the original seismic section in Figure 12, the details in the de-noising seismic section are more clear.

Edge detection of seismic images using LWD
The study on the edge in seismic images has very important significance.Crack, fracture, fault, small rupture, river channel sand boundary, and the boundary of the other special lithology body, and so on all show the edge characteristics in seismic images.In this section, we mainly applied VMD-based method for edge detection in seismic images.The details of 2D-VMD algorithm can be found in reference [26].
Here, the broadband migrated stacked seismic data from a gas field located in western Sichuan Depression, China is collected for analysis (Figure 15).If we apply Canny edge detection method to the original seismic section, the output shows more discontinuity and mess (Figure 16).

Model analysis
To evaluate the LWD method for suppressing the strong seismic reflection amplitude of coal seams and recovery of the target reflection near coal seams, one model is designed to simulate the seismic response based on the seismic data and reservoir logging parameters of one tight sandstone reservoir from the Ordos Basin, China.
There are five formations in the geological model 1 (Figure 19a).The parameters of each layer are shown in Table 1.The layer marked ③ is gas-bearing layer.The layer marked ④ is coal seams.Sampling frequency is 500 Hz.The frequency of the wavelet is 50 Hz.The corresponding seismic response of the geological model is shown in Figure 19(b).Figure 24 shows the instantaneous frequencies comparison between the original seismic section and the seismic section after VMD-based strong amplitude suppression method.We can find that more frequency details are enhanced in the coal seams area around 1900 ms.

Conclusion
Analysis on synthetic and real data shows that the LWD method is more robust to noise and has stronger local decomposition ability than the traditional Fourier-based and Waveletbased methods.Comparing with the short-time Fourier transform or wavelet transform, LWD-based time-frequency spectrum promises higher temporal and spatial resolution.Application of the LWD on field data demonstrates that the LWD method can effectively be used in seismic signal de-noising and edge detection.Also, LMD-based method targets the thickness variation of coal seams sensitively and suppresses the strong amplitude in coal seams effectively and highlights the fine details that might escape unnoticed near the area where coal seams are located.The LWD-based technique is more promising for seismic signal processing and interpretation.

Figure 1 .
It is comprised of an initial 20 Hz cosine wave in which 90 Hz Ricker wavelets at 0.2 s, 75 Hz Ricker wavelets at 0.4 s, and two 40 Hz Ricker wavelets at 0.68 and 0.72 s are respectively superposed.

Figure 1 .
Figure 1.The synthetic signal.It is comprised of an initial 20 Hz cosine wave in which 90 Hz Ricker wavelets at 0.2 s, 75 Hz Ricker wavelets at 0.4 s, and two 40 Hz Ricker wavelets at 0.68 and 0.72 s are, respectively, superposed.

Figure 2 .
Figure 2. The EMD result of the synthetic signal.Mode mixing is clearly seen in the IMF1.It makes the IMF1 lack the physical meaning.Furthermore, the following IMFs are all distorted.

Figure 3 .
Figure 3.The SST result of the synthetic signal.The reconstructed components of the SST better represent the background cosine wave and the Ricker wavelets, respectively.Mode mixing is better suppressed in SST than in EMD.

Figure 4
Figure 4 shows the VMD result of the synthetic signal.As shown in Figure 4, the background cosine wave is first extracted in IMF1.Then the 90 Hz Ricker wavelets at 0.2 s, 75 Hz Ricker wavelets at 0.4 s, and two 40 Hz Ricker wavelets at 0.68 and 0.72 s are mainly retrieved in IMF2.The IMFs of VMD can better represent the individual component of the original signal than the EMD.It shows that the IMFs have more physical meaning.

Figure 4 .
Figure 4.The VMD result of the synthetic signal.The background cosine wave is first extracted in IMF1.Then the Ricker wavelets are mainly retrieved in IMF2.Mode mixing is suppressed the most.The IMFs have more physical meaning.

Figure 5 .
Figure 5.The spectrums of the EMD output (a), the SST output (b), and the VDM output (c).An overlap of half bandwidth is occurred between the two adjacent IMFs for the IMFs in EMD.

Figure 6 .
Figure 6.The time-frequency spectrum comparison of the EMD-, SST-, and VMD-based method with short-time Fourier transform (STFT) and wavelet transform for the synthetic data.(a) The synthetic data; (b) STFT; a Hamming window with the length 61 is used.(c) Wavelet transform; continuous wavelet transform with Morlet wavelet is used.(d) EMDbased method; (e) SST-based method; (f) VMD-based method.

Figure 7 .
Figure 7.A noisy signal (b) by adding the Gaussian noise only distributed within the time 0.45-0.55s (a) to the synthetic signal in Figure 1.

Figure 8 .
Figure 8. EMD output of the noisy signal.Mode mixing occurred in the first IMF.And influenced by the first IMF, the following IMFs are all disturbed.

Figures 9 8 ,
Figures 9 and 10 respectively show the SST output and the VMD output of the noisy signal.The background cosine wave is extracted in IMF1 in both the first reconstructed component of SST and the first IMF of the VMD output.The Ricker wavelets are mainly reflected in IMF2 with a slight influence of noise.And the noise is mainly reflected in IMF3.Compared with Figure 8, SST and VMD show the noise robustness more than the EMD method.

Figure 9 .
Figure 9. SST output of the noisy signal.The background cosine wave is retrieved in IMF1 in both the first reconstructed component of SST and the first IMF of the VMD output.The noise is mainly reflected in IMF3.The Ricker wavelets are mainly retrieved in IMF2 with a slight influence of noise.

Figure 11
Figure11shows the spectrum of the EMD output, the SST output, and the VMD output for the noisy signal.Note that only the spectrum of the first three IMFs are shown in EMD (Figure11a).As shown in Figure10(a), the bandwidths of IMFs in EMD output are overlapped.Due to the noise influence, the bandwidth of IMF1 becomes larger.But for the spectrums of the SST output and the VMD output (Figures11b and c), the bandwidth of IMF1 to IMF2 is similar to that in Figures5(b) and (c).The noise is separated and mainly reflected in IMF3 for the spectrum of the SST output and the VMD output.Compared Figures 11(b) and (c) with EMD output in Figure 10(a), SST and VMD show better noise robustness.

Figure 10 .
Figure 10.VMD output of the noisy signal.The background cosine wave is mainly extracted in IMF1 in both the first reconstructed component of SST and the first IMF of the VMD output.The noise is mainly reflected in IMF3.The Ricker wavelets are mainly retrieved in IMF2 with a slight influence of noise.

Figure 11 .
Figure 11.The spectrum of the EMD output (a), the SST output (b), and the VMD (c) output for the noisy signal.SST and VMD show better noise robustness than the EMD does.

Figure 12 .
Figure 12.The 2D prestack seismic section with random noise interference from Ordos Basin, China.By analysis, the main frequency of the seismic section ranges from 0 to 60 Hz.The dominant frequency is 33 Hz.For SST, we use the frequency range [0, 33] to reconstruct the first component and the frequency range [34, 60] for the second component.And the residue frequencies are used to reconstruct the third component.The reconstructed components of the SST for the seismic section are shown in Figure13.Since the random noise always distributes in the highfrequency component, we apply soft thresholding on wavelet transform to the second reconstructed components and abandon the third component which mainly reflects the noise.Then, the de-noising component and the first component are used to generate the de-noising seismic section (Figure14).Compared with the original seismic section in Figure12, the details in the de-noising seismic section are more clear.

Figure 13 .
Figure 13.The reconstructed components of the SST for the seismic section.The frequency range [0, 33] is used to reconstruct the first component (a) and the frequency range [34, 60] for the second component (b).The residue frequencies are used to reconstruct the third component (c).

Figure 14 .
Figure 14.The de-noising seismic section by using SST.Compared with the original seismic section in Figure 12, the details in the de-noising seismic section are more clear.

Figure 16 .
Figure 16.The output of the Canny edge detection for the seismic section.After 2D-VMD process, the four IMFs are obtained as shown in Figure17.Compare the four IMFs in Figure17with the original seismic section in Figure15, we can find that the second IMF (Figure17b) targets the main reflection information and reflects more details.Then, we use the second IMF for edge detection.

Figure 18 .
Figure 18.The output of the Canny edge detection for the second IMF.

Figure 19 .
Figure 19.The geological model (a) and the corresponding seismic response (b).

Figure 20 .
Figure 20.The result of using the VMD method to the model for suppressing the strong amplitudes of the coal seams.

Figure 21 . 22 .
Figure 21.The broadband migrated, stacked seismic data.The data is sampled at 2 ms.There exist strong reflection amplitudes where coal seams are located around 1900 ms.The sandstone information is with weak amplitudes.The result of using the VMD-based strong amplitude suppression method is shown in Figure 22.As shown in Figure 22, the strong amplitudes in coal seams are suppressed and the weak reflections near coal seams are enhanced.

Figure 22 .
Figure22.The result of using the VMD method to the seismic section for suppressing the strong amplitudes of the coal seams.

Figure 23
Figure23shows the instantaneous amplitudes comparison between the original seismic section and the seismic section after VMD-based strong amplitude suppression method.From Figure23(b), we can find that the strong amplitudes in coal seam are suppressed most and the details near the coal seams are highlighted.

Figure 23 .
Figure 23.The instantaneous amplitudes comparison between the original seismic section (a) and the seismic section after VMD-based strong amplitude suppression method (b).

Figure 24 .
Figure 24.The instantaneous frequencies comparison between the original seismic section (a) and the seismic section after VMD-based strong amplitude suppression method (b).