Application of Local Wave Decomposition in Seismic Signal Processing

Abstract

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.

1. 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 time‐frequency methods including short‐time Fourier transform, wavelet transform, Wigner‐Ville distribution, and so on can appropriately describe the time variation of the nonstationary signal to some extent and show advantages over the Fourier transform. But overall, they are still in the scope of global wave.

Local wave decomposition (LWD) is a class of local wave processing methods that suit for the analysis of the nonlinear and nonstationary signals. LWD can decompose a complex multifrequency component signal into a number of monofrequency or narrow‐band signals, that is, intrinsic mode function (IMF), which has the better Hilbert transform characteristics and the calculation of its instantaneous frequency makes sense. LWD has evolved from empirical mode decomposition (EMD), which is first introduced by Huang et al. [1], to the recent works including ensemble empirical mode decomposition (EEMD) [2], complete ensemble empirical mode decomposition (CEEMD) [3], the synchrosqueezing transform (SST) [4], variational mode decomposition [5], and so on. Among these LWD methods, here we mainly study three typical methods: EMD, SST, and VMD. EMD is the method, which has already been widely used in time‐frequency applications in various signal processing such as climate analysis [6], seismic signal processing [7], mechanical fault diagnosis [8], speech signal processing [9], medicine and biology signal processing [10, 11], and so on. Due to the lack of mathematical understanding and some other obvious shortcomings including end effects [12], the problem of stopping criterion [13, 14], the influence of sampling [15], spline problems and mode mixing [16, 17], and so on in EMD algorithm, the other tools and improved methods pursuing the same goal are developed. SST is a wavelet‐based EMD‐like tool which has a firm theoretical foundation [18]. It introduces precise mathematical definition of a class of functions constructed by several limited approximate harmonic components. SST is a combination of wavelet analysis and frequency reassignment. While VMD is another newly adaptive and fully intrinsic method that has a firm theoretical foundation. It decomposes a nonlinear and nonstationary signal with quasi‐orthogonal, nonrecursive manner into a series of band‐limited subsignals.

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.

2. 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.

2.1. 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 zero‐crossings 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 $m(t)$ of the original signal $x(t)$ for the upper envelope $u(t)$ and the lower envelope $v(t)$, which are respectively fitted by all the maxima and minima of $x(t)$:

$m(t)=[u(t)+v(t)]/2$E1

in which, $v(t)\le x(t)\le u(t)$.

2. Subtract $m(t)$ from the original signal $x(t)$ and let $h_1(t)=x(t)-m(t)$.

3. Return to step (a) and replace $x(t)$ with $h_1(t)$. For $h_1(t)$, the corresponded upper and lower envelopes are $u_1(t)$ and $v_1(t)$, respectively. Repeat the above process until the resultant $h_k(t)$ meets the IMF definition:

$m_1(t)=[u_1(t)+v_1(t)]/2$E101

$h_2(t)=h_1(t)-m_1(t)$E102

$m_{k-1}(t)=[u_{k-1}(t)+v_{k-1}(t)]/2$E103

$h_k(t)=h_{k-1}(t)-m_{k-1}(t)$E104

Then the first IMF $c_1(t)=h_k(t)$. The residual part of the signal is: $r_1(t)=x(t)-c_1(t)$.

4. When there are more than one extremes (neither the constant nor the trend term) in $x(t)$, 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 $x(t)$ can be expressed as the sum of the IMFs and the margin:

$x(t)=c_1(t)+c_2(t)+\cdots +c_n(t)+r_n(t)={\displaystyle \sum _{i=1}^n\mathrm{Re}[a{}_i}(t)\exp (j{\theta }_i(t))]+r_n(t)$E2

where $c_i(t)$ is the ith IMF($i=1\sim n$ )and $r_n(t)$ is the residue. $a_i(t)$ and ${\theta }_i(t)$ are respectively the instantaneous amplitude and phase of the ith IMF $c_i(t)$.

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].

2.2. 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 $s(t)$ is [22]:

$W_s(a,b)=\frac{1}{\sqrt{a}}{\displaystyle \int s(t){\psi }^{*}}(\frac{t-b}{a})dt$E3

where a is the scale factor, $b$ is the time shift factor, and ${\psi }^{*}$ is the complex conjugate of the mother wavelet. The coefficients $W_s$(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:

$W_s(a,b)=\frac{1}{2\pi }{\displaystyle \int \frac{1}{\sqrt{a}}\widehat{s}(\xi ){\widehat{\psi }}^{*}}(a\xi )e^{jb\xi }d\xi$E4

where $\xi$ is the angular frequency. $\widehat{\psi }(\xi )$ and $\widehat{s}(\xi )$ are respectively the Fourier transform of $\psi (t)$ and $s(t)$. $j=\sqrt{-1}$.

Considering a single harmonic signal with the form:

$s(t)=A\cos (\omega t)$E5

Its Fourier pair can be expressed as

$\widehat{s}(\xi )=\pi A[\delta (\xi -\omega )+\delta (\xi +\omega )]$E6

Then Eq. (4) can be transformed into

$\begin{array}{l}{W}_s(a,b)=\frac{A}{2}{\displaystyle \int \frac{1}{\sqrt{a}}[\delta (\xi -\omega )+\delta (\xi +\omega )]{\widehat{\psi }}^{*}}(a\xi )e^{jb\xi }d\xi \\ =\frac{A}{2\sqrt{a}}{\widehat{\psi }}^{*}(a\omega )e^{jb\omega }\end{array}$E7

Since the wavelet ${\widehat{\psi }}^{*}(\xi )$ is condensed around its central frequency ${\omega }_0$, $W_s$(a,b) will be condensed around the horizontal line $a={\omega }_0/\omega$, which represents the ratio of the central frequency of the wavelet to the central frequency of the signal. But in fact $W_s$(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 $b$. And little smearing in dimension $b$ along the scale axis can be found. It is proved that the instantaneous frequency ${\omega }_s$(a,b) can be computed using the following equation when the smear in the b‐dimension can be neglected:

${\omega }_s(a,b)=\frac{-j}{W_s(a,b)}\frac{\partial W_s(a,b)}{\partial b}$E8

where for any point (a,b), $W_s\mathrm{}$(a,b) ≠ 0.

Then map the information from the time‐scale plane to the time‐frequency plane and convert every point $(b,a)$ to $(b,{\omega }_s(a,b))$, which is named as synchrosqueezing operation. Since $a$ and $b$ are discrete values, a scaling step $\Delta a_k=a_{k-1}-a_k$ can be used for any $a_k$, where $W_s(a,b)$ is computed. Then the SST $T_s(w,b)$ is determined only at the centers ${\omega }_l$ with the frequency range $[{\omega }_l-\Delta \omega /2,{\omega }_l+\Delta \omega /2]$:

$T_s({\omega }_l,b)=\frac{1}{\Delta \omega }{\displaystyle \sum _{a_k:\left|\omega (a_k,b)-{\omega }_l\right|\le \Delta \omega /2}{W}_s(a_k,b)a^{-3/2}}\Delta a_k$E9

where $\Delta \omega ={\omega }_l-{\omega }_{l-1}$.

Eq. (9) shows that the new time‐frequency representation of the signal $T_s(w,b)$ 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 $s_k$ can be reconstructed by the discrete SST ${\tilde{T}}_{\tilde{s}}$:

$s_k(t_m)=2C_{\phi }^{-1}\mathrm{Re}({\displaystyle \sum _{l\in L_k(t_m)}{\tilde{T}}_{\tilde{s}}}(w_l,t_m))$E10

where $C_{\phi }^{}$ is a constant dependent on the selected wavelet. $\mathrm{Re}(•)$ takes the real part of the representation. ${\tilde{T}}_{\tilde{s}}$ is the discretized version of $T_s({\omega }_l,b)$, that is, ${\tilde{T}}_{\tilde{s}}(w_l,t_m)$. $t_m$ is the discrete time, $t_m=t_0+m\Delta t$, $m=0,1,\dots ,n-1$. $n$ is the total number of samples in the discrete signal ${\tilde{s}}_m$. $\Delta t$ is the sampling rate.

2.3. 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]:

$u(t)=A(t)\cos \left(\phi (t)\right)$E11

where $u(t)$ resembles an IMF which is a narrow‐band amplitude‐modulated, frequency‐modulated signal. The phase $\phi (t)$ is a nondecreasing function. The instantaneous frequency $\omega (t)={\phi }^{\prime }(t)$ and the envelope $A(t)$ are nonnegative. $\omega (t)$ and $A(t)$ vary much slower than $\phi (t)$ itself.

VMD is a constrained variational problem as shown by Eq. (12) [5, 25]:

$\underset{\left\{u_k\right\},\left\{{\omega }_k\right\}}{\min }\left\{{\displaystyle \sum _k^{}{\Vert {\partial }_t\left[(\delta (t)+\frac{j}{\pi t})\ast u_k(t)\right]e^{-j{\omega }_{k}t}\Vert }_2^2}\right\}$E12

subject to ${\displaystyle \sum _k^{}{u}_k}=f$,

where $u_k$ and ${\omega }_k$ are respectively the kth mode and their center frequency. For each mode, its bandwidth is estimated through the H¹ 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:

$\begin{array}{l}L\left(u_k,{\omega }_k,\lambda \right)=\alpha {{\displaystyle \sum _k\Vert {\partial }_t\left[\left(\delta (t)+\frac{j}{\pi t}\right)\ast u_k(t)\right]e^{-j{\omega }_{k}t}\Vert }}_2^2\\ +{\Vert f-{\displaystyle \sum u_k}\Vert }_2^2+⟨\lambda ,f-{\displaystyle \sum u_k}⟩\end{array}$E13

where $\alpha$ represents the balancing parameter of the data‐fidelity constraint.

The main steps of VMD include the following:

1. Modes $u_{{}_k}^{}$ update. Eq. (14) shows how the modes ${\widehat{u}}_{{}_k}^{n+1}\left(\omega \right)$ are updated:

   ${\widehat{u}}_{{}_k}^{n+1}\left(\omega \right)=\frac{\widehat{f}\left(\omega \right)-{\displaystyle {\sum }_{i<k}{\widehat{u}}_{{}_i}^{n+1}\left(\omega \right)-{\displaystyle {\sum }_{i>k}{\widehat{u}}_{{}_i}^n\left(\omega \right)+}\left({\widehat{\lambda }}^n\left(\omega \right)/2\right)}}{1+2\alpha {\left(\omega -{\omega }_k^n\right)}^2}$E14

   Eq. (14) clearly demonstrates that the modes ${\widehat{u}}_{{}_k}^{n+1}\left(\omega \right)$ acts as a Wiener filtering of the current residual with signal prior $1/{(\omega -{\omega }_k^n)}^2$. The mode in time domain can be obtained by taking the real part of the inverse Fourier transform of this filtered analytic signal.

2. Center frequencies ${\omega }_k^{n+1}$ update. The new center frequencies ${\omega }_k^{n+1}$ are put at the center of gravity of the corresponding mode's power spectrum which are shown in Eq. (15):

   ${\omega }_k^{n+1}=\frac{{\displaystyle {\int }_0^{\infty }\omega {\left|{\widehat{u}}_k^{n+1}\left(\omega \right)\right|}^{2}d\omega }}{{\displaystyle {\int }_0^{\infty }{\left|{\widehat{u}}_k^{n+1}\left(\omega \right)\right|}^{2}d\omega }}$E15

3. Dual ascent update. The Lagrangian multiplier ${\widehat{\lambda }}^{n+1}$ is used as a dual ascent to enforce exact signal reconstruction. For all $\omega \ge 0$, ${\widehat{\lambda }}^{n+1}$ is updated as shown in Eq. (16) until convergence ${{\displaystyle {\sum }_k\Vert {\widehat{u}}_k^{n+1}-{\widehat{u}}_k^n\Vert }}_2^2/{\Vert {\widehat{u}}_k^n\Vert }_2^2<\epsilon$

   ${\widehat{\lambda }}^{n+1}={\widehat{\lambda }}^n+\tau \left(\widehat{f}-{\displaystyle \sum _k{\widehat{u}}_k^{n+1}}\right)$E16

The detailed complete algorithm of VMD can be found in [5].

3. Applications

3.1. 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 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.

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.

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.

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.

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.

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 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.

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.

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).

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.

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.

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) EMD‐based method; (e) SST‐based method; (f) VMD‐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.55 s (Figure 7a) to the synthetic signal in Figure 1.

Figure 7.

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

When EMD is applied to this noisy signal, seven IMFs and one residue are obtained (Figure 8). As Figure 8 shows, the first IMF shows mode mixing phenomenon and it does not represent the either component of the signal. Influenced by the first IMF, the following IMFs are all disturbed and lose their physical meaning.

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 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.

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 shows 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 (Figure 11a). As shown in Figure 10(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 (Figures 11b and c), the bandwidth of IMF1 to IMF2 is similar to that in Figures 5(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.

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.

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.

3.2. 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.

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 Figure 13. Since the random noise always distributes in the high‐frequency 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 (Figure 14). Compared with the original seismic section in Figure 12, the details in the de‐noising seismic section are more clear.

3.3. 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).

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.

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 15.

The broadband migrated, stacked seismic section.

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 Figure 17. Compare the four IMFs in Figure 17 with the original seismic section in Figure 15, we can find that the second IMF (Figure 17b) targets the main reflection information and reflects more details. Then, we use the second IMF for edge detection.

Figure 17.

The four IMFs after 2D‐VMD for the broadband migrated stacked seismic section. (a) IMF1; (b) IMF2; (c) IMF3; and (d) IMF4.

Next, Canny edge detection is carried out to the second IMF. The output of the Canny edge detection is shown in Figure 18. Compared with the results in Figure 16, the output of the Canny edge detection applied to the second IMF shows more continuity and the changes in detail and also can be helpful for the automatic seismic event pickup.

Figure 18.

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

3.4. Application of LWD on the recovery of the target reflection near coal seams

3.4.1. 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 19.

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

	①	②	③	④	⑤
V (m/s)	3960	4393	4375	4557	3608
Den (kg/m^3)	2520	2549	2495	2434	2296

Table 1.

The parameters of each layer for the model.

Note: V represents the compressional velocity. Den is the density.

VMD is a data‐driven adaptive band‐limited filtering method. After VMD process, the original seismic data will be decomposed into a series of IMFs. For the coal‐bearing strata coexisted situation, seam stratigraphic features are obvious on one IMF section and weak reflection formation signals are not well represented. Here, the main IMF contributors to coal seams are selected by using the combination of logging information and the maximum correlation. For each selected IMFs, we calculate the energy trace by trace to find out the maximum energy point $t_{\max }$. Then, estimate the main frequency $f_d$ of the IMF of the seismic trace. Time thickness of the thin layer of coal seams is determined by $[t_{\max }-k_1{T}_d,t_{\max }-k_2{T}_d]$, where $k_1,k_2$ are constant factors, respectively, and are estimated by the well‐seismic calibration, $T_d=1/f_d$.

Finally, let $index=\frac{E_{ave}}{E_s}$. Using $index$ to suppress the strong amplitude to be consistent with the average energy of the corresponding seismic trace, in which $E_s$ is the energy of the seismic data within the range of $[t_{\max }-k_1{T}_d,t_{\max }-k_2{T}_d]$, $E_{ave}$ represents the average energy of the seismic trace.

When all the IMFs that contributed to the main coal seams are processed, sum up the processed IMFs and the remaining IMFs unchanging to reconstruct the seismic trace. Repeat the above process trace by trace, and the strong amplitude suppression section is obtained. This method only suppresses the strong amplitudes of the main coal seams, IMF contributors and the other IMFs, which reflect the remaining information. So the VMD‐based method enhances the weak signal components of the coexistence coal‐bearing strata.

The result of using the VMD method is shown in Figure 20. As shown in Figure 20, the strong amplitudes in coal seams are suppressed and the weak reflections near coal seams are enhanced.

3.4.2. Real data

Next, the broadband migrated stacked seismic data from the Ordos Basin, China is collected for analysis, as shown in Figure 21. The data is sampled at 2 ms. Strong reflection amplitudes exist in the area where the coal seam is located at around 1900 ms. The sandstone information is with weak amplitudes. The details near coal seams are very weak and evenly hidden by the influence of the strong amplitudes in coal seams.

Figure 20.

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

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.

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

Figure 23 shows the instantaneous amplitudes comparison between the original seismic section and the seismic section after VMD‐based strong amplitude suppression method. From Figure 23(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.

The instantaneous amplitudes comparison between the original seismic section (a) and the seismic section after VMD‐based strong amplitude suppression method (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.

Figure 24.

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

4. 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 Wavelet‐based 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.

Acknowledgments

This work was supported in part by National Natural Science Foundation of China (Grant Nos. 41404102, 41430323, and 41274128) and Sichuan Youth Science and Technology Foundation (Grant No. 2016JQ0012) and in part by Key Project of Sichuan provincial Education Department (No.16ZA0218) and the 2015 annual young Academic Leaders Scientific Research Foundation of CUIT (No. J201507) and the Project of the Scientific Research Foundation of CUIT (No. KYTZ201503). The authors also thank for the supportion of SINOPEC Key Laboratory of Geophysics.