For a time-domain signal x t , the time–frequency spectrum S t , f based on the ST is expressed as Equation 1 (Stockwell et al., 1996) S τ , f = ∫ − ∞ + ∞ x t g t − τ , f e − i 2 π f t d t , (1) where g t is the Gaussian window of the specific form f 2 π e − t 2 f 2 2 , τ is a translation factor to control the position of the Gaussian window on the time axis t , f represents the frequency, and i is imaginary units.

Then, the depth–wavenumber spectrum of a depth-domain signal based on the ST can be expressed as Equation 2 S η , k = ∫ − ∞ + ∞ x h g h − η , k e − i 2 π k h d h , (2) where g h represents the depth-domain Gaussian window of the specific form k 2 π e − h 2 k 2 2 , η is a translation factor to control the position of the Gaussian window on the depth axis h , and k represents the wavenumber.

To reduce the bias of the frequency components in the time–frequency spectrum based on the ST, the time-domain UST removes the linear frequency-dependent term f of g t . Then, the depth-domain UST is derived as S 1 η , k = ∫ − ∞ + ∞ x h 1 2 π e − h − η 2 k 2 2 e − i 2 π k h d h . (3)

Then, the expression to the right of the equal sign in Equation 3 in the wavenumber domain (i.e., WUST) can be written as S 1 η , k = ∫ − ∞ + ∞ X α + k 1 k e − 2 π 2 α 2 k 2 e i 2 π α η d α , k ≠ 0 , (4) where X α + k is the Fourier transform of x h e − i 2 π k h . α denotes the translating wavenumber (Wang, 2016).

However, the UST sacrifices time/depth resolution at low frequencies/low wavenumbers in order to obtain a reliable frequency/wavenumber distribution. To overcome this shortcoming, slope ( A ) and intercept ( B ) parameters are introduced to obtain the desired depth–wavenumber resolution. The equation of the modified UST in the depth domain is S 2 η , k = ∫ − ∞ + ∞ x h 1 2 π e − h − η 2 A k + B 2 2 e − i 2 π k h d h . (5)

Based on the convolution theorem, the expression to the right of the equal sign in Equation 5 in the wavenumber domain (i.e., MWUST) can be written as S 2 η , k = ∫ − ∞ + ∞ X α + k 1 A k + B e − 2 π 2 α 2 A k + B 2 e i 2 π α η d α , k ≠ 0 , (6) where X α + k is the Fourier transform of x h e − i 2 π k h . The constant average of the signal x h is put into zero wavenumber, which ensures the feasibility of the inverse transform. By integrating and inverse Fourier-transforming along different axes, we can reconstruct the original signal. It is worth noting that Equation 6 degenerates to Equation 4 when A = 1 and B = 0 .

There are two parameters (i.e., A and B) in the proposed method that need to be determined in advance. For different tasks and frequencies, the values of A and B may vary, which requires manual adjustment. One can use trial-and-error methods to obtain relatively optimal results. However, before that, we can use the full-window spatial or wavenumber width at half-maximum (FWHM) to roughly estimate parameters A and B (Li et al., 2016; Zhang et al., 2022b). For the depth and wavenumber domains, there are different expressions for A and B.

For the depth domain, A and B have the following expressions (Equation 7) (George et al., 2009): A = 2.355 1 Δ h 2 − 1 Δ h 1 k 2 − k 1 B = 2.355 Δ h 1 − A k 1 , (7) where ∆ h 1 and ∆ h 2 represent specified spatial (or depth) FWHM resolutions at wavenumbers k 1 and k 2 , respectively.

For the wavenumber domain, A and B have the following expressions (Equation 8): A = 2.668 Δ α 2 − Δ α 1 k 2 − k 1 B = 2.668 Δ α 1 k 2 − Δ α 2 k 2 k 2 − k 1 , (8) where ∆ α 1 and ∆ α 2 represent specified wavenumber FWHM resolutions at wavenumbers k 1 and k 2 , respectively.

Figure 1 shows relationships between FWHM amplitudes and wavenumber for different A and B expressions. In Figure 1A, the width of the MWUST is narrower than that of the WUST in the low-wavenumber region when A > 1 (i.e., orange-, green-, and gray-solid curves). It indicates that the MWUST has higher resolution in the low-wavenumber region, and the difference decreases as the wavenumber increases. When A < 1 , shown as a purple solid curve, the width of the MWUST is broader than that of the WUST in the high-wavenumber region, indicated by larger depth FWHM values. The wavenumber FWHM values increase linearly as the wavenumber increases, as shown in Figure 1B. From Figure 1, it can also be noticed that A plays a greater role for FWHM values than B, indicating that A has a great influence on the depth–wavenumber resolution. B controls the starting depth resolution.

The depth-domain non-stationary convolution model is expressed as u x = ∫ − ∞ + ∞ w x − η , η r η d η , (9) where u x represents depth-domain synthetic seismic records, w x , η represents the non-stationary depth-domain seismic wavelet, and r η represents the depth-domain reflectivity series. Equation 9 is transformed into the wavenumber domain with the following expression (Equation 10): U k = ∫ − ∞ + ∞ R ξ W k , k − ξ d ξ , (10) where W k , ξ is the 2D Fourier transformation of w x , η . U k and R ξ are the Fourier transformation of u x and r η , respectively. If r is the random reflectivity, we have (Yilmaz, 2001) U k = R ξ ∫ − ∞ + ∞ W k , ξ d ξ , (11) where R ξ is a constant value, which can be determined from the well-side seismic traces. We then perform the inverse Fourier transform to Equation 11 to reproduce the original signal, as follows (Equation 12): u x = ∫ − ∞ + ∞ U k e 2 π i k x d k = R ξ ∫ − ∞ + ∞ ∫ − ∞ + ∞ W k , ξ e 2 π i k x d ξ d k = R ξ ∫ − ∞ + ∞ w x , ξ d ξ . (12)

Combined with the original signal reconstruction method, the following expression is obtained: R ξ ∫ − ∞ + ∞ w x , ξ d ξ = ∫ − ∞ + ∞ S x , η d η , (13) where w x , ξ represents the depth-domain seismic wavelet and S x , η represents the results after inverse Fourier transform along the depth of the depth–wavenumber spectrum (e.g., S 2 ). Then, the estimation of the depth-domain seismic wavelet can be achieved based on the depth–wavenumber spectrum calculated by the method proposed in Section 2.1.

We first experiment with the proposed method using a simple three-wavelet model, referencing the model used in Zhang and Deng (2018). Figure 2A shows a 40-Hz zero-phase time-domain Ricker wavelet, and its corresponding amplitude spectrum is shown in Figure 2B. Taking this wavelet as the source, the corresponding depth-domain wavelets as they propagate through the different velocity strata are shown in Figure 2C, and their corresponding amplitude spectrums are shown in Figure 2D. It can be seen that as the velocity changes, the depth-domain wavelet is significantly stretched or compressed, i.e., the depth-domain wavelet is velocity-dependent. Since subsurface velocities are spatially variable, the interpretation and inversion of depth-domain seismic data can result in unreliable results if a constant wavelet is utilized.

Figure 3 shows the single-wavelet (i.e., sky-blue line shown in Figure 2C) depth–wavenumber spectrum experiments for different methods. It can be seen that WUST (Figure 3C) obtains a more accurate wavenumber spectrum than ST (Figure 3B) by removing the linear wavenumber-dependent term, but it performs poorly in the low-wavenumber region. By introducing parameters A and B, MWUST ensures that the wavenumber spectrum is accurate while performing well in the low-wavenumber region (Figures 3D–F). Figures 3D–F show that A plays a greater role in the depth–wavenumber spectrum than B, which is consistent with the conclusion reached in Figure 1. The amplitude spectra are extracted from the red dashed lines in Figure 3 to describe the wavenumber distribution in detail, as shown in Figure 4. Here, 1/3 and 10 are used for A and B in MWUST, respectively. It can be seen that the central frequency of the amplitude spectrum obtained based on the ST deviates from the reference frequency. However, MWUST matches the FT curve well, which further validates the advantages and effectiveness of the proposed method (i.e., MWUST) over ST and WUST.

Subsequently, a simple model in which the three reflection coefficients (Figure 5A) are related to the wavelet velocities mentioned in Figure 2C is used to perform wavelet estimation. Figures 6A–C show the depth–wavenumber spectrum of ST, WUST, and MWUST, where the black and red asterisks indicate the reference primary wavenumber of the signal at the position of the reflection coefficient and the primary wavenumber corresponding to the depth–wavenumber spectrum calculated by the different methods, respectively. It can be seen that the wavenumber distributions of WUST and MWUST match the reference values more closely than those of ST, but the depth–wavenumber spectrum of MWUST has a higher depth resolution in the low-wavenumber region (white arrows).

Figures 7A–C show the depth-variant wavelets extracted from the depth–wavenumber spectrum obtained using different methods. Based on the known reflection coefficients, the depth-domain wavelets extracted by the different methods are convolved with them to obtain reconstructed seismic records, as shown by the red line in Figures 5B–D. The seismic records reconstructed by the proposed method (i.e., MWUST) are best matched to the reference seismogram (black line in Figure 5D). Figure 8 shows the normalized errors of the reconstructed seismic records of the wavelet extracted by the different methods. The white dashed line indicates the exact solution, i.e., the closer the focus is to the white dashed line, the more accurate the reconstructed result is. As expected, the more accurate depth–wavenumber spectrum obtained by MWUST resulted in more accurate extracted depth-domain wavelets, thus minimizing the error between the reconstructed seismogram and the reference seismogram.

A more realistic example is further used to test the effectiveness of the method. The well–seismic ties between velocity curves and depth-domain seismic data (black lines) are displayed in Figure 9. Figures 10A–C show the depth–wavenumber spectrum of depth-domain seismic data using ST, WUST, and MWUST, respectively. The depth-variant wavelets extracted using the different methods are shown in Figures 11A–C. The reconstructed seismograms (red lines) obtained by combining known velocity well-log data and depth-domain wavelets extracted using the different methods are shown in Figures 9B–D. In this test, the density is assumed to be a constant density that does not vary with depth. The seismic records reconstructed (red lines) by the three methods roughly match the reference seismic records (black lines). However, the reconstructed seismic records of the proposed method are more consistent with the reference seismic records in some aspects compared to the reconstructed records of the other two methods (e.g., green arrows).

We perform further inversion tests to demonstrate the impact of the extracted wavelet on subsequent seismic inversion and interpretation. The algorithm used here is a Bayesian-based inversion method, the details of which can be found in Zhang et al. (2021a). In order to ensure a fair comparison of the inversion results, we perform the inversion by changing only the input wavelet, leaving the other inversion parameters unchanged. Figure 12 shows the inversion results obtained from the wavelet extracted by the different methods (i.e., ST, WUST, and MWUST) as input. It can be seen that the results obtained from the inversion of the wavelet extracted based on the ST show significant oscillations compared to those obtained from the inversion of the wavelets extracted by the other two methods. Although small differences are exhibited in the reconstructed seismic records (Figure 9B), the impact of the small errors on the inversion results is significant. The accuracy of the inversion results (green arrows in Figure 12) of the proposed method is greater due to the extraction of a more accurate depth-domain wavelet.