Reflection waveform inversion based on full-band seismic data reconstruction for salt structure inversion

SUMMARY

Salt structures are high potential targets for oil and gas exploration. However, large-scale salt domes with irregular surfaces pose significant challenges for velocity model building. For full waveform inversion, in the absence of a high-fidelity initial model, the success of the inversion depends on low-frequency seismic data, which are scarce in the exploration data sets. This paper presents a new idea to solve the problem of salt structure velocity modelling. First, we propose an envelope-based full-band seismic data reconstruction algorithm. The smoothness of envelope is used to segment the events in seismic data, and the phase independence of envelope is used for the identification of the seismic event's arrival-time to obtain the apparent reflection sequences of the subsurface. Full-band seismic data are obtained by convolving the apparent reflection sequence with full-band source. Window averaging function and threshold strategy are used to ensure the accuracy of seismic event segmentation and the stability of the algorithm when dealing with noisy data. Then the multiscale reflection waveform inversion based on reconstructed data is proposed for salt structure velocity building. The numerical experiment results of the Sigbee2A model demonstrate the performance of the inversion algorithm in the case where the seismic data lack low-frequency components and contain noise. The limitations of the algorithm have also been analysed and studied.

Inverse theory, Waveform inversion, Acoustic properties, Seismic tomography, Wave propagation

1 INTRODUCTION

Salt structures are promising hydrocarbon reservoirs. Many well-known hydrocarbon reservoirs, such as those in Gulf of Mexico, Angola, Gabon and Tupi, are related to salt structures (Lerche 2013). Seismic exploration requires high-resolution subsurface imaging, which requires high-quality velocity models. However, the velocity building of salt structure is challenging due to the large-scale salt domes that have strong velocity contrast with sedimentary environment. The ‘Migrate-Pick-Flood’ workflow, which combining migration iterations, manual picking of salt boundary and salt dome velocity filling, is the most commonly used technology in industry for salt structure velocity reconstruction (Dellinger et al.2017). However, this technology involves extensive human intervention and the quality of the results cannot be guaranteed. Full waveform inversion (FWI) uses the optimization algorithm to minimize the difference between the observed and simulated data to reconstruct the model parameters (Tarantola 1984; Brossier et al. 2009; Virieux & Operto 2009). Recently, there have been many attempts to build the salt structure velocity model using FWI: Level set method (Lewis et al. 2012), total variation (Esser et al. 2015; Kalita et al. 2019), optimal transport distance (Engquist et al. 2016; Métivier et al.2016a, b, c, 2018) and variance-based model interpolation (Ovcharenko et al. 2018a) have already made promising progress. However, in order to obtain the global optimal solution, these methods more or less depend on low-frequency information in seismic data, which is always missing in field data sets.

Several data reconstruction methods have been proposed to solve the problem of missing low-frequency information. Xie (2013) proposed converting high-frequency waveform into low-frequency waveform using linear phase approximation based on the property that the seismic data are weakly dispersed at low frequencies. Li & Demanet (2015, 2016) used a phase tracking method to synthesize seismic data with low-frequency components by solving a nonlinear least-squares optimization problem. Inspired by interference beat tone, Hu et al. proposed the beat-tone FWI that uses the interference between the two relatively high-frequency seismic data to retrieve the low-wavenumber information of the subsurface (Hu 2014; Hu et al. 2018). Machine learning has also been applied to the reconstruction of full-band seismic data (Jin et al. 2018; Sun & Demanet 2018; Ovcharenko et al. 2018b).

Seismic data reconstruction needs to solve two problems: The first one is the definition of seismic data composition. Different definitions will significantly affect the tangibility and effectiveness of data reconstruction algorithms. Seismic data contains amplitude and phase information. Since both amplitude and phase are time-dependent and instantaneous concepts, this means that the dimensions of the solution space are huge, especially when the amount of seismic data is large, the computational cost and complexity of the reconstruction algorithm will be enormous. Based on this definition, the huge computational cost of the phase and amplitude tracking method (Li & Demanet 2015, 2016) verifies this conclusion. Another problem is the data reconstruction method. The current data reconstruction algorithms are either computationally intensive or unstable and need strict preconditions for the data. For example, the data reconstruction algorithm proposed by Li & Demanet (2015, 2016) needs to assume that the seismic data is sparse and has continuity in both the space and frequency domains. De-convolution can be used to obtain the reflection sequence of the subsurface (Kail et al. 2012; Kazemi & Sacchi 2014). However, de-convolution operation requires seismic data containing low-frequency components, otherwise the accuracy cannot be ensured (Yang & Zhu 2018).

The framework of this paper is as follows: First, we define the seismic data as the convolution of source wavelets and the apparent reflection sequence inspired by the post-stack migration. Then based on the definition of seismic data and the characteristics of observation data on salt structure, an envelope-based seismic data reconstruction algorithm is proposed. We use reconstructed seismic data for multiscale FWI. To improve the subsalt illumination, we introduce the wavefield decomposition into the inversion. Then we propose a multiscale reflection waveform inversion based on reconstructed data (MSRIRD). We use the Sigsbee2A model to verify the antinoise ability of the algorithm and its insensitivity to low-frequency data. Finally, in the discussion section, we analyse the essential reasons and preconditions for the success of this proposed method, and the significance and influence of the method on the improvement of FWI. We also discussed the limitations of this method and the corresponding solutions to provide research directions for future work.

2 ENVELOPE-BASED SEISMIC DATA RECONSTRUCTION

Low-frequency seismic data reconstruction can be regarded as an inverse problem with existing seismic data as a priori condition. The ill-posedness of the inverse problem is the main obstacle to obtain the global optimal solution. When thinking of how to reduce the multisolution of the inverse problem, compressed sensing in signal processing is a good inspiration. Compressed sensing reduces the ill-posedness of the inverse problem by reducing the number of parameters that need to be solved by sparse constraints (Donoho 2006), which has been widely used in the geophysical data processing (Herrmann & Hennenfent 2008; Chen et al. 2013; Zhou et al. 2016). Inspired by the principle of compressed sensing, we can seek sparse representation of seismic data to reduce the computational cost and ill-posedness of the data reconstruction.

Convolution mode in post-stack data attribute analysis can be used as a sparse representation of seismic data:

$$\begin{eqnarray*} u({x_r},{x_S},t) = s(t)*G({x_r},{x_S},t) = s(t)*\sum\limits_i {{\gamma _i}\delta (t - {t_i}} ), \end{eqnarray*}$$

(1)

where s(t) is the source wavelet, G(x_r, x_s, t) is the Green's function, which represents the unit impulse response of the subsurface medium, or the reflection sequence |${\gamma _i}$|of the medium. |$\delta (t)$| is a Dirac delta function and * represents the convolution operation. If the source wavelet is known, the problem of waveform matching between the observed and simulated data is simplified to the matching problem of the reflection sequence. The number of parameters in the inverse problem will be greatly reduced, especially when the subsurface is a salt structure, huge and relatively homogeneous salt domes generate less internal reflections, which causes the reflection events in seismic data to be sparse.

Similar to post-stack data, pre-stack seismic data can also be expressed using the convolution mode in case where the earth is regarded as a linear system. Each event in the pre-stack data can be regarded as the reflection of the source wavelet on a certain interface, and the time that the event is generated determines when the reflection wavelet reaches the reflection interface, that is, the traveltime, corresponding to the reflection interface position. We define an apparent reflection sequence as the traveltime mark of the source wavelet, and the magnitude is proportional to the interface reflectivity. The pre-stack data are the convolution of the source wavelet and the apparent reflection sequence:

$$\begin{eqnarray*} u({x_r},{x_S},t) = s(t)*\sum\limits_i {{{\tilde{\gamma }}_i}\delta (t - {t_i}} ), \end{eqnarray*}$$

(2)

where |$\tilde{\gamma }$| is the apparent reflection sequence.

The reconstruction of the apparent reflection sequence needs to solve two problems, namely the segmentation of seismic events and the identification of arrival-times of individual reflection event. For the former, the conventional approach is to use a window to manually divide the seismic data, rather than segmenting based on the actual start and end times of the events in the seismic data. For the latter, we should ensure that the criteria for the arrival time of each event are uniform and the arrival time mark of each event based on this criterion is unique. However, in a generally heterogeneous medium, complicated wave phenomena such as reflection, refraction or diffraction could happen, and the phase of the source will shift at the reflective interface. The determination of the arrival-time of each event in a constantly changing waveform is a challenge.

Envelope was first applied for picking seismic events in earthquake seismology (Baer & Kradolfer 1987). In view of the linear behaviour of envelope in the subsurface parameters, envelope was used to construct an objective function in envelope inversion (Bozdağ et al. 2011; Wu et al. 2014; Yuan et al. 2015; Gao et al. 2017). In order to reduce the artificial noise generated by the envelope and improve the convergence of the inversion algorithm, the auxiliary equations were introduced into the objective function of envelope inversion (Bharadwaj et al. 2016; Liu & Zhang 2017). Envelope and global correlation norm were combined together to simultaneously take advantage of envelope and global correlation norm to further reduce the possibility of FWI occurring cycle skipping (Oh & Alkhalifah 2018; Wu et al. 2019). In order to solve the problem caused by the use of the waveform Fréchet derivative in envelope inversion, multiscale direct envelope inversion based on a direct envelope Fréchet derivative was proposed for salt structure inversion (Chen et al.2017, 2018a, 2019). In addition to the above mentioned applications in inversion, envelope also has two properties that can be used to solve the problem of this paper.

Given seismic data |$f(t)$|⁠, envelope |$e(t)$|of |$f(t)$| is:

$$\begin{eqnarray*} e(t) = \sqrt {{f^2}(t) + f_H^2(t)} , \end{eqnarray*}$$

(3)

where H is the Hilbert transform. As can be seen from eq. (3), envelope can be considered as the instantaneous energy. We define a complete seismic event as a process that begins with the wavelet's wavefront reaching the reflection interface and ends with the entire wavelet passing the reflection interface. Correspondingly, the envelope's amplitude of a seismic event starts from zero and then reaches the peak and finally returns to zero. Picking the start and end points of events from the intricate seismic data involves extensive human intervention and is difficult to control quality due to the seismic data's oscillation. But for envelope, the difficulty of segmenting the events will be greatly reduced. First, envelope is an energy flow without positive and negative polarity, which is smoother than seismic data. Secondly, as analysed above, the start and end points of an event correspond to envelope are the start and end points of energy, according to the extreme value theorem, that is, the minima point of envelope. The other case is where there are overlapping parts of the two events. However, due to the smooth property of envelope, the inconsistency between different events will cause envelope to be discontinuous. It is also possible to use the minima point of envelope to separate the two events. In this way, the separation of the event becomes the calculation of the minima point of envelope, which is a simple problem that does not require human intervention.

After the segmentation of the events, the next step is to identify the event's traveltime. Another property of envelope helps us solve this problem. For seismic data, we can represent it as

$$\begin{eqnarray*} f(t) = e(t) * cos(\omega t + \varphi ). \end{eqnarray*}$$

(4)

Where wis frequency and |$\varphi $|is the initial phase. It can be seen from eq. (4) that envelope is independent of the seismic data phase, in other words, the phase variation of the source wavelet during propagation does not change the amplitude of envelope. In Fig. 1, we give the zero-phase Ricker wavelet and the wavelets with different phase shift. Although the phase of these wavelets are different, their envelope are the same. In addition to the Ricker wavelet, other types of wavelets such as Morlet, Kaluder wavelets also have such property. In fact, all wavelets whose phase and amplitude information are independent of each other have such property. This property of envelope can be used to identify the reflector. To guarantee uniqueness, we choose the maxima point of envelope as the arrival time mark. After that, the amplitude at the maxima point of envelope is used as the reflectivity of the reflector. In addition, the importance of the polarity information has been verified in signed envelope inversion (Chen et al. 2018b). We use the phase threshold to define the polarity of the reflection sequence, that is, compare with the original source wavelet, if the phase shift of the seismic data is less than 180°, it is positive reflection, and if it is greater than 180°, it is negative reflection, corresponding to the positive and negative of the reflection coefficient.

Figure 1.

Envelope's amplitude is independent of phase. Zero phase Ricker wavelet(blue line); Ricker wavelet with pi/4 phase shift (cyan-blue line); Pi/2 phase shift (pink line); 3*Pi/4 phase shift (red line); Pi phase shift (black line); green line is the envelope of the wavelets.

The method mentioned above is suitable for the case where envelope of the source wavelet is good in continuity and the seismic data do not contain noise. There are also multiple discontinuities in envelope of one single event due to the source wavelet being too turbulent or the seismic data containing noise. In this paper, we use the window function to obtain multiscale envelope to deal with the problem of source wavelet discontinuity.

$$\begin{eqnarray*} {e_L}(t) = \frac{1}{L}\int_{{ - L/2}}^{{L/2}}{{dt{\prime}W(t - t^{\prime})[e(t^{\prime})]}}, \end{eqnarray*}$$

(5)

where W(t) is the window averaging function and L is the window length. The |${e_L}$| indicates the window-averaged envelope (WAE). Window averaging function is a linear filter operator, which will not change the position of envelope's maxima point. In order to demonstrate the effect of window averaging function on increasing envelope smoothness, we apply a bandpass filter with a frequency band of 4–20 Hz to the Ricker wavelet to obtain a band-pass wavelet (Fig. 2). Because of the lack of low-frequency components, high-frequency oscillations occur in the source wavelet, resulting in multiple local minima points in the envelope, which hinders the segmentation of seismic events. Compared with the conventional envelope, the smoothness of the WAE is significantly improved, which meets the requirements of envelope smoothness for events segmentation. When the seismic data contains noise, the noise will create false events in places that do not contain effective seismic events. In order to solve the problem caused by noise, we set amplitude threshold for the WAE |${e_L}$|to eliminate false events. Table 1 shows the detailed workflow of the data reconstruction algorithm. The choice of window length and amplitude threshold depends on the continuity of the source wavelet and the signal-to-noise ratio [SNR = log₁₀ (signal power/noise power)\ of the seismic data, so there is no strict formula to follow. A basic principle is that the window width and amplitude thresholds should be as small as possible when the conditions for accurate segmentation of events are met. An oversized window and amplitude threshold will oversmooth the seismic data and neglect some useful reflections with smaller amplitude, resulting in data reconstruction errors.

Figure 2.

The effect of window averaging function on increasing envelope smoothness. Red line is the source wavelet after being filtered by the band-pass filter, which has a frequency band of 4–20 Hz. Black triangle are the local minima points of envelope. (a) Conventional envelope of the source wavelet (green line); (b) The window-averaged envelope (WAE) of the source wavelet for window length L = 51(blue line). The conventional envelope has eight minima points, while the WAE has only two minima points.

Table 1.

The workflow of envelope-based data reconstruction algorithm.

(1)	Calculate envelope e of the seismic data u, of which the total recorded time is T.
(2)	Set up the window length Laccording to continuity of the source S and the SNR of seismic data, calculate the window-averaged envelope|${e_L}$|⁠.
(3)	Calculate the local minima point coordinates of the window-averaged envelope |${P_k}(k = 2,3.....m - 2,m - 1)$|⁠. |${P_1}$|is the starting point and |${P_m}$| is the end point of the seismic data.
(4)	Set threshold|$\varepsilon {\rm{ = }}\lambda \sum\limits_{it = {P_1}}^{{P_m}} {{e_L}(it)} $|⁠, where|$\lambda $| is the threshold control parameter and|$0 < \lambda < 1$|⁠, if|$\sum\limits_{it = {P_k}}^{{P_{k + 1}}} {{e_L}(it)} < \varepsilon $|⁠, remove|${P_{k + 1}}$|from|$\{ {{P_1}} ,P{}_2......{P_{m - 1}}, {{P_m}} \}$|⁠.
(5)	Define|${u^k}(it) = u(it),e_L^k(it) = {e_L}(it)(it = {P_k},{P_k} + 1,{P_k} + 2,...{P_{k + 1}} - 1)$|⁠. Define the maxima point coordinates of the absolute value of |$e_L^k$| as Pmaxk.
(6)	Initialize the apparent reflection sequence|$\tilde{\gamma }(it)(it = 1,2,3...T)$|⁠. Define the phase shift between the event |${u^k}$| and the source wavelet S as|$PhaseS$|⁠. $$\begin{eqnarray*} \tilde{\gamma }({P_{{{\max }_k}}}) = \{ \begin{array}{@{}l@{}} {e_L}({P_{{{\max }_k}}})forPhaseS \lt \pi \\ - {e_L}({P_{{{\max }_k}}})forPhaseS \gt \pi \end{array} \end{eqnarray*}$$
(7)	The reconstructed seismic data |$\tilde{u} = \tilde{\gamma } * \tilde{S}$|⁠, where |$\tilde{S}$|is full-band source.

In order to facilitate the understanding of the workflow above, we present a schematic diagram (Figs 3a–d) showing the results of reconstructing seismic data in different cases. The properties of envelope, window averaging function and threshold strategy make the envelope-based data reconstruction algorithm can deal with situations where seismic data lacks low frequencies, phase shifts in the source and seismic data contain noise. Figs 4(a–d) shows the reconstruction results of four traces synthetic data that are obtained from Sigsbee2A (Fig. 6a). We add random noise to the data so that the SNR of the seismic data is 1. Window length Land threshold parameter |$\lambda $| are set to 50 and 5 per cent, respectively. Fig. 5 shows the shot profile of the original noisy data (Fig. 5a) and reconstructed data (Fig. 5b). Comparing the two shot profiles shows that the reconstruction method can preserve local coherency along the offset axis. The results in Figs 4 and 5 demonstrate the stable performance of envelope-based reconstruction algorithm.

Figure 3.

The results of seismic data reconstruction in four different cases: (a) Source with full-frequency band; (b) Source without low-frequency components below 4 Hz; (c) Source without low-frequency components below 4 Hz and has a phase shift of 45°; (d) Source without low-frequency components below 4 Hz and has a phase shift of 45° and contains random noise, the signal-to-noise ratio (SNR) is 1.Red line is original seismic data, green line is envelope and blue line is reconstructed data. Black triangle and blue star are the local minima and maxima points of envelope, respectively.

Figure 4.

(a-d) Reconstruction results of four traces synthetic data that are obtained from the Sigsbee2A model. The seismic data contains random noise, which the SNR is 1. The frequency components below 4 Hz are truncated.

Figure 5.

(a) The shot profile of original noisy data, whose SNR is 1; (b) The shot profile of reconstructed data.

3 MULTISCALE REFLECTION WAVEFORM INVERSION BASED ON RECONSTRUCTED DATA (MSRIRD)

Multiscale FWI can be achieved after obtaining reconstructed seismic data with low-frequency components. In order to avoid the truncation effect caused by linear filtering, we implement a multiscale inversion strategy (Bunks et al. 1995) that uses reconstructed seismic data with different dominant frequencies. We convolve the reconstructed apparent reflection sequence with the source wavelets with different dominant frequencies. It should be noted that when we reconstruct the observed data, the phase of the observed data are adjusted, so the synthetic data should be reconstructed in the same way to ensure that the observed and synthetic data can be perfectly matched. The objective function of multiscale FWI based on reconstructed data (MSIRD) can be written as

$$\begin{eqnarray*} \sigma = \sum\limits_S {\sum\limits_R {\sum\limits_T {\left( {{{\tilde{u}}_{obs}} - {{\tilde{u}}_{syn}}} \right)} } } , \end{eqnarray*}$$

(6)

where |${\tilde{u}_{obs}},{\tilde{u}_{syn}}$| are the reconstructed observed and synthetic data, respectively. σ is the summed residual between |${\tilde{u}_{obs}}$|and |${\tilde{u}_{syn}}$|with respect to the source S, the receiver R and the total recording time T. Using the adjoint state method, the gradient of the multiscale inversion can be calculated by:

$$\begin{eqnarray*} g = \sum\limits_{SRT} {\frac{2}{{{\bf v}_0^3(x,z)}}\frac{{{\partial ^2}{{\left( {{{{\bf \tilde{u}}}_S}\left( {x,z,t} \right)} \right)}_W}}}{{{\partial ^2}t}}{G^*}({{\tilde{u}}_{obs}} - {{\tilde{u}}_{syn}})} , \end{eqnarray*}$$

(7)

where |${{\bf \tilde{u}}_S}( {x,z,t} )$|is the reconstructed source wavefield. G is a Green's operator and |${G^*}$| is the conjugate transformation ofG, which represents the wavefield backpropagation, and we can calculate |${G^*}({\tilde{u}_{obs}} - {\tilde{u}_{syn}})$| to obtain the reconstructed residual wavefield |${{\bf \tilde{u}}_R}( {x,z,t} )$|⁠.

In addition to the velocity inversion of the salt dome, the velocity of subsalt structure is critical to the determination of the reservoir position. However, the reflected energy from subsalt region is weak due to the occlusion of the large-scale salt domes with complex shapes. That is to say, the low-wavenumber components that can reflect the subsalt structural information in the model gradient are few. Reflection waveform inversion (RWI, Xu et al. 2012) provides a good solution to this problem. This method decomposes the short-wavelength and long-wavelength components of the gradient by performing wavefield decomposition on the source and residual wavefield. We introduce this method into MSIRD. We use the 2-D fast Fourier transform (2-D FFT, Liu et al.2011) to separate the reconstructed source wavefield |${{\bf \tilde{u}}_S}$| and residual wavefield |${{\bf \tilde{u}}_R}$| into up- and downgoing parts:

$$\begin{eqnarray*} \left\{ \begin{array}{@{}l@{}} {{{\bf \tilde{u}}}_S} = {\bf \tilde{u}}_S^ + + {\bf \tilde{u}}_S^ - \\ {{{\bf \tilde{u}}}_R} = {\bf \tilde{u}}_R^ + + {\bf \tilde{u}}_R^ - \end{array} \right., \end{eqnarray*}$$

(8)

where |${\bf \tilde{u}}_S^ + $|and |${\bf \tilde{u}}_S^ - $| are the direct and reflected forward wavefield, respectively. |${\bf \tilde{u}}_R^ + $| is the direct residual wavefield and |${\bf \tilde{u}}_R^ - $| is the reflected residual wavefield. When the source and residual wavefield propagate in opposite directions, the cross-correlation between them can obtain the high-wavenumber update of the model, otherwise the low-wavenumber update is obtained (Tang & Lee 2013; Yao et al. 2018). We define |${g_h}$|and|${g_l}$| as the high and low-wavenumber updating gradients, respectively:

$$\begin{eqnarray*} \begin{array}{@{}l@{}} {g_h} = {\bf \tilde{u}}_s^ + * {\bf \tilde{u}}_R^ - + {\bf \tilde{u}}_s^ - * {\bf \tilde{u}}_R^ + .\\ {g_l} = {\bf \tilde{u}}_s^ + * {\bf \tilde{u}}_R^ + + {\bf \tilde{u}}_s^ - * {\bf \tilde{u}}_R^ - . \end{array} \end{eqnarray*}$$

(9)

The workflow of MSRIRD is to use the reconstructed data with low dominant frequency to obtain the forward wavefield and residual wavefield according to eq. (7), and then to decompose the gradient into low- and high-wavenumber components according to eq. (9). Then the alternating iteration method (Irabor & Warner 2016) is used to update the model. We use |${g_h}$| to generate reflectors, then use |${g_l}$| to update the subsurface background velocity. When the long-wavelength components of the model have been well reconstructed, for example, when the iteration meets the iteration termination condition:

$$\begin{eqnarray*} \frac{{\left| {{\sigma _n} - {\sigma _{(n - 1)}}} \right|}}{{\left| {{\sigma _{(n - 1)}}} \right|}} < \eta , \end{eqnarray*}$$

(10)

where |${\sigma _n}$| is the least square error of the misfit function for the nth iteration,|$\eta $| is a positive number far less than 1.0. We gradually increase the dominant frequency of the reconstructed data to reconstruct the short-wavelength components of the model in the same way. Finally, we use the conventional FWI based on the original seismic data to correct the inversion error in high-wavenumber components, which is caused by the difference between the reconstructed and true seismic data. Table 2 is the workflow of MSRIRD.

Table 2.

The workflow of multiscale reflection waveform inversion based on reconstructed data.

Step1: Set up an initial model |${m_0}$|

Step2: Set up the dominant frequency f of the reconstructed data. If |$f \to \tilde{f}$|⁠, where |$\tilde{f}$| is the dominant frequency of the original seismic data, go to Step8.

Step3: Applying the envelope-based data reconstruction algorithm to the original source wavefield |${{\bf u}_s}$| to obtain the reconstructed source wavefield |${{\bf \tilde{u}}_s}$|with the dominant frequencyf.

Step4: Applying the envelope-based data reconstruction algorithm to the original observed data|${u_{obs}}$|and synthetic data |${u_{syn}}$|to obtain the reconstructed observed data|${\tilde{u}_{obs}}$|and synthetic data|${\tilde{u}_{syn}}$|⁠, which dominant frequency isf. Then use |${G^*}({\tilde{u}_{obs}} - {\tilde{u}_{syn}})$| to calculate residual wavefield |${{\bf \tilde{u}}_R}$|⁠,.

Step5: Calculate the high-wavenumber gradient|${g_h}$| to update the model's high-wavenumber components.

Step6: Use the model that has been updated in Step 5 as the initial model, calculate the low-wavenumber gradient|${g_l}$| to update the model's low-wavenumber components.

Step7: If the iteration termination condition is satisfied, go to Step 2, or go to Step 3, repeat all the iterations represented by Steps 3 to 6.

Step8: Do FWI using original seismic data to obtain the final inversion result m.

4 NUMERICAL EXPERIMENTS ON 2-D Sigsbee2A MODEL

As the numerical example for the proposed inversion method, we use the Sigsbee2A benchmark test model (Fig. 6a), which is a typical salt structure (Stoughton et al. 2001). One of the major difficulties in the inversion for this model is the large-scale salt dome, which extends for 18 km horizontally and nearly 4 km vertically. In addition, the shape of the salt dome is complicated, which disturbs the observation of reflected waves from the subsalt. Another difficulty is the inversion of the horizontal salt layer in the model bottom. To accurately reconstruct these deep structures, high-fidelity velocity in the upper part of the model is required. The difficulties mentioned above are typical in current salt structure inversions. To start the inversion, we use a 1-D initial velocity model with velocity increasing linearly with depth except for the water portion, which remains the same as the true model (Fig. 6b). An 8 Hz Ricker wavelet is used as the source and the frequency components below 4 Hz are truncated. The regular-grid finite-difference method is used to do forward modelling in time domain and perfect match layer method as boundary condition for boundary reflection. The misfit function is minimized by the steepest descent method and the adjoint state method is used to calculate the gradient. 107 shots are evenly distribute with a horizontal interval of 228.6 m, and each shot is recorded by 640 receivers located every 38.1 m.

Figure 6.

(a) Sigsbee2A model; (b) initial model.

As a comparison, we give the inversion result of conventional FWI (Fig. 7a) after 200 iterations. As expected, conventional FWI is most sensitive to the sharp upper boundary of the salt dome. Shown in Fig. 7(b) is the inversion result using RWI after 150 iterations, while the inversion depth is better than conventional FWI, in the absence of low-frequency information, it still falls into local minima. Fig. 8 is the inversion result of MSIRD. In MSIRD, the reconstructed data with dominant frequencies of 0.5, 2 and 4 Hz are used successively and the number of iterations are 20, 35 and 50, respectively. The inversion result of the previous stage is taken as the initial model for the inversion of the next stage and the inversion results are shown in Figs 8(a)–(c), respectively. Then the original seismic data are used to retrieve the high-wavenumber components of the model (Fig. 8d). It can be seen from the inversion results that the accuracy of the reconstructed data can satisfy the inversion of the salt dome. Not only the velocity of the salt dome, but also the top and bottom boundaries, are close to the true model. One imperfection that needs to be improved is the subsalt structure in the large dip area, and artificial noise is generated, which also leads to the failure of the horizontal salt layer reconstruction.

Figure 7.

Inversion result of (a) Conventional FWI; (b) RWI.

Figure 8.

(a)–(c): Inversion results of different stages of MSIRD. Reconstructed data that dominant frequency is (a) 0.5 Hz; (b) 2 Hz; (c) 4 Hz; (d) inversion result of conventional FWI using original seismic data, the initial model is Fig. 8(c).

To demonstrate the effect of the RWI on the inversion results, we use the same parameters as the previous experiment, and the corresponding inversion results of MSRIRD data are shown in Figs 9 (a)–(c) are the inversion results using the reconstructed data with dominant frequencies of 0.5, 2 and 4 Hz, respectively. Due to the contribution of the reflection waveform inversion strategy, the convergence rate of the inversion is accelerated, and the number of iterations in the three stages are 10, 15 and 20, respectively. Fig. 9(d) is the finial inversion result using conventional FWI after 80 iterations, which the initial model is Fig. 9(c). After the introduction of the RWI, the inversion accuracy of subsalt structure in the steep dip area is improved, and the inversion effect of the horizontal salt layer at the bottom is also improved, although there is still large deviation from the true model (Fig. 6a).

Figure 9.

(a)–(c): Inversion results of different stages of MSRIRD. Reconstructed data that dominant frequency is (a) 0.5 Hz; (b) 2 Hz; (c) 4 Hz; (d) inversion result of conventional FWI using original seismic data, the initial model is Fig. 9(c).

To show the inversion results difference among the different inversion methods, we select two velocity curves from the model that are located at 12.9 km (Fig. 10a) and 20.5 km (Fig. 10b) from the left side of the model. True velocity, initial velocity and the inversion results of FWI, RWI, MSIRD and MSRIRD are shown in one figure. Comparing the velocity curves of FWI and RWI, it can be seen that RWI is better than FWI in depth and fidelity for the salt dome inversion. However, RWI does not reach the inversion requirement of salt domes in the absence of low-frequency components. Both MSIRD and MSRIRD can meet the fidelity requirements for the salt domes inversion. The main difference between these two methods is the inversion accuracy of the subsalt velocity. The subsalt inversion error in MSIRD is well corrected after using the RWI strategy in MSRIRD. Fig. 11 shows the reverse time migration (RTM) results of the true velocity model and inversion results of conventional FWI, RWI, MSIRD and MSRIRD. These RTM images verify the difference in accuracy of the inversion results.

Figure 10.

Comparison of inversion results among the FWI, RWI, MSIRD and MSRIRD. Velocity of true model (cyan-blue line), initial model (black line), FWI (red line), RWI (green line), MSIRD (pink line) and MSRIRD (blue line): (a) Velocity profiles located at the distance 12.9 km; (b) Velocity profiles located at the distance 20.5 km.

Figure 11.

Migration results using different velocity model: (a) True model; (b) Inversion result of Conventional FWI; (c) Conventional RWI; (d) MSIRD + FWI; (e) MSRIRD + FWI.

Finally, we add random noise to the observed data to generate noisy data to test the antinoise ability of MSRIRD. Fig. 12(a) is one shot profile from the Sigsbee2A model that without noise, Figs 12(b) and (c) are noisy data whose SNR is 5 and 1, respectively. In order to deal with the discontinuity caused by noise to seismic data, we set the window length L as 25 and 50, and the threshold parameter |$\lambda $| is set as 2.5 and 5 per cent, respectively. Figs 13(a) and (b) show the inversion results using data whose SNR are 5 and 1, respectively.With the help of the window averaging function and the threshold strategy, there is not much difference in the inversion results using noisy data compared to noise-free data.

Figure 12.

One shot profile from the Sigsbee2A model: (a) without noise; (b) noisy data, which SNR is 5; (c) noisy data, which SNR is 1.

Figure 13.

The inversion results of MSRIRD using noisy data (a) SNR is 5;(b) SNR is 1.

5 DISCUSSION

The method proposed in this paper is essentially an inversion method based on traveltime information, which is not affected by the frequency band of seismic data. We use envelope to identify and locate the traveltime nodes of seismic waves to obtain the apparent reflection sequence of the subsurface. Full-band seismic data are obtained by convolving the apparent reflection sequence with full-band source wavelet, which avoids both the high-frequency approximation and waveform matching that needs to be solved in conventional traveltime inversion. Another factor that contributes to the success of this method is the sparsity of the seismic data obtained from the salt structure due to huge and relatively homogeneous salt domes generate less internal reflections. Sparseness of reflected data causes the data misfit function to be insensitive to the velocity changes in salt dome, which is a serious problem for conventional FWI in salt structure inversion. However, for the envelope-based data reconstruction algorithm, it is a favourable condition.

In salt structure, in addition to salt dome, there are also several small scale scatterers or thin layer structures, which will introduce overlapping events into the seismic records. The method in this article does not guarantee complete separation of two events when the time interval between them is too small. In Figs 14 and 15, we use the Ricker wavelet as an example to show the situation of two overlapping reflection waves generated by a thin layer or a small reflector. Due to the space interval between the two waves are small, we do not need to consider the amplitude attenuation and can assume that the amplitude of two reflection waves are same. So the situation can be simply reduced into two cases: two adjacent reflection waves are either positive reflections (or negative) or one is positive and the other is negative reflection. Fig. 14 shows two adjacent reflection waves that are both positive reflections. In Figs 14(a)–(c), the time interval between two events is 3T/8, T/4, T/8, respectively. T is the period of Ricker wavelet. It can be seen from Fig. 14(c) that when the time interval between two events is less than T/8, the two events cannot be segmented by the method of this paper. In Fig. 15 we show two adjacent reflection waves that one is positive and the other is negative reflections. In this case, T/4 is the minimum time interval at which the segmentation method can work. When the subsurface contains a large number of thin layers or small scatterers, the overlapping events are the dominant components in the seismic data, which is a serious challenge for envelope-based data reconstruction algorithms. Fig. 16 is one shot profile from the SEG/EAGE Overthrust model and the source is Ricker wavelet that the dominant frequency is 8 Hz. The source is located 5 km away from the left end of the model. A comparison of the waveform and spectrum between the reconstructed and original seismic data is shown in Figs 17(a) and (b), respectively. When the seismic data interval does not meet the conditions shown in Figs 14 and 15, the waveform and spectrum of the reconstructed data will be greatly deviated from the true data. Correspondingly, errors generated in seismic data reconstruction can affect the accuracy of inversions, even affect the reconstruction of salt dome. For example, in the Sigsbee2B model (Fig. 18), multiple becomes the dominant components of seismic data because of a more realistic water to floor boundary. The multiple affects the inversion fidelity, especially the bottom of the salt dome. Fig. 19 shows the inversion results of the Sigsbe2B model for different stages of MSRIRD, which uses the same parameters as the numerical experiments on the Sigsbee2A model. Comparing the inversion results of the Sisgsbee2A model (Fig. 9), the importance of seismic data reconstruction accuracy for inversion is verified.

Figure 14.

Two adjacent reflection waves that are both positive reflections. The time interval between two events: (a) 3T/8; (b) T/4; (c) T/8. Red line is the seismic data and green line is the envelope, blue pulses are reflection sequences and black triangle are the local minima points of envelope.

Figure 15.

Two adjacent reflected waves that have opposite polarity. The time interval between two events: (a) 3T/4; (b) T/2; (c) T/4. Red line is the seismic data and green line is the envelope, blue pulses are reflection sequences and black triangle are the local minima points of envelope.

Figure 16.

One shot profile from the SEG/EAGE Overthrust model. The source is located 5 km away from the left end of the model.

Figure 17.

(a) The 760th trace seismic data (red line) and the corresponding reconstructed data (blue line). (b) The spectra of the seismic data (red line) and reconstructed data (blue line).

Figure 18.

Sigsbee2B velocity model.

Figure 19.

(a)–(c): Inversion results of different stages of MSRIRD on Sigsbee2B model. Reconstructed data that dominant frequency is (a) 0.5 Hz; (b) 2 Hz; (c) 4 Hz; (d) Inversion result of conventional FWI using original seismic data, the initial model is Fig. 19(c).

As discussed above, the effectiveness of the envelope-based data reconstruction algorithm is limited by the sparseness of the seismic data. The sparseness, in other words, is the resolution of the seismic data. The higher the resolution of the original seismic data, the higher the fidelity of the data reconstruction algorithm. Correspondingly, the inversion results based on the reconstructed data are guaranteed. Figs 20(a)–(c) show the inversion results of MSRIRD on Sigsbee2B model using dominant frequencies 4, 8 and 16 Hz, respectively. Comparing the three inversion results, it can be seen that the inversion quality is getting higher with the improvement of the resolution of the seismic data. When the dominant frequency of the seismic data reaches 16 Hz, the inversion accuracy of Sigsbee2B is similar to that of Sigsbee2A. The experimental results in Fig. 20 indicate that for the envelope-based data reconstruction algorithm, the demand for high resolution seismic data with high-fidelity is more urgent than low-frequency seismic data. Studying the data processing method to obtain the high resolution seismic data is the most effective way to further improve the accuracy of the inversion algorithm that proposed in this paper. Deconvolution method based on sparse constraints is a commonly used method to improve the resolution of seismic data (Taylor et al. 1979; Sacchi 1997; Kail et al. 2012). Chen & Wang (2018) proposed a wavelet compression method that utilizes the scale characteristic in the Fourier transform to extend the high- and low-frequency band simultaneously. For visco-acoustic cases, inverse Q-filters are used to enhance the seismic data resolution (Wang 2002, 2006, 2008; Cavalca et al. 2011).

Figure 20.

Inversion results of MSRIRD on Sigsbee2B model using different dominant frequencies: (a) 4 Hz; (b) 8 Hz; (c) 16 Hz.

In addition to reducing the dependence of inversion on low-frequency seismic data, the proposed method has another potential advantage. In this paper the source is assumed to be known. In the case where the source wavelet is unknown, the estimation of the source wavelet is also an important problem. If the phase of the estimated source wavelet is different from the true source wavelet, the accuracy of the inversion will be severely affected. However, the envelope-based data reconstruction algorithm in this paper is not affected by the phase of the source wavelet. That is to say, as long as the source wavelet is roughly estimated and the appropriate window length and threshold parameter are set according to the requirements of the event segmentation, we can select any suitable wavelet as the source and convolve with the apparent reflection sequence to obtain full-band seismic data. So the inversion algorithm proposed in this paper can be regarded as a source-independent inversion method.

6 CONCLUSIONS

This study demonstrates a multiscale reflection waveform inversion based on reconstructed seismic data for salt structure velocity building. The phase independence and smoothness of envelope and the sparsity of reflection seismic data from salt structures guarantee the success of data reconstruction. Window averaging function and threshold strategy guarantee the stability and antinoise ability of the data reconstruction algorithm. Reflection waveform inversion strategy provides a good solution to solve the inversion problem of subsalt structure in the steep dip area. The numerical results of the synthetic data set for the Sigsbee2A model demonstrate the accuracy of the reconstructed seismic data, as well as the robustness of the algorithm. The accuracy in reconstructing seismic data plays a critical role in the velocity inversion.

ACKNOWLEDGEMENTS

This work is financially supported by the National Natural Science Foundation of China (Grant 41904100, 41874133 and 41674123), China Postdoctoral Science Foundation funded project (2018M640559). The authors would like to thank Prof W. Hu and other two anonymous reviewers for their fruitful suggestions and comments. The author also expresses sincere gratitude to the editor, Prof Jean Virieux for coordinating the review of this paper and approving it for publication.