Elastic versus acoustic inversion for marine surveys

SUMMARY

Full wavefield inversion (FWI) is a powerful and elegant approach for seismic imaging that is on the way to becoming the method of choice when processing exploration or global seismic data. In the case of processing marine survey data, one may be tempted to assume that acoustic FWI is sufficient given that only pressure waves exist in the water layer. In this paper, we pose the question as to whether or not in theory—at least for a hard waterbottom case—it should be possible to resolve the shear modulus or S-wave velocity in a marine setting using large offset data. We, therefore, conduct numerical experiments with idealized marine data calculated with the elastic wave equation. We study two cases, FWI of data due to a diffractor model, and FWI of data due to a fault model. We find that at least in idealized situation, elastic FWI of hard waterbottom data is capable of resolving between the two Lamé parameters λ and μ. Another numerical experiment with a soft waterbottom layer gives the same result. In contrast, acoustic FWI of the synthetic elastic data results in a single image of the first Lamé parameter λ which contains severe artefacts for diffraction data and notable artefacts for layer reflection data. Based on these results, it would appear that at least the inversions of large offset marine data should be fully elastic rather than acoustic, unless it has been demonstrated that for the specific case in question (offsets, model and water depth, practical issues such as soft sediment attenuation of shear waves or computational time), an acoustic-only inversion provides a reasonably good quality of image comparable to that of an elastic inversion. Further research with real data is required to determine the degree to which practical issues such as shear wave attenuation in soft sediments may affect this result.

Waveform inversion, Inverse theory, Image processing, Computational seismology, Controlled source seismology

1 INTRODUCTION

Full wavefield inversion (FWI) in the acoustic approximation was first proposed in the early 1980s (Lailly 1983; Tarantola 1984). During the late 1980s, FWI for the fully elastic case was developed and demonstrated (Mora 1987b; Tarantola 1988) but has not been commonly applied to exploration data due to significantly higher computational cost.

Much of the research is done in the acoustic approximation to decrease computational cost, or based on the assumption that most of the energy is in the compressional waves, particularly in the case of marine data where the waves must convert to pressure waves in the water layer. However, if one does decide to perform the inversion in the acoustic approximation, one must pose the question as to what are the consequences of this approximation. Clearly, the earth is elastic to first order, so both compressional and shear waves exist in the solid layers. We may wonder whether the acoustic approximation is appropriate. What is the consequence on the quality of the image of using the acoustic approximation, at least under ideal circumstances? And what extra information can be gained by using the elastic wave equation?

A number of papers has been published that aimed to study the consequences of applying acoustic wavefield inversion to elastic data on the P-wave velocity result (Barnes & Charara 2009; Vigh et al. 2009; Marelli et al. 2012; He & Plessix 2016). These papers showed inversion examples of 2-D and 3-D models and noted that artefacts exist on the P-wave velocity inversion result when acoustic inversion is carried out using elastic data. Several authors have proposed methods to help mitigate the artefacts on the P-wave image that resulted from applying acoustic inversion to elastic data (Agudo et al. 2016; Albertin et al. 2016; Willemsen et al. 2016). Work towards the goal of mitigating artefacts of using the acoustic approximation includes research that improves amplitude versus incidence angle and use of multiple propagations of acoustic waves by adjusting the density (Mulder & Plessix 2008; Chapman et al. 2014; Hobro et al. 2014). Such research is focussed on the artefacts resulting from incorrect amplitude versus incidence angle of the P–P reflections, rather than artefacts due to the acoustic inversion attempting to fit mode converted arrivals. Our goal is not to mitigate the effects of artefacts. Rather, we aim to study the nature of all artefacts, particularly those related to fitting mode converted energy, and to demonstrate what extra information can be gained by fully elastic inversion, for the case of large offset (incidence angle) data when mode conversions become significant. We do not study the small offset case where mode conversions can be neglected to first order and the acoustic approximation is expected to be adequate. Based on the work of Shipp & Singh (2002) in which 2-D elastic FWI was applied to wide aperture marine seismic data to invert for the P-wave velocity, mode conversions play an important role in improving the P-wave velocity image, even for marine data.

The goals of this paper are: (1) to conduct numerical experiments for ‘hard waterbottom’ and ‘soft waterbottom’ models to enable the quality of full elastic wavefield inversion results to be compared to full acoustic wavefield inversion results of elastic data, and (2) to determine whether elastic inversion can resolve a second image (e.g. the S-wave velocity or shear modulus) compared to acoustic inversion which is only capable of resolving a P-wave velocity image (or P-wave modulus). We limit this study to the case of a simulated marine survey and hence, the uppermost layer is acoustic and can support only pressure waves. This case was chosen because some researchers have argued that the lack of shear waves in the water layer justifies the application of the acoustic approach.

The following numerical experiments involve FWI in the elastic case and in the acoustic approximation. We apply the classical time domain FWI formulation described in Tarantola (1984, 1988) and Mora (1987b) and use a relatively broad-band source wavelet (a first derivative of a Gaussian). As such, the inversion results presented are FWIs in which both the amplitudes and phases of the recorded wavefields must be matched. We assume that a smooth background velocity field is known, and that this is provided as a starting model so we will not be considering issues of cycle-skipping that can occur if the smooth initial model is not close enough to the true model. The numerical experiments aim to investigate the assumption that acoustic inversion provides adequate results in the case of large offset data, in which offsets are comparable to or larger than the depths of interest. Namely, data in which mode conversions will be present. We do not try to mitigate the effects of using acoustic inversion. We rather focus on inversions using a simple diffractors model, and a simple anticline and thrust fault model to study acoustic inversion relative to elastic inversion. Hence, our goal is to clarify, for simple and easy to understand models, the kind of artefacts that may occur when applying FWI in the acoustic approximation, and whether elastic FWI is able to resolve added information such as on the shear wave modulus.

2 FORWARD MODELLING AND INVERSION METHODOLOGY

In the following, we utilize a staggered finite-difference grid formulation to solve the elastic wave equation as described by Igel et al. (1995). This method is very briefly reviewed in Appendix  A for completeness.

The FWI method used in this paper is the time domain formulation as described by Tarantola (1988) and Mora (1987b). Again, the basic equations and the conjugate gradient algorithm used in this paper are briefly reviewed in Appendix  B for completeness. The inversion parameters in this paper are λ, μ and ρ. A careful choice of parameters (He et al. 2018) that considers radiation patterns can potentially do an even better job. However, such considerations are outside the scope of this paper.

2.1 Conversion of model units to physical units

In this paper, we have chosen dimensionless units (i.e. scaled to Δx = Δz = 1). To see what size of model the water layer model would correspond to, we note that simulations were done at about eight gridpoints per fundamental wavelength in the water layer. Hence, the conversion to physical units can be achieved as follows. For a model that is 500 × 200 units, this equates to a model that is 500/8 × 200/8 = 62.5 × 25 wavelengths. For V_p(water) = 1500 m s⁻¹ and a source at 10–30 Hz with a fundamental frequency of 20 Hz, we have a fundamental wavelength of 1500/20 = 75 m, so the model size is 75 × 62.5 m in x, and 25 × 75 m in z, namely, 4687 × 1875 m². Hence, the 500 × 200 model is roughly 5 km × 2 km in size at 20 Hz. All the examples in this paper can be similarly scaled. Specifically, at 20 Hz, the grid spacing is approximately Δx = Δz = 75/8 ∼ 10 m, and hence, a model size of nx × nz is nx × 10 m by nz × 10 m in size.

The calculations were done with a value of Δt* = 0.175 where here, the superscript* implies model units. We can convert from the model time units to time as follows. In the water layer, we have λ* = 1 and ρ* = 1 so we have a model velocity of |$V^*_p = \sqrt{\lambda ^* / \rho ^*} = 1$|⁠. So for the 20 Hz example above where we have Δx = 10 m and V_p(water) = 1500 m s⁻¹, we can write

\begin{equation*} V_p {\Delta t \over \Delta x} \ = \ 1500 {\Delta t \over 10} \ = \ 150 \Delta t \ = \ V_p^* {\Delta t^* \over \Delta _x^*} \ = \ 0.175 , \end{equation*}

(1)

and hence, the time step in seconds, Δt, is given by

\begin{equation*} \Delta t \ = \ 0.175/150 \ = \ 0.0011667 \ \mathrm{ s}. \end{equation*}

(2)

Hence, in our example for modelling a 20 Hz fundamental frequency, we have a grid spacing of Δx = 10 m and a time step Δt = 0.0011667s. In this case, a time axis that is 1500 points long is equivalent to a time axis of 1500×0.0011667 ≈ 1.75 s. Therefore, in all the following plots that are labelled in time steps and grid spacing, a simple conversion of model units to physical units, for a 20 Hz fundamental frequency, is to multiply the x and z ordinates by 10 m, and to multiply the time axis by 0.0011667 s. Thus, time axes that are, respectively, 1200 and 1500 model time units long are equivalent to time axes that are 1.4 and 1.75 s long.

2.2 Simulation of a water layer model

Fig. 1 shows the model consisting of a water layer above two solid layers. And Fig. 2 shows snapshots of a simulation of a pressure source located centrally at the top of the model. The simulation was made using the staggered grid finite-difference scheme of Igel et al. (1995) to solve eqs (A5)–(A7). These snapshots show only pressure waves in the uppermost acoustic water layer, and significant mode converted energy in the solid layers below the waterbottom. For example, one observes the P−P and P−S mode conversion when the pressure wave enters the solid layer at a depth of 60 units (Fig. 2b), and one observes the P−P, P−S, S−P and S−S reflected and transmitted waves when the waves intersect the halfspace boundary at a depth of 120 units (Fig. 2c). This illustrates the ability of the staggered finite-difference scheme to model a water layer above solid layers.

Figure 1.

A 2-D layer model. The uppermost layer represents water and has properties λ_w = 1, μ_w = 0 and ρ_w = 1. The middle layer has properties λ = μ = ρ = 1. The lower half-space has properties λ = 0.1, μ = 0.5 and ρ = 1. Hence, the P-wave velocities of the different layers, given by |$v_\mathrm{ p} = \sqrt{\lambda + 2 \mu /\rho }$|⁠, are, repectively, 1 for the water layer, |$\sqrt{3}$| for the solid layer below the waterbottom and |$\sqrt{1.1} \approx 1.05$| for the lower half-space. Note that the S-wave velocity in the middle solid layer is equal to unity (i.e. it has the same speed as the P-wave in water).

Figure 2.

Snapshots of a simulation of a model consisting of a water layer with a depth of 60 above a layer over a half-space where the half-space interface is at a depth of 120. Mode conversions such as P−P, P−S, S−P and S−S are visible in the solid layers below z = 60 whereas only P-waves exist in the water layer above z = 60.

3 DIFFRACTORS MODEL

In the following simulations, we define a model that is 256 × 174 gridpoints in size, with the lower part of the model below 128 gridpoints deep being an absorbing band. Hence, the domain of the model spans from x_left = 0 through to x_right = 255Δx = 255. We define three sets of four diffractors in this model for ρ, λ and μ located, respectively, at x_ρ = x_left + 0.3(x_right − x_left) = 76.5Δx, x_λ = x_left + 0.5(x_right − x_left) = 128Δx and at x_μ = x_left + 0.7(x_right − x_left) = 178.5Δx as shown in Fig. 3, where each set of diffractors is located at four depths (40, 60, 80 and 100 gridpoints deep). The model has a narrow water layer that is 10 grid spacings deep, and the model is 128 gridpoints deep before an absorbing layer is specified. In contrast, the horizontal extent of the model is 256 gridpoints wide and an absorbing band that is 30 gridpoints wide is specified on either side of this model. As such, for a source that is located centrally, the maximum offset is 128 gridpoints wide less than the absorbing layer (30 gridpoints) = 98 gridpoints. Hence, we have offset to depth ratios spanning (128 − 30)/40 = 2.45 through to (128 − 30)/100 ≈ 1. These can be considered as moderate to large offsets relative to depth, a realistic assumption for typical seismic surveys.

Figure 3.

2-D diffractors model: first Lamé parameter λ (a), second Lamé parameter μ (b) and density ρ (c).

3.1 Elastic inversion of diffractors—single shot

In this example, a single shot record with the pressure source located centrally at the surface of the model is inverted. The residuals after 100 iterations were small and hence, the inversion was stopped after 100 iterations. The inversion results after 100 iterations are shown in Fig. 4. These indicate that the two Lamé parameters λ and μ are well resolved from one another for this example which involves moderate to large offsets relative to the depths of the diffractors (offset to depth ratios of around 2.5:1 for the shallow diffractors, and 1:1 for the deep diffractors). Namely, one can see clearly the four diffractors in λ at x_λ in Fig. 4 and four diffractors in μ at x_μ with only minor artefacts elsewhere. Similarly, the inversion result for density ρ has four strong perturbations at the correct locations of the density diffractors x_ρ, and also smaller ghost perturbations at the locations of the diffractors in λ and μ. This indicates an imperfect resolvability between density and the Lamé parameter, and a relatively good resolvability between the two Lamé parameters λ and μ.

Figure 4.

The diffractors model single shot elastic inversion result of elastic marine data after 100 iterations: first Lamé parameter result λ₁₀₀ (a), second Lamé parameter result μ₁₀₀ (b) and density result ρ₁₀₀ (c).

Fig. 5 shows the error versus iteration for this example and indicates that after 100 iterations, almost all of the energy is modelled (i.e. the error approaches zero). The initial residuals and the residuals after 100 iterations are shown at the same scale in Fig. 6. After 100 iterations, there is no visible energy left in the residuals again indicating that the simulated data through the inversion result nearly perfectly matched the simulated data through the diffractors model.

Figure 5.

The error versus iteration for the elastic inversion of elastic marine data for the diffractors model.

Figure 6.

The diffractors model initial residuals (a) and the residuals after 100 iterations for the elastic inversion of elastic marine data.

3.2 Acoustic inversion of diffractors—single shot

In this section, we investigate application of an acoustic inversion of the elastic wavefield data, which is the same as previous section. We found that when we applied an acoustic inversion to the fully elastic data, most of the effort was spent trying to futilely correct the seafloor bottom reflection. As such, we stripped the seafloor bottom reflection from the simulated data prior to doing the acoustic inversion. This was achieved by doing an elastic and acoustic simulation in a model with the water layer only, and subtracting this waterbottom reflection data from the simulated elastic dataset, and from the acoustic simulations carried out during the inversion. In this way, the acoustic inversion only tries to match the data due to diffractions and reflections from within solid layers below the seafloor, thereby providing a good chance for an acoustic inversion to succeed when inverting for the subsurface below the waterbottom. In contrast, the waterbottom reflection was not stripped from the elastic inversion example as it was automatically modelled correctly because the starting model had the correct water depth and elastic properties.

The diffractors and results are shown in Figs 7 and 8. The first iteration results show the four diffractors in λ in the correct positions, albeit a lower frequency image than obtained in the elastic inversion. The density image also shows four perturbations located at the four positions of the λ diffractors. Hence, after one iteration, the acoustic inversion was only able to correctly resolve the diffractors in λ, but was unable to resolve the density diffractors. Furthermore, after 100 iterations, the initial result for λ was degraded with large artefacts appearing throughout the model. This is due to the inversion trying to match the mode converted energy which is not possible using acoustic simulations. As such, the inversion placed artificial diffractors throughout the model such that the overall error, shown in Fig. 9, was reduced. One can see from the plot of error versus iteration that the error only decreased by an additional 25 per cent from the 5th to the 100th iteration (i.e. from 0.017 to 0.013). In contrast, the residual in the elastic inversion was already lower than the acoustic case after only two iterations, and the residual continued to drop close to 0 after the 100 iterations.

Figure 7.

A 2-D diffractors model of the first Lamé modulus λ of the subsurface (a), and the acoustic inversion result of elastic marine data after 1 iteration (b), and after 100 iterations (c). The water layer has a value of μ_water = 0 and hence, can be seen as a black band at the top of the model in (b). The diffractors are point perturbations in the model at four different depths (z = 40, 60, 80 and 100). The locations of the diffractors are as follows: density diffractors at x = 76.5, λ diffractors at x = 128 and μ diffractors at x = 178.5.

Figure 8.

A 2-D diffractors model of the density ρ of the subsurface (a), and the acoustic inversion result of elastic marine data after 1 iteration (b), and after 100 iterations (c).

Figure 9.

The error denoted E versus iteration for the acoustic inversion of elastic marine data for the diffractors model.

If one considers the residual for the acoustic case, it is clear why the error could not be reduced much after the first few iterations. One can see from the plot of residuals shown in Fig. 10 that the residuals correspond to mode conversions from the μ diffractors. Such mode conversions are centred at x_μ ≈ 180 and have a lower velocity (i.e. rate of moveout per unit time at large offset → S-velocity) than the diffractions centred at x_λ. However, the acoustic inversion does not model such mode conversions, and therefore, it cannot resolve this different rate of moveout with time as a diffractor. Rather, it attempts to create λ and ρ perturbations that emulate the effect of the mode-converted diffractions.

Figure 10.

The diffractors model initial residuals (a) and the residuals after 100 iterations for the acoustic inversion of elastic marine data.

Note that a test was also conducted of a supershot gather consisting of 18 different shot locations. The result was similar to the single shot result, and for this reason, has not been included here. And worse, the more iterations that are made, the worse the result becomes due to the inversion adding artificial perturbations in the model to try to emulate the effect of the μ diffractors.

Note that the absence of any perturbations in the inversion images above 10 units deep is due to energy in the gradient being muted out in the water layer which is considered as a known region of the model. In the case of the acoustic inversion of elastic data, to give the best possible chance to the acoustic inversion, we stripped the waterbottom reflection, and hence, we can consider that in our case, the model is known perfectly down to the start of the solid layer below the waterbottom.

4 ELASTIC INVERSION OF ELASTIC MARINE DATA (FAULT MODEL)

A model was defined consisting of a simple Gaussian anticline with two thrust faults as shown in Figs 11(a) to 12(a). The layer properties for μ were independent to those of λ. For example, the low λ layers were layers 3 (dark blue) and 7 (blue between two yellow layers), cf. the low μ layers were layers 4 (bright blue) and 8 (aqua-green below yellow). The layer properties are shown in Table 1. Note that from the layer properties, one can see that the maximum values of the moduli are 1.4 × the amplitude of the lowest values of unity. This corresponds to V_p and V_s ranging up by 20 per cent from the shallowest to deepest layer. This is a somewhat smaller range than normal. However, this small range was necessary due to limited computer resources (a home workstation). A supershot consisting of 28 pressure sources equally spaced along the surface of the model was simulated using the staggered finite-difference scheme. The start times of each shot were randomized to eliminate artificial correlations in the model due to pre-summing of the shots into a supershot. This randomization of source times was re-done after every two iterations of the conjugate gradient algorithm to minimize the correlations between the shots (see also Schiemenz & Igel 2013). An example of a single-shot profile for data simulated in the model is shown in Fig. 13. One can see seafloor multiples and strong P−P reflections at zero offset and extending to large offsets. In addition, one can see the slow direct water P-wave, and following behind the direct P-wave, one can observe weak reflections with a similar velocity to the water velocity, and with zero amplitude at zero offset. These reflections are P−S−S−P modes and are what allow the shear modulus to be resolved. Fig. 14 shows a 28 shot supershot gather for the fault model which represents the dataset used to perform the inversion.

Figure 11.

A 2-D fault model of the first Lamé modulus λ of the subsurface (a), initial model (b) and the elastic inversion result of elastic marine data after 200 iterations (c).

Figure 12.

A 2-D fault model of the second Lamé modulus μ of the subsurface (a), initial model (b) and the elastic inversion result of elastic marine data after 200 iterations (c). The value of μ in the upper water layer is μ_water = 0 (see Table 1).

Figure 13.

A single shot gather due to a source located in the centre of the fault model shown in Figs 11(a) and 12(a) simulated using absorbing boundaries. The supershot gather consisted of 28 such shots.

Figure 14.

The 28 shot supershot gather for the fault model shown in Figs 11(a) and 12(a). The supershot was calculated as a single simulation with 28 sources located at a regular increment along the top of the model and triggered at random start times. The boundary conditions were circular in the x-direction to efficiently simulate the equivalent case of summing 28 large offset gathers without any edge effect.

Figure 15.

The elastic inversion result for density ρ of the elastic marine data after 200 iterations. The true and starting models for density were both constant with ρ = 1 and hence, are not shown here.

Figure 16.

The error versus iteration for the elastic inversion of elastic marine data.

Table 1.

Lamé parameters and density for the elastic 2-D fault model and its equivalent acoustic 2-D fault model. Values of the elastic model Lamé parameters λ = λ_e and μ = μ_e (first two columns), and the acoustic model Lamé parameters λ_a, μ_a (last two columns). The density ρ (middle column) is the same for both the elastic and acoustic models. The column denoted M = λ + 2μ represents the modulus for pressure waves, which is what can potentially be resolved by an acoustic inversion. Note that the values of λ and μ differ so these models are independent of one another. For example, layers 3 and 7 have a low value of λ sandwiched between higher λ layers, whereas layers 4 and 8 have a low value of μ sandwiched between higher μ layers. Note that the P-wave modulus denoted M = λ + 2μ has a different character than λ and μ, with layers 3 and 4 having somewhat lower M than the surrounding layers.

Layer	λe	μe	ρe = ρa	M = λa = λe + 2μe	μa
water	1	0	1	1	0
1	1	1	1	3	0
2	1.06	1.07	1	3.2	0
3	0.95	1	1	2.95	0
4	1.04	0.95	1	2.94	0
5	1.08	1.04	1	3.16	0
6	1.22	1.08	1	3.38	0
7	1.05	1.22	1	3.49	0
8	1.25	1.05	1	3.35	0
9	1.4	1.4	1	4.2	0
10	1.32	1.32	1	3.96	0

Our goal here is to compare elastic inversion of elastic data with acoustic inversion of elastic data under ideal circumstances for both. As such, periodic boundary conditions were used in the x-direction and 1500 time steps were calculated of the finite-difference simulation of elastic waves on a 502×232 sized finite-difference mesh. This is enough to record seismic waves to large offsets with ratios of offset to depth in the range 1:1 to 3:1. Also, as the aim is to identify artefacts due to using an acoustic inversion of elastic data, and has nothing to do with quantification of acoustic or elastic inversion results, no colour scales are provided in the following sections as these would distract the reader from this aim. Furthermore, due to the large contrast in elastic parameters in the water layer versus the rest of the model, a standard rainbow colour scheme would mainly visualize the water layer, with almost no dynamic range left to visualize the elastic properties in the solid layers, or the soft waterbottom layer. In the following, dark blue is used to visualize the water layer and bright blue is used to visualize the soft sediments layer, whereas a standard colour rainbow is used to visualize the elastic parameters in the solid layers. As is typically done for idealized test case inversions, we set the initial models to be smoothed versions of the exact models as shown in Figs 11(b) and 12(b).

Figs 11(c) to 15 show the results of an inversion after 200 iterations of a conjugate gradient algorithm (Note: the residuals after 200 iterations were small and hence, the inversion was stopped at 200 iterations). These results show that λ and μ are well resolved from one another, but the density ρ, which should be constant, has a ghost image of the Lamé parameter (see Fig. 15). This is due to the imperfect resolvability of ρ from λ and μ as was also seen in the inversion of the diffractors.

Figure 17.

The initial residuals (a) and the final residuals after 200 iterations (b) for the elastic inversion of elastic marine data.

Figure 18.

A 2-D fault model of the P-wave modulus = first Lamé modulus M = λ_a of the subsurface (a), initial model (b) and the acoustic inversion result of elastic marine data after 200 iterations (c). This result has a lower clarity than the result of elastic inversion of elastic data shown in Fig. 11, with a cross-hatched appearance that are mode-conversion artefacts.

Figure 19.

The acoustic inversion result for density ρ of the elastic marine data after 200 iterations. The true and starting models for density were both constant with ρ = 1 and hence, are not shown here.

Fig. 16 shows the error as a function of iterations and demonstrates that most of the energy has been matched after the 200 iterations. Note that the zig-zag pattern of error up and down is due to the source randomization procedure after every second iteration. While this randomization procedure seems to slow the convergence relative to a case where source times are left constant, tests have shown that this randomization leads to clearer images. Namely, the time randomization of shots when summing the shots into supershots yields images without artefacts due to source correlations (i.e. from source cross-talk) which occur if one sums the shots into a supershot with a fixed timing. Hence, rather than slowing down an inversion, the pre-summing of shots into a supershot speeds up the inversion by a factor approaching the number of shots when compared to an inversion of individual shots.

Fig. 17 shows the initial residual, the initial supershot gather with the seafloor bottom reflection stripped, compared to the residual after 200 iterations. This plot shows that the final residual is close to zero and hence, the inversion has matched most of the data in the supershot record.

5 ACOUSTIC INVERSION OF ELASTIC MARINE DATA (FAULT MODEL)

In this section, we wish to apply an acoustic inversion to the data calculated in the elastic model. That is, using the same elastic data as was calculated in the previous section, we perform an acoustic inversion. In this case, to ensure the physics of the inversion is comparable to the elastic case, we have initialized the acoustic model to the P-wave modulus denoted byM = λ + 2μ, and performed finite-difference calculations using this P-wave modulus in place of λ, and with a shear modulus μ = 0. Namely, we define λ_a and μ_a as the acoustic Lamé parameters. These are related to the elastic Lamé parameters λ = λ_e and μ = μ_e as follows:

\begin{equation*} \begin{array}{ll}\lambda _\mathrm{ a} \ & = \ M \ = \ \lambda _\mathrm{ e} + 2 \mu _\mathrm{ e} \\\mu _\mathrm{ a} \ & = 0 . \end{array} \end{equation*}

(3)

Hence, P-waves will propagate with the same speeds in the acoustic and elastic models, but no shear waves will propagate in the acoustic model.

The elastic model used to calculate the data is shown in Figs 11(a) to 15. However, what is potentially resolvable by the acoustic inversion is the P-wave modulus M. As such, the true acoustic model that is potentially resolvable by an acoustic inversion is shown in Figs 18(a) to 19. The scales are different for λ_a = M compared to λ_e = λ as per Table 1 (i.e. λ_a = M ≈ 3λ_e). Again, the initial model is set to be a smoothed version of what can potentially be resolved, which in this case, is the exact model for M = λ_a as shown in Fig. 18(b).

Figs 18(c) and 19 show the results of an acoustic inversion after 200 iterations of a conjugate gradient algorithm applied to the waterbottom reflection-stripped data. These results show that λ_a is fairly well resolved although it contains cross hatching artefacts and is less clean than the elastic inversion result. Again, the density ρ, which should be constant, has a ghost image of the Lamé modulus λ_a due to imperfect resolvability of ρ from λ as was also seen in the inversion of the diffractors, and the elastic inversions.

Fig. 20 shows the error as a function of iterations and demonstrates that only a fraction of the energy has been matched after the 200 iterations by the acoustic inversion. Again, the zig-zag pattern of error up and down is due to the source-randomization procedure after every second iteration. Fig. 21 shows the initial residual, the initial supershot gather with the seafloor bottom reflection stripped, compared to the residual after 200 iterations. This plot shows that the residual is far from zero and hence, the inversion has not matched the data in the supershot record well. The high residual is due to mode converted energy, as was the case for the diffractors inversion result. Again, the S-wave mode conversions cannot be well matched. As such, artefacts are visible on the inversion result for λ_a.

Figure 20.

The error versus iteration for the acoustic inversion of elastic marine data.

Figure 21.

The initial residuals (a) and the final residuals after 200 iterations (b) for the acoustic inversion of elastic marine data.

6 ACOUSTIC INVERSION OF ACOUSTIC MARINE DATA (FAULT MODEL)

In this section, we apply an acoustic inversion to data calculated in an acoustic model, cf. the last section was an acoustic inversion of elastic data and hence, was intended as an idealized acoustic inversion applied to real data, which by definition, is elastic. The reason for doing this last artificial case is for completeness, and to resolve whether the acoustic inversion has some intrinsic limitation. Specifically, are the artefacts due to doing an acoustic inversion such as seen in previous sections when the data is calculated with the same equation (acoustic wave equation) as is used to do the inversion. Or are the artefacts due in some way relatedto the use of acoustic waves?

Figs 22–25 show results of an acoustic inversion of acoustic data after 200 iterations. Namely, acoustic waves were simulated in the P-wave modulus version of the fault model (see Table 1), and this acoustic data was then inverted using acoustic inversion (i.e. acoustic simulations were used to calculate the gradient in M = λ_a and ρ). Fig. 22 shows that the image for λ_a is now sharp and is of similar quality to the λ_e result of the elastic inversion of elastic data (see Fig. 11c). Hence, the artefacts and lower quality images for the acoustic inversions of elastic data are entirely due to using the wrong wave equation to invert the data. When we say wrong, we mean the approximate acoustic wave equation that is unable to simulate all elastic wave phenomena in the data such as P−S mode conversions.

Fig. 24 shows the error versus iteration for the acoustic inversion of acoustic data, and Fig. 25 shows the initial and final data residuals. Whereas in the acoustic inversion of elastic data, a significant amount of energy could not be matched, the error of the acoustic inversion of acoustic data converged to almost zero after the 200 iterations (Fig. 24), and the data residuals became small (of similar amplitude as for the case of elastic inversion of elastic data—Fig. 17).

Based on these results, the lower quality inversion results and artefacts seen in previous sections on acoustic inversion of elastic data are entirely due to the presence of mode-converted energy that cannot be modelled by the acoustic wave equation. It is not in any way due to any intrinsic limitation of acoustic inversion.

Figure 22.

The result of an acoustic inversion of acoustic data in the 2-D fault model showing the P-wave modulus M = λ_a after 200 iterations. This result is of comparable quality to the result of elastic inversion of elastic data shown in Fig. 11.

Figure 23.

The acoustic inversion result for density ρ of the acoustic marine data after 200 iterations. The true and starting models for density were both constant with ρ = 1 and hence, are not shown here.

Figure 24.

The error versus iteration for the acoustic inversion of acoustic marine data.

Figure 25.

The initial residuals (a) and the final residuals after 200 iterations (b) for the acoustic inversion of acoustic marine data.

7 SOFT SEDIMENTS EXAMPLE

This section aims to study the effect of a layer of soft sediments on hard waterbottom examples studied so far. The question being posed was ‘Will the results be the same when a soft sediment layer exists where potentially less converted S-wave energy occurs at the seafloor?’. The model used was identical to the fault model of previous sections except that a layer under the acoustic (water) layer was added with values of λ = 1 and μ = 0.25², compared to the upper rock layer with λ = μ = 1. Hence, the S-wave velocity of the layer representing soft sediments was 0.25 × V_s(rock). This layer was seven gridpoints thick or about one wavelength, and hence, is intended to represent a narrow layer of soft sediments with a low shear wave velocity. While this may be an adequate representation of a soft sediment layer in terms of P- and S-wave velocities, it does not have attenuation that may be expected of soft sediments, and is therefore potentially not accurately representative.

Fig. 26 shows the inversion results after 200 iterations for the three cases run. Namely, elastic inversion of elastic data, acoustic inversion of elastic data and acoustic inversion of acoustic data. The three sets of results are almost identical to the results of the previous sections for the hard waterbottom case. Namely, the elastic inversion of elastic data yields two clean images of the λ and μ whereas the acoustic inversion of elastic data yields one less clean image of λ with some artefacts. The final image shown is for acoustic inversion of acoustic data and this image is shown to be clean, demonstrating that the cause of the artefacts in the case of acoustic inversion of elastic data is the mode-converted waves.

Figure 26.

The inversion results after 200 iterations for the soft sediments case. Namely, the result for λ (a) and μ (b) for the elastic inversion of elastic data, the result of acoustic inversion of elastic data (c) and the result for acoustic inversion of acoustic data (d).

8 DISCUSSION

The results presented in the previous sections demonstrate that when applying acoustic inversion to elastic data, artefacts are severe for diffraction data, whereas artefacts for layer reflections are minor. We interpret this as follows. Consider a layer result as the sum of diffractor results where diffractors are placed at all points in a layer. Hence, consider a thin reflector result as the sum of the diffractor results where the diffractors are shifted in x. In this case, one will obtain a coherent signal at the depth of the thin reflector. However, the artefacts, which are oscillatory events that are similar to migration ‘smiles’, are summed but these largely cancel out due to having a zero mean. A layer result would then be the sum of thin reflector results where the depth of the thin reflector varies from z_reflector to z_{lower layer interface}. The result would be a coherent signal at the reflector depth. Based on the results of this paper and the above interpretation, one can expect that artefacts due to applying acoustic inversion to elastic data of a model containing layer interfaces will be relatively minor. However, when diffractions are present, such as at fault boundaries and other strong lateral heterogeneities in the model, severe artefacts will be present.

9 CONCLUSIONS

Several numerical experiments have been conducted to study the implications of applying acoustic inversion to elastic marine data. We found that it was necessary to strip the seafloor reflection from the simulation data in order for acoustic inversion to converge to a reasonable quality subsurface image. In cases where the seafloor was not stripped, all effort of the inversion went in trying to match the seafloor reflection which was not possible. Namely, because the seismic waves entering the solid model underwent mode conversions, the seafloor reflection could never be matched by an acoustic simulation. However, by ‘stripping’ the first seafloor reflection, it was possible to obtain convergence for the subseafloor layers. The stripping step means that first, the seafloor must be accurately modelled by an elastic simulation, and then subtracted from the data. Similarly, the acoustic simulations of the inversion must also subtract the acoustic seafloor reflection from the acoustic simulated data. The following results apply to comparisons of an elastic inversion to a seafloor reflection-stripped acoustic inversion of elastic data. It can be expected that clear and strong mode-converted signals occur.

In the case of diffraction signals, an acoustic inversion of elastic data is likely to give a very poor quality results consisting mainly of artefacts after several iterations, even for the case of a marine survey where an upper water layer is present. For the case of data consisting largely of reflections from layer interfaces, we have shown that acoustic inversion of elastic data of a marine survey can give adequate results for a single parameter inversion, although a significant fraction of the data may be left unmatched corresponding to modes such as P−S−S−P (i.e. a P-wave in the water converting to a downgoing S-wave in the subseafloor solid layer, reflecting back as an S-wave to be converted back to a P-wave in the water layer). However, a much higher quality image, and of at least two parameters, λ and μ or equivalent, can be obtained using fully elastic wavefield inversion, which matches the remaining data consisting of such mode-converted subseafloor reflections.

In summary, for the case of large offset marine data in the case of both a hard and soft waterbottom, there is a significant amount of mode-converted energy. This leads to artefacts in the acoustic inversions due to the inversion attempting to fit the mode-converted energy. These artefacts are minor for the case of layer reflections, and are severe for diffractions. In contrast, elastic inversion leads to a clean image, and the resolution of two elastic parameters such as λ and μ or V_p and V_s.

Based on these results, we recommend that elastic inversion need not be carried out unless it has been shown for the specific survey that acoustic inversion gives satisfactory results (i.e. when artefacts are low enough in amplitude and when multiparameter inversion is not beneficial), or that practical issues such as high attenuation of shear waves in soft sediments remove the benefits that may be obtained by elastic inversion, or when elastic computations are too expensive despite the benefits. These conclusions are based on idealized data without noise for both the case of a hard and soft waterbottom. And we do not take into consideration issues that may arise in practice such as the possibility of high attenuation of shear waves in seafloor soft sediments. Nonetheless, the results suggest that one cannot assume an acoustic inversion is adequate due to it having matched most of the data as it may have matched the data by introducing artefacts into the solution, rather than having found the correct solution. However, if there are residuals left after an acoustic inversion corresponding to lower than P-wave velocity events, then it is possible or likely that these residuals originate from data due to mode conversions. In such cases, elastic inversion may be able to match these residuals, and in so doing, obtain a second image for the shear modulus μ or equivalent parameter such as v_s. As these results pertain to idealized theoretical data for the cases of a hard and soft waterbottom, subsequent research is required to validate these results for the practical constraints of real data, and to study how the conclusions of this study apply to more complex and realistic models.

In summary, the synthetic results show that acoustic inversion of large offset marine data leads to artefacts that are severe for diffractions and notable but more minor for layer reflections. The artefacts are the result of the inversion attempting to match the mode converted energy at large offsets with the acoustic wave equation. In contrast, elastic inversion provides two clear images of the two elastic parameters, this being an improvement in both quality and quantity, for a relatively modest (less than about an order of magnitude) increase in computational cost.

ACKNOWLEDGEMENTS

We would like to thank John Etgen for reading the manuscript and helpful suggestions which greatly improved the quality of the paper. We also thank the two reviewers (Milena Marjanovic and anonymous) and the editor (Rene-Edouard Plessix) whose helpful comments led to many improvements in the paper.

APPENDIX A: NUMERICAL SIMULATION METHODOLOGY

In the calculations in this paper, we solve the elastic wave equation on a finite-difference mesh using a staggered grid formulation as described in Igel et al. (1995). The use of a staggered grid allows us to model pressure waves in an acoustic water layer using the same numerical scheme that is used for the elastic waves in the solid layers below the waterbottom. Specifically, we need only set the second Lamé parameter μ = 0 in order to simulate an acoustic medium such as a water layer in which shear stresses are null and hence, shear waves do not propagate.

In this study, we define the grid spacing to be unity, and the Lamé parameters are also set to unity, or variations relative to unity. The use of such unitary parameters has no effect on interpretation of results which only depends on relative scales and dimensionless numbers such as the ratio of offset relative to depth and ratio of the model width to model depth.

A1 Elastic wave equation

We require the ability to forward model seismic waves in a 2-D medium. The equation for the conservation of linear momentum (i.e. F = ma) in a continuum is given by

\begin{equation*} \rho \ddot{u}_i \ = \ {{\partial \sigma _{ij}} \over {\partial x_j}} \ + \ f_i , \end{equation*}

(A1)

where σ_ij is the internal stress tensor in the medium, u_i is the displacement and f_i is the body force source term. In an elastic medium, the stress tensor is related to the strain tensor ε_ij through the generalized Hooke’s Law given by

\begin{equation*} \sigma _{ij} \ = \ c_{ijk\ell } \epsilon _{k\ell } , \end{equation*}

(A2)

where the Einstein convention is used for repeated indices in the above equations (i.e. repeated indices are summed). For the case of infinitesimal strains, the strain tensor is given by

\begin{equation*} \epsilon _{ij} \ = \ {1 \over 2} \left( {{\partial u_i} \over {\partial x_j}} + {{\partial u_j} \over {\partial x_i}}\right) . \end{equation*}

(A3)

For the case of an isotropic elastic medium, we have the following expression for the Hooke tensor (Aki & Richards 2002)

\begin{equation*} c_{ijk\ell } \ = \ \delta _{ij} \delta _{k\ell } \lambda \ + \ (\delta _{ik} \delta _{j\ell } + \delta _{i\ell } \delta _{jk}) \mu , \end{equation*}

(A4)

where δ_ij is the Kronecker delta, that is, δ_ij = 1 if i = j, and δ_ij = 0 if i ≠ j.

To propagate elastic waves in a 2-D medium, one must first calculate the three independent components of the strain tensor using eq. (A3), namely,

\begin{equation*} \begin{aligned} \epsilon _{00} \ & = \ {{\partial u_0} \over {\partial x_0}} , \ \ \ \\\epsilon _{11} \ & = \ {{\partial u_1} \over {\partial x_1}} , \ \ \ \\\epsilon _{01} \ = \ \epsilon _{10} \ & = \ {1 \over 2} \left( {{\partial u_0} \over {\partial x_1}} + {{\partial u_1} \over {\partial x_0}} \right) , \end{aligned} \end{equation*}

(A5)

where the indices 0 and 1 are, respectively, used to denote the x-axis and the z-axis. Subsequently, one must calculate the stress tensor σ_ij using eqs (A2) and (A4), namely,

\begin{equation*} \begin{aligned} \sigma _{00} \ & = \ (\lambda + 2 \mu ) \epsilon _{00} \ + \ \lambda \epsilon _{11} , \ \ \ \\\sigma _{11} \ & = \ (\lambda + 2 \mu ) \epsilon _{11} \ + \ \lambda \epsilon _{00} , \ \ \ \\\sigma _{01} \ & = \sigma _{10} \ = \ \mu ( \epsilon _{01} + \epsilon _{10} ) \ = \ 2 \mu \epsilon _{01} . \end{aligned} \end{equation*}

(A6)

We can then expand using Taylor series, the expression given by

\begin{equation*} \begin{aligned} {{u_i (t+ \Delta t ) - 2 u_i (t) + u_i (t- \Delta t ) } \over {\Delta t^2 }} \ & = \left( {1 \over \Delta t^2} \right) \left[ \sum _{n=0}^N \left({ (\Delta t)^n \over {n!}} {{\partial ^n u_i} \over {\partial t^n}}\right) \ - \ 2 u_i ( t ) \ + \ \sum _{n=0}^N \left({ (-\Delta t)^n \over {n!}} {{\partial ^n u_i} \over {\partial t^n}}\right) \right] \\& = \left( {2 \over {\Delta t^2}} \right) \sum _{m=1}^M \left({ \Delta t^{(2m)} \over {(2m)!}} {{\partial ^{(2m)} u_i} \over {\partial t^{(2m)}}}\right) . \\\end{aligned} \end{equation*}

(A7)

One can then rearrange eq. (A7) to obtain the explicit finite-difference scheme to order 2M given by

\begin{equation*} u_i ( t + \Delta t ) \ = \ 2 u_i (t) \ - \ u_i (t- \Delta t ) \ + 2 \sum _{m=1}^M \left({ \Delta t^{(2m)} \over {(2m)!}} {{\partial ^{(2m)} u_i} \over {\partial t^{(2m)}}}\right) , \end{equation*}

(A8)

where the (2m)-th time derivative of the wavefield u_i can be calculated by applying the wave equation given by eq. (A1) m times. Namely, one first calculates the strains ε_ij from the displacements using eq. (A3). Subsequently, one applies Hooke’s Law given by eq. (A2) to compute the stress tensor σ_ij. Finally, one can calculate the second time derivative of the displacement via the wave equation given by eq. (A1).

To calculate the (2m)th time derivative of the displacement, one repeats this sequence (m − 1) additional times. Namely, one feeds the second time derivative of the displacement into the strain-displacement relation given by eq. (A3) to compute the second time derivative of the strains. Then, application of Hooke’s Law given by eq. (A2) yields the second time derivative of the stress tensor |$\ddot{\sigma }_{ij}$|⁠. Finally, one can calculate the 4-th time derivative of displacement by applying the wave equation given by eq. (A1) to the second time derivative of the stress tensor |$\ddot{\sigma }_{ij}$|⁠. Hence, the even time derivatives to any order can be computed to the precision of the spatial derivative operator. If we take this series to second order, one obtains the second order explicit finite-difference extrapolation equation given by

\begin{equation*} \begin{aligned} u_0 ( t+\Delta t) \ & = \ 2 u_0 ( t) \ - \ u_0 ( t - \Delta t) \ + \ {{\Delta t^2} \over {\rho }} \left( {{\partial \sigma _{00}} \over {\partial x_0}} \ + \ {{\partial \sigma _{01}} \over {\partial x_1}} \right) , \\u_1 ( t+\Delta t) \ & = \ 2 u_1 ( t) \ - \ u_1 ( t - \Delta t) \ + \ {{\Delta t^2} \over {\rho }} \left( {{\partial \sigma _{11}} \over {\partial x_1}} \ + \ {{\partial \sigma _{01}} \over {\partial x_0}} \right) . \\\end{aligned} \end{equation*}

(A9)

In this paper, we use a 12-point convolutional operator to compute space derivatives (Mora 1986; Igel 2017) and we calculate the time derivatives to second order (i.e. M = 1 in eq. A8). Namely, the spatial derivatives on the staggered grid are given by the pseudo-spectral operator ik_xexp (ik_xΔx/2) multiplied by a Gaussian function (see Mora 1986; Igel et al. 1995). For efficiency, these convolutional spatial derivatives are computed by Fourier transforming the displacement wavefields, followed by multiplication by the staggered derivative operator and inverse Fourier transforming the result, where the Fourier transforms required by this pseudo-spectral spatial derivative operator are calculated using the parallel fftw library called from a python source code (see https://pypi.python.org/pypi/pyFFTW for fftw and www.python.org for python). All calculations were done on a small home workstation, and these limited computational resources led to restrictions of the size of the calculations to 2-D finite-difference meshes not larger than around 502 × 232 units.

APPENDIX B: INVERSION FORMULATION

The inverse problem is expressed in terms of finding the most probable model m for a given set of data observations d₀ (e.g. Tarantola 2005; Menke 2012). Consider the Gaussian probability density function

\begin{equation*} P \ \propto \ \exp \left[ - {1 \over 2} \left( \Delta {\bf d}^T {\bf C_d^{-1}} \Delta {\bf d} + \Delta {\bf m}^T {\bf C_m^{-1}} \Delta {\bf m} \right) \right] \ = \ \exp \left( - E \right) , \end{equation*}

(B1)

where C_d is the data covariance matrix, C_m is the model covariance matrix and

\begin{equation*} \Delta {\bf d} = {\bf d} ( {\bf m} ) - {\bf d}_0 , \ \ \ \ \Delta {\bf m} = {\bf m} - {\bf m}_0 , \end{equation*}

(B2)

are, respectively, the data mismatch or residual defined as the difference between the calculated data d(m) and observed data d₀, and the model mismatch between the model m and an a priori model denoted m₀. The most probable or maximum likelihood solution, when the probability P is maximal, occurs when the derivative of the square error functional

\begin{equation*} E \ = \ {1 \over 2} \left( \Delta {\bf d}^T {\bf C_d^{-1}} \Delta {\bf d} + \Delta {\bf m}^T {\bf C_m^{-1}} \Delta {\bf m} \right) , \end{equation*}

(B3)

is minimized. Hence, we seek a solution when the derivative of the error functional or gradient vector g = ∂E/∂m^T = 0. Namely, we seek the solution to

\begin{equation*} {\bf g} \ = \ {{\partial E} \over {\partial {\bf m}^T}} \ = \ {\bf D}^T {\bf C_d^{-1}} \Delta {\bf d} + {\bf C_m^{-1}} \Delta {\bf m} \ = \ 0 , \end{equation*}

(B4)

where D = ∂d/∂m.

In the following, for simplicity in the subsequent idealized inversions, we assume that the model parameters are unconstrained so |${\bf C}_{\bf m}^{-1} = 0$|⁠, and the data errors are uncorrelated and constant so |${\bf C}_{\bf d}^{-1} = {\bf I}/\sigma _\mathrm{ d}^2$|⁠. In this case, we need to solve

\begin{equation*} {\bf g} \ = \ {\bf D}^T \Delta {\bf d} \ = \ 0 . \end{equation*}

(B5)

This can be achieved by a conjugate gradient algorithm, provided the gradient can be calculated.

For isotropic elastic media, the equations for the gradient are given by (see Mora 1987b; Tarantola 1988)

\begin{equation*} \begin{aligned} \delta \hat{\rho }\ & = \ \ \ \ \sum _S \int \mathrm{ d}t \ \ \ {{\partial u_i} \over {\partial t}} {{\partial \psi _i \over \partial t}} \ \ \ \ \ = \ \ \ \ \sum _S \int \mathrm{ d}t \ \ \ \dot{u}_i \dot{\psi }_i , \\\delta \hat{\lambda }\ & = \ - \sum _S \int \mathrm{ d}t \ \ \ {{\partial u_i} \over {\partial x_i}}{{\partial \psi _j} \over {\partial x_j}} \ \ \ \ = \ - \sum _S \int \mathrm{ d}t \ \ ( \nabla _i u_i ) ( \nabla _j u_j ) \ \ \ \ = \ \ \ - \sum _S \int \mathrm{ d}t \ \ \ \epsilon _{ii} \epsilon _{jj}^\dagger , \\\delta \hat{\mu }\ & = \ - \sum _S \int \mathrm{ d}t \ \ {1 \over \sqrt{2}} \left({{\partial u_i} \over {\partial x_j}} + {{\partial u_j} \over {\partial x_i}} \right) {1 \over \sqrt{2}} \left({{\partial \psi _i} \over {\partial x_j}} + {{\partial \psi _j} \over {\partial x_i}} \right) \ \ \ \ \ = \ \ \ - \sum _S \int \mathrm{ d}t \ \ 2 \epsilon _{ij} \epsilon _{ij}^\dagger , \end{aligned} \end{equation*}

(B6)

where the wavefield ψ_i(x, t) is the back propagated data residuals, ε_ij denotes the strain tensor of wavefield u_i, |$\epsilon _{ij}^\dagger$| denotes the strain tensor of the adjoint wavefield ψ_i and the Einstein convention is used for repeated indices (i.e. repeated indices are summed). Note that the equations for the general case in terms of the density ρ, the Hooke tensor coefficients c_ijkℓ and the attenuation are provided by Tarantola (1988). In the following, we denote the gradient vector at the n-th iteration as g_n where

\begin{equation*} {\bf g} \ = \ \left[ \delta \hat{\rho }, \delta \hat{\lambda }, \delta \hat{\mu }\right]^T . \end{equation*}

(B7)

B1 Conjugate gradient method

The conjugate gradient method involves iterations updating the model in the conjugate gradient direction c_n as follows

\begin{equation*} {\bf m}_{n+1} \ = \ {\bf m}_n \ - \ \alpha _n {\bf c}_n , \end{equation*}

(B8)

where α_n is a scalar step-length and c_n is used to denote the conjugate direction at the n-th iteration which is calculated as a linear combination of the current gradient vector and the last conjugate direction. Namely,

\begin{equation*} {\bf c}_n \ = \ {\bf g}_n + \beta _n {\bf c}_{n-1} , \end{equation*}

(B9)

where β_n is a scalar which can be calculated according to a number of different formulae which will be specified in the following. The benefit of the conjugate gradient method is that it avoids the zig-zag behaviour suffered by the steepest descent method as it progresses down a banana-shaped objective function.

The steplength α_n can be calculated by a line search, with an initial estimate of α_n calculated using the following formula

\begin{equation*} \alpha \ = \ { {{\bf c}^T {\bf c} } \over { {\bf c}^T {\bf D}^T {\bf C}_{\bf d}^{-1} {\bf D} {\bf c} + {\bf c}^T {\bf C}_{\bf m}^{-1} {\bf c} }} , \end{equation*}

(B10)

which can be derived by solving ∂ε/∂α_n = 0 (see Mora 1987a).

The conjugate gradient method starts by an iteration in the direction of steepest descent (i.e. β₁ = 0) and subsequently, iterations are calculated according to eq. (B9). After the first iteration, the value of the scalar β_n can be calculated by various formulae which are named after their developers. Four popular choices are the Fletcher–Reeves, Polak–Ribiére, Hestenes–Steifl and Dau–Yuan methods (denoted as FR, PR, HS and DY). Each of these formulae yield the same result for quadratic objective functions (linear theories), and result in the solution after exactly n iterations where n is the number of unknowns or model parameters being solved for. However, the performance of these different methods to calculate β_n varies. We have tested each of these popular choices for the examples of FWI shown in this paper, and found that for the cases studied that the performance (in order of descending performance) is PR, FR, DY and lastly, HS. The PR and FR methods were similar in performance with PR being slightly better, and these were both significantly better in performance than the DY and SH methods which were also similar in performance. For this reason, we use the PR method given by

\begin{equation*} \beta _n^{PR} \ = \ { {{\bf g}_n^T ( {\bf g}_n - {\bf g}_{n-1} ) } \over {\bf g}_{n-1}^T {\bf g}_{n-1}} . \ \ \ \end{equation*}

(B11)

At some stage, progress using the conjugate gradient method may slow or stop due to loss of conjugacy for nonlinear functions (i.e. non-quadratic objective functions). If this occurs, a reset to the direction of steepest descent will allow subsequent iterations to regain conjugacy. A popular choice for β that provides an automatic direction reset and, provides the best performance for the examples in this paper is given by

\begin{equation*} \beta \ = \ \max \left\lbrace 0,\beta ^{PR} \right\rbrace . \end{equation*}

(B12)