Incorporating attenuation effects into frequency-domain full waveform inversion from zero-offset VSP data from the Stokes equation

For seismic modelling, the Stokes equation can be obtained by adding the absorbing term to the classical inhomogeneous elastic equation (Ricker 1977) in the time domain expressed as

$$\begin{eqnarray} &&\fl\rho ( z )\frac{{\partial ^2 u( {z,t} )}}{{\partial t^2 }} = \frac{\partial }{{\partial z}}\left[ {\rho ( z )\upsilon ^2 ( z )\frac{{\partial u( {z,t} )}}{{\partial z}}} \right]\nonumber\\ &&+ \,\frac{\partial }{{\partial z}}\left[ {\eta ( z )\frac{\partial }{{\partial z}}\left( {\frac{{\partial u( {z,t} )}}{{\partial t}}} \right)} \right] + f( {z,t} ). \end{eqnarray} \tag{ 1 }$$

Considering the first-order hyperbolic system with the perfectly matched layer (PML) boundary conditions, we can obtain

$$\begin{eqnarray} &&\fl\frac{{\partial u( {z,t} )}}{{\partial t}} + \gamma ( z )u( {z,t} ) = \frac{1}{{\rho ( z )}}\frac{{\partial q( {z,t} )}}{{\partial z}} + f(z,t), \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} &&\fl\frac{{\partial q( {z,t} )}}{{\partial t}} + \gamma ( z )q( {z,t} ) = \rho ( z )\upsilon ^2 (z)\frac{{\partial u( {z,t} )}}{{\partial z}}\nonumber\\ &&+\, \eta ( z )\frac{\partial }{{\partial z}}\left( {\frac{{\partial u( {z,t} )}}{{\partial t}}} \right). \end{eqnarray} \tag{ 3 }$$

where u(z, t) is the time domain wavefield and f is the time domain source function; ρ(z), υ(z), and η(z) are the density, velocity, and viscosity coefficient, respectively; γ(z) is the damping function in the PML boundary conditions; and q(z, t) is the stress.

Then, we transform (2) and (3) into frequency domains by Fourier transformation, and after some calculations we can obtain

$$\begin{equation} {\mathop{\rm i}\nolimits} \omega \xi ( z )U( {z,\omega } ) = \frac{1}{{\rho ( z )}}\frac{{\partial Q( {z,\omega } )}}{{\partial z}} + F(z,\omega ), \end{equation} \tag{ 4 }$$

$$\begin{equation} Q( {z,\omega } ) = \frac{1}{{{\mathop{\rm i}\nolimits} \omega }}\frac{1}{{\xi ( z )}}( {\rho ( z )\upsilon ^2 (z){\rm + }{\mathop{\rm i}\nolimits} \omega \eta ( z )} )\frac{{\partial U( {z,\omega } )}}{{\partial z}}. \end{equation} \tag{ 5 }$$

where $\xi ( z ) = 1 + \frac{{\gamma ( z )}}{{{\mathop{\rm i}\nolimits} \omega }}$, ω is the angular frequency, U(z, ω) is the frequency domain wavefield, and Q(z, ω) is the frequency domain stress.

Equations (4) and (5) are discretized by the fourth-order space finite-difference scheme based on the staggered-grid formulation (Levander 1988). The PML method is used for absorbing boundary conditions along the edges of the model. The discretization format leads to solving the linear system

$$\begin{equation} {\bf AU = F}, \end{equation} \tag{ 6 }$$

where A is the impedance matrix that depends on the angular frequency and medium properties, F is the source term, and U represents wavefields. Equation (6) is solved by using the LU factorization method.

2.2.1. Inverse algorithm

FWI is an optimization problem that can be recast as a linearized least-squares problem that attempts to minimize the misfit between the recorded and modelled wavefields (Tarantola 1984). Only the main features of frequency-domain FWI used in our study are reviewed. A detailed description of the algorithm can be found in Ravaut et al (2004) and Sourbier et al (2009).

The associated objective function to be minimized is defined by

$$\begin{equation} C( {\bf m} ) = \case{1}{2}\delta {\bf U}^t \delta {\bf U}^{\rm *} , \end{equation} \tag{ 7 }$$

where C(m) is the misfit function; δU is the difference between the modelled data U_cal and observed data U_obs; t is the transposition; * is the complex conjugate; and ${\bf m} = [{\bf \rho },{\bf \upsilon },{\bf \eta }]$ refers to parameters where ρ = [ρ₁, ρ₂, ..., ρ_N], ${\bf \upsilon } = [\upsilon _1 ,\upsilon _2 , \ldots ,\upsilon _N ]$, and η = [η₁, η₂, ..., η_N].

We use the steepest-descent method to solve this optimization problem. The relationship between the model perturbation and data residuals is expressed as

$$\begin{equation} \Delta {\bf m} = - \nabla _m C = - {\mathop{\rm Re}\nolimits} \{ {{\bf J}^t \delta {\bf U}^* } \}, \end{equation} \tag{ 8 }$$

where Δm is the model perturbation, ∇_mC is the gradient of the misfit function with respect to m, J is the Jacobian matrix, δU* is the conjugate of the data residual, and Re is the real part.

To derive the expression of the Jacobian matrix, we calculate the partial derivative wavefield by the adjoint method (Pratt et al 1998). After injecting the expression of the partial derivative wavefields into (8) and applying some scaling and regularization to the gradient, we can obtain the reliable perturbation models

$$\begin{equation} \Delta m_i = - ( {{\rm diag}\,{\bf H}_a + \varepsilon {\bf I}} )^{ - 1} {\bf G}_m {\mathop{\rm Re}\nolimits} \left\{ {U^t \left( {\frac{{\partial {\bf A}}}{{\partial m_i }}} \right)^t {\bf A}^{ - 1} \delta {\bf U}^* } \right\}, \end{equation} \tag{ 9 }$$

where ${\rm diag}{\bf H}_a = {\mathop{\rm Re}\nolimits} \{ {{\rm diag}({\bf J}^t {\bf J}^* )} \}$ denotes the diagonal elements of the approximate Hessian H_a and G_m is the smoothing operator.

Finally, the model is updated with the perturbation model

$$\begin{equation} {\bf m}^{( {k + 1} )} = {\bf m}^{( k )} + \alpha ^{( k )} \Delta {\bf m}, \end{equation} \tag{ 10 }$$

where α^(k) is the step length, k is the iteration number, and m^(k) and m^{(k + 1)} are the kth and (k+1)th iterative model.

2.2.2. Calculation of the step length

The parabolic fitting method (Vigh et al 2009) is applied to calculate the optimal step length, which means that we need to find three step lengths. The misfit values corresponding to the three step lengths are referred to as fcost0, fcost1, and fcost2, which should satisfy fcost0 > fcost1 < fcost2. The first step length is set to zero, and the second and the third step lengths are calculated by using the following formulas:

$$\begin{equation} \alpha _{\rho i}^{( k )} {\rm = }\beta _i \cdot \frac{{\rho _{\max }^{( k )} }}{{\Delta \rho _{\max }^{( k )} }},\alpha _{\upsilon i}^{( k )} {\rm = }\beta _i \cdot \frac{{\upsilon _{\max }^{( k )} }}{{\Delta \upsilon _{\max }^{( k )} }},\alpha _{\eta i}^{( k )} {\rm = }\beta _i \cdot \frac{{\eta _{\max }^{( k )} }}{{\Delta \eta _{\max }^{( k )} }}, \end{equation} \tag{ 11 }$$

$$\begin{equation} \alpha _i ^{( k )} = \max \big(\alpha _{\rho i}^{( k )} ,\alpha _{\upsilon i}^{( k )} ,\alpha _{\eta i}^{( k )} \big),\quad i = 1,2, \end{equation} \tag{ 12 }$$

where $\Delta \rho _{\max }^{( k )}$, $\Delta \upsilon _{\max }^{( k )}$, and $\Delta \eta _{\max }^{( k)}$ are the maximum values of $\Delta \rho _{}^{( k )}$, $\Delta \upsilon _{}^{( k )}$, and $\Delta \eta _{}^{( k )}$, respectively. The coordinates of $\upsilon _{\max }^{( k )}$ and $\eta _{\max }^{( k )}$ are same as the ones of $\Delta \rho _{\max }^{( k )}$, $\Delta \upsilon _{\max }^{( k )}$, and $\Delta \eta _{\max }^{( k )}$ in the kth iteration. β_i is a coefficient where β₁ is for the second step length and β₂ is for the third step length. $\alpha _i ^{( k )}$ is the maximum of $\alpha _{\rho i}^{( k )}$, $\alpha _{\upsilon i}^{( k )}$, and $\alpha _{\eta i}^{( k )}$; and $\alpha _1 ^{( k )} ,\alpha _2 ^{( k )}$ are the second and third step lengths, respectively.

When Fourier transformation is applied to (1), we can obtain

$$\begin{eqnarray} &&\fl\omega ^2 U( {z,\omega } ) + \frac{\partial }{{\partial z}}\left[ {( {\rho ( z )\upsilon( z )^2 + {\mathop{\rm i}\nolimits} \omega \eta( z )} )\frac{{\partial U( {z,\omega } )}}{{\partial z}}} \right]\nonumber\\ &&= F( {z,\omega } ). \end{eqnarray} \tag{ 13 }$$

From (13), we can find that the derivatives of impedance matrix to ρ, υ, and η can be approximated by υ(z)², 2ρ(z)υ(z), and iω. The magnitude of the derivatives with respect to ρ and υ are the same; and the derivative with respect to η is dominated by the angular frequency. Given that FWI works well in low frequencies, the magnitude of the derivative to η is less than 3. For example, for a given angular frequency (such as 35) and grid point (such as the fourth one), the magnitudes of the partial derivative wavefield for ρ, υ, and η are –12, –13, and –18, respectively. From this point of view, ρ and υ are more sensitive than η to the misfit function.

From the term $\rho ( z )\upsilon ^2 (z){\rm + }{\mathop{\rm i}\nolimits} \omega \eta ( z )$, we can see that η is related to the angular frequency. For the given $\rho ( z )\upsilon ^2 (z)$ and angular frequency, the absorbing term increases with increased η. When η increases to some extent, η is dominant in the misfit function. Then, absorption becomes so intensive that no wave may be recorded, which results in difficult reconstruction of the models. To obtain reliable reconstructed results, an appropriate η model should be selected. Our tests (not shown) reveal that when the average magnitude of $\rho ( z )\upsilon ^2 (z)$ is close to that of ωη(z), the three parameters can be well rebuilt.

Parameters ρ and υ are set to satisfy Gardner's formula (Ma and Jie 2005) during inversion, so we only analyse the variation in the misfit functions as a function of υ and η. Figure 1(a) shows a seven-layer model geometry, and figures 1(b)–(d) are the true models for density, velocity, and viscosity coefficient. Table 1 displays the model values for each layer. Shot is set at a depth of 200 m and 120 receivers are located at an interval of 10 m in the borehole.

Zoom In Zoom Out Reset image size

We test the variation in the misfit function on the seven-layer model by changing only the values of the sixth layer. Figure 2(a) is the misfit function curve of the relative variation to the true υ, and figure 2(b) is the misfit function curve of the relative variation to the true η. The relative variation is obtained by

$$\begin{equation} {\rm Mvar = }\frac{{m - m_{{\rm true}} }}{{m_{{\rm true}} }} \times 100\% , \end{equation} \tag{ 14 }$$

where m is the changing model value and m_true is the true model value.

Zoom In Zoom Out Reset image size

Figure 2(a) shows that the misfit function increases with increased difference between the true and changed velocities, which means that velocity has a high sensitivity in the data, as shown by the well-convex shape. For parameter η, the misfit function slowly varies when η is smaller than the true model, but significantly varies when η is higher than the true ones. For the same relative variation, the magnitude of the misfit function to υ is larger than that to η by a magnitude of 2–3, such as the relative variation 50% with a magnitude of 2. Thus, υ is much more sensitive than η to the misfit function.

In this section, we compare two inversion strategies, and the better one is used in subsequent numerical tests. As mentioned above, υ is more sensitive so we first try a one-by-one parameter inversion strategy. In the first step, we reconstruct only ρ and υ while fixing η as the initial model. In the second step, only η is rebuilt while ρ and υ are fixed as the reconstructed models in the first step. The seven-layer model shown in figures 1(b)–(d) is tested by using this inversion strategy. The initial models are homogeneous models with 2400 kg m⁻³, 3600 m s⁻¹, and 3.4 × 10⁷ Pa s for density, velocity, and viscosity, respectively. Nine frequencies between 3.5 and 20 Hz are used, and the iteration number is 1000 for each frequency. Figure 3 shows the reconstructed results where ρ and υ well match the true models while η is a poor match. The reason is that the reconstructed ρ and υ do not well match the true models, and the errors caused by ρ and υ change into noise to η, which results in the poor convergence of η. The inverted models in figure 3 show that ρ and υ can be well reconstructed with η being the initial model, whereas η cannot be rebuilt with ρ and υ having some errors to the true models. The result reveals that η is more sensitive to noise than ρ and υ.

Zoom In Zoom Out Reset image size

Then, we try the second inversion strategy, namely, the joint inversion strategy. We reconstruct the three parameters jointly and obtain the step lengths by the parabolic fitting with (11)–(12). In the early inversion stage, υ is more sensitive than η, so ρ and υ are dominant and the step length is more suitable for them. As iteration is continued, ρ and υ become closer to the true models, and then η dominates and the step length becomes more appropriate forη. Consequently, ρ and υ can be well reconstructed in the early stage, whereas η can be well rebuilt in the later stage of inversion. We also test the joint inversion strategy on the seven-layer model in figures 1(b)–(d), and the parameters used in the inversion are the same as the ones in figure 3. Figure 4 shows that the reconstructed models are in agreement with the true models. In the following numerical tests, we adopt the joint inversion strategy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To obtain global minimum solutions and improve the efficiency of waveform inversion, we compare three frequency-selection strategies. The first strategy, referred to as the single-frequency method (Brossier et al 2009, Kim et al 2011), is implemented one frequency by one frequency and from low to high frequency. The second strategy, referred to as the Bunks approach (Brossier et al 2009, Kim et al 2011), is an adaptation in the frequency domain of the time-domain approach of Bunks et al (1995). It consists of successive inversions of overlapping frequency groups. The third strategy, referred to as the group inversion approach (Brossier et al 2009), consists of successive inversions of slightly overlapping frequency groups.

In our study, we test and compare the three frequency selection strategies on the seven-layer model in figures 1(b)–(d). The initial models are homogeneous models with 2400 kg m⁻³, 3600 m s⁻¹, and 3.4 × 10⁷ Pa s for density, velocity, and viscosity coefficient, respectively. The inverted frequencies for each frequency selection strategy are empirically selected as shown in table 2. The main frequency of the Ricker wavelet is 15 Hz, and the space interval is 10 m; 1 shot and 120 receivers are used for the zero-offset VSP test. The iteration number is 1000 for each frequency or frequency group. We compare the inversion times and accuracies of the final results of each strategy. The accuracy of the reconstructed model is assessed by the normalized model error (NME):

$$\begin{equation} {\rm NME} = \sqrt {\frac{1}{N}\left\| {\frac{{M_{{\rm cal}} - M_{{\rm true}} }}{{M_{{\rm true}} }}} \right\|_2^2 } , \end{equation} \tag{ 15 }$$

where M represents the parameters ρ, υ, and η; M_cal and M_true are the inverted and true models, respectively; and N is the total number of each parameter.

Figures 4–6 show the reconstructed models for the single-frequency approach, Bunks approach, and group inversion approach, respectively. Table 2 shows the accuracies of the final models and calculation time for three frequency-selection strategies. Given the sufficiently large iterations, the reconstructed models of ρ, υ, and η well match the true models. Table 2 shows that the highest accuracies of the final results are 8.82 × 10⁻⁴, 0.0035, and 0.0187 for ρ, υ, and η by using the Bunks approach. However, this approach is the most time consuming one at 197.5 minutes, which is about four times the inversion time for the single-frequency method (46.07 minutes) and about three times the calculation time for the frequency group approach (67.73 minutes). The accuracies of the final models by the single-frequency method (i.e. 9.88 × 10⁻⁴, 0.004, and 0.0208 for ρ, υ, and η, respectively) are close to the accuracies by the Bunks approach but less time consuming. The accuracies (0.0014, 0.0056, and 0.0289 for ρ, υ, and η, respectively) and inversion time of the group inversion approach lies between the Bunks approach and the single-frequency approach. Considering the accuracies and calculation time without parallel implementation, we believe the single-frequency approach is preferred.

Zoom In Zoom Out Reset image size

To demonstrate further that the inverted models match the true models, we present the responses to the true, initial, and final models in figure 3. The time forward modelling is computed by finite difference in time domain, where the discrete pattern and PML boundary are the same as frequency domain modelling. The seismographs are recorded in 0.7 s, as shown in figure 7. Figure 7(a) is the response to the true models without considering the attenuation, whereas figure 7(b) is for the true model with attenuation. Clearly, the wave energy with attenuation is weaker than that without attenuation. Figures 7(b), (d) and (f) show that the seismograph for the final models is close to the one for the true models, which means that the reconstructed models are reliable and the inversion strategy used is feasible.

Zoom In Zoom Out Reset image size

To assess the influence of the initial models on the final models, we compare two initial models. One is the 70-point smoothing of the true models in figures 1(b)–(d) and the other is the homogeneous models with 2600 kg m⁻³, 5000 m s⁻¹, and 5.0 × 10⁷ Pa s for density, velocity, and viscosity coefficient, respectively. The single-frequency approach is used, and the inverted frequencies are the same as the ones in table 2 for this approach. The iteration number is 500 for each frequency. The true model is shown in figures 1(b)–(d). Figures 8 and 9 are the final models for the 70-point smoothing initial models and the homogeneous initial ones. The accuracies of ρ, υ, and η are 0.0014, 0.0055, and 0.0236 for the 70-point smoothing initial models. For the homogeneous initial models, the accuracies are 0.0014, 0.0055, and 0.0226, respectively. We can observe that the differences between the final models for the 70-point smoothing initial models and homogeneous initial ones are so small that they can be ignored. We can conclude that the influences of the initial models on the final results are small.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To determine the influences of inverted frequencies and the total number of iterations on the final results, we compare the final results with different frequencies and iterations. In this comparison, the true models of ρ and υ are the same as the ones in figures 1(b) and (c). However, the true model of η differs from the one in figure 1(d), which is modified from the VSP example of SEISCOPE (Operto et al 2007). The initial models are homogeneous models with 2400 kg m⁻³, 3600 m s⁻¹, and 2.7 × 10⁷ Pa s for ρ, υ, and η, respectively.

We first compare the results of 500 iterations (figure 10) with 1000 iterations (figure 11), where four inverted frequencies (1.57, 3.53, 5.49, and 7.84 Hz) are used. The accuracies of ρ, υ, and η with 500 iterations are 0.0015, 0.0061, and 0.0204, respectively, whereas the accuracies with 1000 iterations are 0.0046, 0.0012, and 0.0153, respectively. The differences in the inverted ρ and υ between the two iterations are small, whereas the difference in η is large, especially in the position of layer stripping. We can observe that larger iterations result in better reconstruction of η, which means that η is more dependent on the iterations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We then compare the final models from four inverted frequencies (figure 10) with the ones obtained by eight inverted frequencies (figure 12), which are the same as the ones in table 2 for the single-frequency strategy. In both cases, the iteration number is 500. The accuracies of ρ, υ, and η with eight frequencies are 0.0013, 0.0051, and 0.0174, respectively. Table 3 compares the accuracies of the two kinds of frequencies, and the results with eight frequencies and 500 iterations are much better. We then compare the final models by 1000 iterations and four inverted frequencies with the ones by 500 iterations and eight inverted frequencies, and results show that they have close inversion accuracies.

Zoom In Zoom Out Reset image size

From table 3 and figures 10–12, we can conclude that the total number of iterations, i.e. the number of inverted frequencies multiplied by the number of iterations, plays an important role in obtaining high-accuracy final models. When fewer inverted frequencies exist, we can improve the final models by increasing the number of iterations, and for more frequencies the final models can be well rebuilt with less iteration.