ABSTRACT

Full-waveform inversion (FWI) is a technique used to obtain high-quality velocity models of the subsurface. Despite the elastic nature of the earth, the anisotropic acoustic wave equation is typically used to model wave propagation in FWI. In part, this simplification is essential for being efficient when inverting large 3D data sets, but it has the adverse effect of reducing the accuracy and resolution of the recovered P-wave velocity models, as well as a loss in potential to constrain other physical properties, such as the S-wave velocity given that amplitude information in the observed data set is not fully used. Here, we first apply conventional acoustic FWI to acoustic and elastic data generated using the same velocity model to investigate the effect of neglecting the elastic component in field data and we find that it leads to a loss in resolution and accuracy in the recovered velocity model. Then, we develop a method to mitigate elastic effects in acoustic FWI using matching filters that transform elastic data into acoustic data and find that it is applicable to marine and land data sets. Tests show that our approach is successful: The imprint of elastic effects on the recovered P-wave models is mitigated, leading to better-resolved models than those obtained after conventional acoustic FWI. Our method requires a guess of VP/VS and is marginally more computationally demanding than acoustic FWI, but much less so than elastic FWI.

INTRODUCTION

Obtaining accurate velocity models of wave propagation in the subsurface is key to obtain high-quality images at the depth of reservoirs. Full-waveform inversion (FWI) is a model building technique that is increasingly being used in the industry to obtain highly accurate subsurface velocity models (Sirgue et al., 2010; Warner et al., 2013). This nonlinear optimization method requires modeling the data at the receiver positions and minimizing its difference with the observed data in the field in an iterative fashion. Hence, wave propagation in the subsurface needs to be computed.

In the ideal case, the full content of the recorded signal is used and a wave equation that accounts for all types of generated waves is used to model wave propagation (Virieux and Operto, 2009). Properly addressing the kinematics and dynamics of the wavefield by accounting for the anisotropic and viscoelastic nature of the earth is expected to lead to high-resolution P-wave velocity subsurface models and the recovery of other subsurface properties such as S-wave velocity, given that the match between observed and modeled data is expected to improve (Warner et al., 2012). But, elastic FWI is not yet widely used due to (1) the high computational cost of elastic modeling (Guasch et al., 2010; Hobro et al., 2014), (2) the nonunique and coupled effect of the different physical properties of the wavefield (Virieux and Operto, 2009), and (3) the difficulty in modeling amplitudes produced by inaccuracies in the source wavelet and often poorly known subsurface density, attenuation, and S-wave velocity (Tarantola, 1986).

In practice, however, several simplifications are made mainly to make the computation more efficient. Typically, the anisotropic acoustic wave equation is used to model wave propagation, especially when inverting large 3D data sets (Warner et al., 2013). Anisotropic acoustic FWI performs well because it accounts for the correct kinematics of the wavefield, but it is adversely affected by the elastic component in the data because it does not account for phase conversions at interfaces (P- to S-wave conversions in an isotropic earth), which affect the amplitude of all arrivals, including first-arriving P-waves (Hobro et al., 2014), leading to a reduction in accuracy and resolution of the recovered P-wave velocity models in FWI.

Some strategies have been suggested to address these problems in FWI. Chapman et al. (2010) propose a method to correct acoustic simulations for some of the effects of elasticity at a lower cost than performing fully elastic simulations by adding a correction term to the acoustic wave equation. Although this method can be used to correct the amplitude of reflected P-waves in an approximate fashion without having to perform a full elastic simulation on multicomponent data, it can degrade the result of acoustic simulations in the presence of strong velocity contrasts (Hobro et al., 2014). Other approaches to mitigate elastic effects consist of using an objective function that focuses more on the phase of the recorded data and less on the amplitudes, in which the observed and modeled data are weighted by their root-mean-square (rms) energy trace by trace (Shen, 2010). This strategy is frequently used in the industry, but it is insufficient in the presence of strong elastic effects.

Here, we propose an alternative method to mitigate a wider range of elastic effects when performing acoustic FWI based on the use of matching filters to match elastic to acoustic data, which requires a guess of the VP/VS model and which is applicable to pressure and multicomponent data. This method does not aim at obtaining an S-wave model of the subsurface, but at improving the quality and resolution of the recovered P-wave models obtained with conventional acoustic FWI. Hence, it does not suffer from crosstalk between different parameter classes because we only invert for P-wave velocity. Nevertheless, an approximate S-wave velocity model can be obtained using this strategy, the quality of which relies on the guess of VP/VS and improvement in the recovered P-wave velocity model. This method can also be combined with the amplitude balancing method of Shen (2010) to further mitigate elastic effects.

Within this paper, we demonstrate how not considering the elasticity of the earth adversely affects the recovered P-wave velocity models obtained with acoustic FWI and we suggest a workflow to mitigate these elastic effects. We show that the suggested method improves the performance of acoustic FWI in marine and land synthetic data sets in a computationally efficient way without the need for an accurate model of VP/VS or a full elastic FWI scheme. Finally, we discuss the benefits and limitations of the current method, as well as further applications of this scheme to account for other effects not typically considered in acoustic FWI.

METHODOLOGY

In acoustic FWI, the misfit between observed and modeled data using the acoustic wave equation is minimized iteratively to recover a velocity model of the subsurface. Given that the acoustic equation is used, S-waves and P- to S-converted waves are not properly modeled, which means that the amplitudes of P-wave arrivals are also incorrect. This is illustrated in Figure 1, in which a snapshot of the wavefield in a horizontally layered model with three layers is shown for an acoustic (Figure 1a) and an elastic media (Figure 1b).

In both media, a pressure source is fired at the center of the model and the P-wavefield initially propagates through the water layer, for which the S-wave velocity is set to zero. After reaching the fluid-solid interface, part of the wavefield is reflected toward the receivers and part of it is transmitted to the second layer. A critically refracted P-wavefield is also generated, which propagates toward the receivers with the P-wave velocity of the second layer. In the elastic case, surface waves at the water-rock interface (Scholte waves) are produced, which propagate within this zone with the S-wave velocity of the second layer and have an amplitude that decreases exponentially away from the interface. The latter are difficult to observe in a single snapshot, but they can be noticed if multiple snapshots are concatenated. Due to the differences in the splitting of the wavefield in the acoustic and elastic media, the relative amplitudes of the reflected and transmitted wavefields are different. At a later time-step, the transmitted P-wave field in the second layer hits the second interface, part of the wavefield is transmitted to the bottom layer, and the remaining part is reflected upward. We note the difference in amplitude variation with offset of the reflected P-wave in both media. Moreover, a converted P- to S-wave is generated at the second interface in the elastic case, which propagates with the S-wave velocity of the second layer and does not exist in the acoustic medium. Thus, the recorded modeled data at the receivers differ in both media. Because FWI works by reducing the misfit between modeled and observed data, its performance will be affected if the full viscoelastic wavefield is not correctly modeled.

To transform elastic to acoustic data when performing acoustic FWI, we suggest the workflow depicted in Figure 2, based on the use of matching filters. The first step is to perform conventional acoustic FWI of the observed (elastic) data to obtain a recovered P-wave model (step 1 of our workflow). We follow the implementation of FWI described in detail by Warner et al. (2013) and we neglect attenuation and anisotropy here for simplicity. The recovered P-wave velocity model and a guess of VP/VS are then used to generate acoustic and elastic modeled data (step 2), and matching filters that match the predicted elastic data to the predicted acoustic data are computed. Provided that the recovered P-wave velocity model and the guess of VP/VS are close enough to the true models, these filters can then be used to mitigate the elastic effects by convolving them with the observed data (step 3). An additional acoustic inversion of the matched observed data then leads to an improved P-wave velocity model in which the effects of elasticity are reduced (step 4). To further improve the result, the process can be repeated again (step 5) by first modeling acoustic and elastic data using the same VP/VS and the improved recovered P-wave velocity model and recomputing the matching filters that match the elastic to the acoustic modeled data. Convolution of these filters with the observed data and the posterior acoustic FWI step should then lead to an improved P-wave velocity model. This process can be repeated several times, with each additional iteration in the workflow incrementing the computational burden and run time. The combination of using an amplitude-balanced objective function during the inversion steps — e.g., based on the method of Shen (2010) — and the current workflow leads to optimal recovered P-wave velocity models in the presence of strong elastic effects.

The computed matching filters are treated as optimum Wiener filters Wiener (1949) given that these can easily be used to convert any input signal into a desired output signal (Karsli and Bayrak, 2004). To match the modeled elastic data (del, input signal) to the modeled acoustic data (dac, output signal), we define a least-squares error function L between the two as follows: 
L1=dac,1del,1*w22=12t=0N+P2(dac,1(t)k=0P1del,1(tk)w(k))2,
(1)
where the subscript 1 indicates the error function is computed for a single trace, w is a Wiener filter, N is the number of time samples per trace, P is the length of the filter, and * denotes convolution. Minimization of the error function with respect to the coefficients w(i) leads to a set of normal equations that can be expressed in terms of an autocorrelation matrix—which is a Toeplitz matrix — and a cross-correlation vector, which need to be solved to find the filters’ coefficients (Karsli and Bayrak, 2004).
To impose smoothness of the matched data and mitigate artifacts of trace-by-trace filtering (e.g., matching two traces with different spectra may lead to an over-amplification of a notch), we constrain the solution so that a single Wiener filter matches several traces from both data sets simultaneously, which are spaced a fixed distance proportional to the receiver spacing. This is done as follows: 
Lj=1MLj=12j=1Mt=0N+P2(dac,j(t)k=0P1del,j(tk)w(k))2,
(2)
where M is set to be an odd number and it represents the number of traces to be matched, being trace (M+1)/2, the central trace in the observed data set to which the matching filter is later on applied in the workflow in Figure 2. The previous equation assumes that variations between neighboring traces are small and, therefore, introduces some degree of smoothness without compromising the computation time — the latter mainly depends on the length of the filter, which is unchanged. Minimization of this error function with respect to the filter coefficients leads to the following system of normal equations: 
C(i)=k=0P1w(k)R(ki),i=0,,P1,
(3)
in which, C and R are the sum of cross correlations of the input and output signals for the different traces and the sum of autocorrelations of the input signal for all traces involved, respectively. The Levinson's recursion method (Levinson, 1947) is used to solve this system of equations. The number of traces used in this computation is linearly reduced close to the edges of the data set; i.e., a trace-by-trace strategy is used for the first and last traces of a shot gather. Additionally, the matched data are smoothed along each trace by using overlapping short filters within a trace weighted and stacked along each trace by using a Blackman window such that the sample at the center of a short filter has a larger weight than the same sample of the next overlapping short filter.

The workflow suggested in Figure 2 also requires a guess of VP/VS. In the next section, we validate our approach on synthetic models and we test the impact of using different VP/VS models.

RESULTS

We now test the performance of our workflow on three different synthetic elastic pressure data sets. For all the examples presented in this manuscript, rock density is linked to P-wave velocity and follows Gardner’s law (Gardner et al., 1974). During the inversions, density is updated at each iteration using Gardner’s law and the recovered P-wave velocity model for that iteration.

In the first set of tests, we use 2D horizontally layered models. First, we test the effect of acoustically inverting elastic subsurface models with different VP/VS, similar to the study performed by Barnes and Charara (2009) for 1D models, but we run 2D inversions and use horizontally layered media with a high level of horizontal smoothing and a larger number of inverted model parameters. Then, we use our workflow to mitigate elastic effects using the recovered models and the true VP/VS and we test the impact on the matched elastic data when matching trace-by-trace or multiple traces simultaneously. We also investigate the influence of the length of the Wiener filters on the obtained matched observed data. Finally, a limitation of the current workflow when elastic effects are too large and starting models are far from the true model is shown and a solution to the latter is successfully implemented.

The second and third tests are built from the well-known Marmousi2 and overthrust velocity models, respectively. These are designed to represent marine and land elastic data sets, respectively. Elasticity has a weaker impact on the former given that elastic effects occur at some distance from the sources and are due to P- to S-wave conversions and surface S-waves generated at the sea bottom. However, elastic effects are stronger on the land data set as sources are in a region with nonzero S-wave velocity, which results in a data set with a stronger elastic component. The performance of acoustic FWI of acoustic data is compared for these synthetic models with the acoustic FWI of elastic data. The former represents the ideal situation, in which elasticity has no effect, whereas the latter represents the typical situation in field data because the observed data contain elastic effects. Our workflow is then applied to mitigate elastic effects in both synthetic data sets, and we show on the Marmousi2 model that a rough estimate of VP/VS is good enough to mitigate most of the elastic effects.

Horizontally layered 2D models

The true models for the first set of tests are 2D horizontally layered models, which are shown in Figure 3. The first model (model A) consists of two layers with constant P- and S-wave velocities and large P- and S-wave impedance contrasts at the sea bottom (450 m) and at 640 m depth. The second model (model B) has a water depth of 340 m, two constant velocity layers, and two layers with velocity gradients. The impedance contrasts between layers in model B are smaller than for model A, and the maximum P- and S-wave velocities are higher. We consider two S-wave velocity models for each model set to test our workflow when different levels of elastic effects are present in the data. The first S-wave velocity model for both data sets (black lines with triangles in Figure 3) leads to weaker elastic effects when modeling elastic data as impedance contrasts are smaller than for the second S-wave velocity models (gray lines with circles in Figure 3).

The reference models in Figure 3 are used to generate observed acoustic and elastic single-component pressure data using finite-difference acoustic and elastic modeling codes. The same source signature, modeling, and inversion parameters are used for both models. The horizontally layered 2D models are 15 km long by 900 m deep, and a grid spacing of 3 m is used to ensure there is no numerical dispersion of the P- or S-waves for either model. Data are generated using 41 sources deployed at a 6 m depth, with a shot spacing of 300 m, and using a 10 Hz Ricker wavelet as a source signature. The wavefield is recorded by 967 hydrophones deployed at 16 m depth. The maximum offset of approximately 7.25 km ensures the recording of turning rays that have traveled in the deeper high-velocity layers. Absorbing boundaries are used on the sides and the bottom of the model, and a free surface is considered at the top of the model to simulate marine seismic data; hence, sea-surface multiples, shot, and receiver ghosts are present in the data.

First, we perform acoustic FWI of the observed acoustic and elastic observed data sets for model set A using a multiscale approach by inverting frequencies from 3 to 15 Hz in a total of 60 iterations and using the conventional FWI objective function. Figure 4a shows the results for model set A as well as the true (the black line) and starting P-wave (the dashed black line) velocity models at the center of the model. We observe that the recovered P-wave velocity model after acoustic FWI of acoustic data (the green line with squares) fits the true model accurately and is only limited by the range of frequencies inverted, i.e., it cannot fit the sudden jump in velocities at 640 m depth perfectly due to limits in bandwidth, but it can fit the jump at the sea bottom because the starting model is exactly the same as the true model in this region. However, the recovered models after acoustic FWI of both observed elastic data sets for model set A (the blue and red lines with triangles and circles, respectively) are overall less accurate due to elastic effects not being considered. We note the ripple introduced in the model at the sea-bottom interface due to the impact of the S-wave velocity on the impedance contrast at the sea bottom, which is neglected when performing acoustic FWI. This has an effect on the goodness of the recovered models within the first layer down to 640 m depth. Deeper in the model, stronger elastic effects introduced by the second elastic model (the red line with circles) reduce the accuracy of the recovered model further than those for the elastic model with reduced elastic effects (the blue line with triangles).

In a second step, we use the recovered P-wave velocity models after acoustic FWI of the observed elastic data and the true VP/VS obtained from Figure 3 to mitigate elastic effects from the data using the workflow in Figure 2. Acoustic and elastic data are modeled for all the shots prior to any inversion (step 2 in the workflow), a set of matching filters that match the elastic to the acoustic modeled data is computed, and these are convolved with the observed data (steps 3 and 4 in the workflow) to mitigate elastic effects, leading to the matched observed data. Figure 4b shows the result of performing acoustic FWI on both matched observed data sets (the blue and red lines with triangles and circles, respectively). We observe there is a clear improvement in terms of accuracy because the resulting models are closer to the result of acoustic FWI of acoustic data (the green line with squares) than those before applying the current workflow (Figure 4a), and hence closer to the true model (the black line). Despite an improvement, a single iteration of our workflow cannot fully mitigate the strong elastic effects due to poor recovery of the impedance contrast at the sea bottom for the second elastic model (the red line with circles). In such cases, we suggest performing a second iteration of the workflow, which further improves the result (cyan line with stars), but it comes with an increase in computational cost. This is mainly due to running an additional acoustic FWI — but we note that this is still less expensive than a full elastic FWI.

In the tests shown in Figure 4b, matching filters have been computed assuming a fixed length of the filters and the number of neighboring traces to be simultaneously matched. The choice of these parameters has an influence on the resulting matched data sets and, hence, also in the recovered P-wave velocity models. Figure 5 shows the impact of a different combination of filtering parameters for the first elastic model for model set A (the black line with triangles in Figure 3a). Figure 5a shows the shot gather for the observed elastic data, and its difference with the observed acoustic data (no elastic effects) is shown in Figure 5b. Figure 5c, 5e, 5g, and 5i shows the shot gathers for the matched elastic data obtained using different data-matching parameters, whereas panels on the right column in Figure 5 correspond to the difference of the panel on the left and the observed acoustic data. The result in Figure 5c has been obtained matching data trace-by-trace, whereas in Figure 5e, 5g, and 5i, 11 neighboring traces have been matched simultaneously. In terms of the length of the matching filters, this is set to 320 ms for Figure 5c and 5e, to 35 ms for Figure 5g, and to 1.58 s for Figure 5i. By comparison of the panels on the right column in Figure 5, the difference with respect to the observed acoustic data is overall reduced in the matched data sets with respect to the observed elastic data (compare Figure 5d, 5f, 5h, and 5j with Figure 5b) except for some localized areas due to small ripples introduced by the matching filters. The ripples are due to the fact that we compute matching filters using modeled data computed with the recovered P-wave velocity model after acoustic FWI of observed elastic data, not the true P-wave velocity models. Throughout this manuscript, all tests indicate that these small artifacts have a small imprint on the improved recovered P-wave velocity models, even when the true VP/VS is unknown and an estimate is used, given that mitigating elastic effects have a bigger effect on the result. When comparing the different matched shot gathers, we also conclude that matching multiple traces simultaneously reduces horizontally varying artifacts and produces a more continuous and plausible data set (compare Figure 5c and 5e). The length of the filters is also key to obtain a reasonable and continuous data set. When comparing Figure 5e, 5g, and 5i, we note that, although short filters reduce the computation cost of the matching filters (Figure 5g), they produce some artifacts at short offsets and also when there are crossing events. Longer filters (Figure 5i) produce better results at short offsets, but they produce worse results at longer offsets as well as an increase in the computation cost. We have observed that when stronger elastic effects are present in the data or when the true VP/VS is unknown, longer filters can produce noisier matched data sets because the portions of the acoustic and elastic model data to be matched have different spectra. Figure 5e shows that an intermediate filter length is a good compromise in terms of computation time and quality of the matched data.

From the latter and several other tests on different data sets for different models, we recommend using enough neighboring traces so that the spacing between the first and the last trace is at least 150 m apart. In terms of the length of the filters, we conclude that a length of between 300 and 350 ms produces good-quality results. These findings will be used later on to mitigate elastic effects on the realistic marine and land synthetic data sets presented below this section.

Next, we perform acoustic FWI of the observed acoustic and elastic observed data sets for the more realistic model set B using the same inversion strategy and inversion parameters as for model set A. Figure 6a and 6c shows the results for two different starting models: (1) a more accurate starting model within the first layer—where the elastic effects are more important due to the difference on the impedance contrast at the sea-bottom — and a less accurate deeper model (Figure 6a) and (2) a starting model that is less accurate within the first layer but more accurate at the bottom of the model (Figure 6c). In both tests, the inverted P-wave velocity model after acoustic FWI of the acoustic data (the green lines with squares) recovers the true model very accurately. In the presence of elastic effects, the result varies depending on the starting model and the level of elasticity. For a better starting model at the top and less accurate at the bottom (Figure 6a), the recovered models after acoustic FWI follow the main trend of the true model down to 750 m depth, where cumulative elastic effects and a poor starting model lead to inaccurate recovered P-wave velocity models for both elastic data sets (the blue and red lines with triangles and circles, respectively). Similar to model set A, we observe that the recovered P-wave velocity models oscillate close to the sea bottom due to the difference in impedance contrasts between the acoustic and elastic scenarios, which translates into an oscillating and less accurate recovered model within the first layer. Overall, the recovered P-wave velocity model for the first elastic model (the blue line with triangles) is closer to the true model than the recovered P-wave velocity model for the second elastic model (the red line with circles). However, for a less accurate starting model at the top but better at the bottom of the model (Figure 6c), we observe that there is a good match of the recovered P-wave velocity model for the first elastic data set (the blue line with triangles) within all layers with respect to the true model. However, the impact of stronger elastic effects for elastic model 2 combined with a poor starting model for the first layer leads to a very poorly recovered P-wave velocity model (the red line with circles). The latter is an effect of the inversion being trapped in a local minimum when we minimize the misfit between the second observed elastic data set and acoustic data initially generated from a poor starting model, given that the misfit between the two is initially larger due to the presence of strong elastic effects. We can minimize these elastic effects for the second elastic model by using the amplitude balancing objective function described by Shen (2010): The result is shown in Figure 6c (the dashed red line with empty circles). We observe an improvement with respect to the result of acoustic FWI of elastic data using the conventional objective function (the solid red line with circles in Figure 6c), especially with depth, but it is still far from the true model. We now compare this result to that obtained with the current method.

Analogously to model A, we now apply the workflow in Figure 2 to mitigate elastic effects using the recovered P-wave velocity models after acoustic FWI of observed elastic data and the true VP/VS obtained from Figure 3b. Figure 6b and 6d shows the results of acoustic FWI of the matched elastic 1 and 2 data sets (the blue and red lines, respectively) for both starting models. For the first starting model (Figure 6b), a single iteration of our workflow is able to mitigate most of the elastic effects effectively and produce an accurate result for both elastic models. We observe a clear improvement at less than 750 m depth after mitigating the elastic effects, although it is not perfect when the elastic effects are stronger (the red line with circles). To further improve this result, we suggest performing a second iteration of the workflow. However, for the second starting model (Figure 6d), we see there is a limitation of the current method to mitigate elastic effects when these are too large (elastic model 2) and the starting model close to the sea bottom is not accurate (the red line with circles), whereas it is able to perform well when elastic effects are weaker (the blue line with triangles). This happens when the recovered P-wave velocity model after acoustic FWI of the observed elastic data is far from the true model (Figure 6c, the red line with circles). However, we suggest a modified workflow to deal with elastic effects in this situation. The latter consists of skipping step 1 of the workflow in Figure 2 and modeling acoustic and elastic data using the starting P-wave velocity model and the estimate of the VP/VS, therefore avoiding the first acoustic inversion in the workflow in Figure 2. As shown in Figure 6d for the elastic model 2 (magenta line with stars), this provides an accurate result because the modeling of acoustic and elastic data is not linked to the acoustic FWI of observed elastic data, which leads to an inverted result that is less accurate than the starting model due to the strong elastic effects and the starting model being poor at the sea bottom. The cost of the modified workflow is slightly higher than the cost of conventional acoustic FWI, and it is about half the cost of a single iteration of the initial proposed workflow in Figure 2. This result is better than that obtained just by amplitude balancing the objective function (Figure 6c, the dashed red line with empty circles). However, after testing on several 2D data sets, we recommend to use the modified workflow only when strong elastic effects are present in the data and when the starting models are not well-constrained. The initially proposed workflow in Figure 2 still provides better results in most of the situations, as demonstrated in the tests in the next sections.

When compared with the results of Barnes and Charara (2009) for 1D models, we also conclude for 2D horizontally layered models that the results of acoustic FWI of elastic data are not reliable in the marine case when large density or S-wave velocity contrasts are present, in which elastic effects need to be taken into account. Additionally, we have observed that the combination of strong elastic effects and a poor starting model may lead to inaccurate inverted models after acoustic FWI of the elastic data, whereas the recovered P-wave velocity model may be accurate for the same poor starting model and weaker elastic effects. In all cases, our workflow in Figure 2—or a simpler modified workflow, as explained above—can be used to mitigate elastic effects and provide higher quality P-wave velocity models at a lower cost than elastic FWI.

A marine synthetic data set

The workflow in Figure 2 is now implemented on a more realistic marine synthetic data set based on the well-known 2D Marmousi2 model. This is an extension of the Marmousi model in offset and depth, which also includes elasticity. We consider the original true P- and S-wave velocity models (Martin et al., 2006), but we limit the top water layer to be 200 m deep to reduce the computation time and enhance the recording of surface S-waves generated at the sea bottom.

Single-component pressure data using acoustic and elastic modeling are generated from the true P- and S-wave velocity models using an elastic finite-difference code and a grid spacing of 2.5 m, which ensures there is no dispersion of either P- or S-waves (the minimum S-wave velocity is very low, approximately 270  m/s). These are treated as the observed data sets in the inversion. Data are generated using 111 sources at 6 m depth with a source spacing of 150 m and a source wavelet with a useful bandwidth that spans from 2 Hz to approximately 20 Hz. The wavefield is recorded on 3361 receivers at 10 m depth with receiver spacing of 5 m, which are spread uniformly along the model. Absorbing boundaries are used on the edges of the model except for the top layer, in which a free surface is used to represent the water-air boundary. Figure 7a and 7b shows a representative shot gather of the true acoustic and elastic data, whereas Figure 7e illustrates the difference of the true elastic data with the true acoustic data in Figure 7a. As expected in a marine environment, differences are small but nevertheless observable and these are mainly seen as variations in amplitude of the P-waves. Converted P- to S-waves are also present in the data, but these are less easy to identify.

Next, acoustic FWI is separately performed on the acoustic and elastic observed data sets using a multiscale approach from 3 to 18 Hz and using all shots per iteration for a total of 145 iterations. The conventional data misfit objective function (no weighting applied) is used here to study the performance of the current workflow as well as the impact of different parameters and noise on the final result. The starting model in Figure 8b is used in all inversions, which is obtained by smoothing the true model (in slowness) with a Gaussian function of 7.5 m correlation length several times depending on the depth level, such that deeper parts of the model are smoother (i.e., the structure is less resolved) than the upper parts. To assess the difference of the various models with respect to the true model in Figure 8a, we compute the average rms difference within the dotted box plotted in this figure. The computed value is shown on the bottom right corner of each model, such that the smaller the value the closer the model within this region is to the true one (i.e., an average rms of 0.0  m/s would indicate that there is no difference between the current model and the true model).

Figure 8c and 8d shows the results of performing acoustic FWI on the acoustic and elastic data, respectively. Whereas we recover most of the layers in both situations, there is a visible loss in resolution of individual layers, and there is a reduction in accuracy in the recovered velocity values when inverting elastic data with an acoustic code not only for deep layers, but also for intermediate and shallow layers. The computed average rms shows that the model in Figure 8d is less resolved than that in Figure 8c. Thus, elasticity has an adverse effect on the recovered P-wave velocity model.

To mitigate these elastic effects, we now apply the suggested workflow. In step 2 of the workflow (Figure 2), we model acoustic and elastic data using the recovered P-wave velocity model after acoustic FWI of the (true) elastic data in Figure 8d and an S-wave velocity model generated from the latter and a guess of VP/VS. For our initial tests, we use a smooth VP/VS (Figure 9b) obtained from smoothing the true VP/VS (Figure 9a). Matching filters that match the elastic to the acoustic data are computed for 35 traces simultaneously, and these are convolved with the true (observed) elastic data in step 3. In step 4 of the workflow, we acoustically invert the matched observed data and we perform a second iteration of the workflow in step 5 to further improve the result. Figure 7c and 7e shows a shot gather of the matched observed data after a second iteration of the workflow and its difference with the true acoustic data, respectively. In comparison with the true elastic data in Figure 7b and 7d, the matched data are closer to the true acoustic data as the overall differences with respect to the latter are reduced, especially in terms of amplitude of the P-wave refractions and reflections. This leads to an improvement in the recovered P-wave model after acoustic inversion (step 4, second iteration), in which we have used the same inversion parameters as the first acoustic inversion, and the result is shown in Figure 8e. We observe an increase in the resolution of the layers with respect to Figure 8d, both visually and also as indicated by the average rms. The improved P-wave velocity model is not as well-resolved as the result after acoustic FWI of acoustic data in Figure 8c, which is expected given that the latter is the ideal case (i.e., no elastic effects). Further iterations on our workflow are expected to improve the recovered P-wave model, even though each repetition increases the computation time. Moreover, these and other tests suggest that the first iteration of the workflow is the most important and improvements during subsequent iterations are smaller.

To further assess the results, we compute the rms as a function of depth and position within the dotted box for the relevant models in Figure 8 and the results are shown in Figure 10a and 10b, respectively. From this result, we observe that (1) the recovered model after acoustic FWI of acoustic data (the red line) is closer to the true model when compared with the other results for all positions and depths, (2) the P-wave model after acoustic FWI of elastic data (the black line) is not as close to the true model as that after acoustic FWI of acoustic data for all positions and depths except for localized anomalies, and (3) there is an improvement in terms of rms difference when using our workflow (the green line).

Overall, Figures 8c8f and 10a show that differences caused by elastic effects are small for the marine Marmousi2 model down to 1.6 km depth, which can be addressed with the current workflow (Figure 10a and 10b). We conclude that elastic effects for marine data sets with not-too-strong P- to S-wave conversions can be initially ignored in acoustic FWI, and they can be taken into account by using the current workflow when needed, leading to small improvements. Alternatively, the amplitude-balancing method of Shen (2010) might be sufficient in such situations. We expect this method to be more beneficial for marine data sets in which VP/VS is larger at the sea bottom but also for land data sets. We illustrate the latter in the following section, but first we test the impact of the VP/VS guess and the impact of noise on the marine Marmousi2 data set and we compare the resulting P-wave velocity models with that obtained after elastic FWI of the elastic data.

Impact of VP/VS guess

The workflow in Figure 2 requires a guess of VP/VS, which is unlikely to be known in most field examples. Thus, we now test the importance and the impact of using an alternative VP/VS model on the marine Marmousi2 data set. The model used in the initial inversions, shown in Figure 9b, was obtained by mildly smoothing the true ratio. This represents the case that we know VP/VS quite well. The second one, much poorer in detail, is a 1D average of the true model (Figure 9c), which is more accurate close to the sea-bottom to capture the relevant elastic effects generated in this region. The former leads to the improved recovered P-wave velocity model already shown in Figure 8e, whereas the latter leads to the model in Figure 8f. We observe that both recovered models are more resolved and closer to the true model than that obtained after acoustically inverting elastic data (Figure 8d), as indicated by the average rms value plotted on the figure, and that the recovered model using a smooth true VP/VS model leads to a slightly better result. However, using the 1D average, VP/VS is sufficient to mitigate the most important elastic effects because we capture the average elastic components at the sea bottom, where the strongest elastic effects are generated (due to the change in VP/VS at the sea bottom). This suggests that, even though VP/VS and VS are not known in detail, it is possible to improve the performance of acoustic FWI of elastic (pressure or multicomponent) data if a reasonably smooth guess of the ratio at the sea bottom is used with the suggested workflow.

Impact of noise

We now study the effect of noise on the results. First, we add white Gaussian noise within the bandwidth of the source to the elastic true data. Figure 11a shows a representative shot gather of the elastic true data with added noise, whereas Figure 11c shows the difference with the true acoustic data in Figure 7a. We invert this data set acoustically keeping the same acquisition geometry and inversion parameters as in the previous tests — which includes a mild smoothing to reduce the ill posedness of the problem — and the result is shown in Figure 12a. The recovered P-wave velocity model is less resolved than that obtained from the noise-free data (Figure 7d), and the noise introduces some artifacts in the recovered P-wave model because the predicted data cannot match the observed data due to the elastic effects but also the added noise.

The workflow is applied again with this recovered P-wave model (step 1 in the workflow in Figure 2). Forward acoustic and elastic data are modeled (step 2) using the latter, and an S-wave model built from the VP/VS model illustrated in Figure 9b, which leads to the P-wave velocity model in Figure 12b after two iterations of the suggested workflow (steps 3−5). Figure 11b shows the corresponding matched elastic data set, in which an automatic mute has been applied prior to the data matching operation to reduce the length of the filters and hence the computation cost, and its difference with the true acoustic data in Figure 7a is shown in Figure 7d. The recovered P-wave velocity model in Figure 12b is not as close to the true model as that obtained from the noise-free data (Figure 7d), but it is better resolved and more accurate than the P-wave velocity model after acoustic FWI of the noisy elastic data (Figure 12a). We note, however, that this improvement is small when compared with the noise-free test: There is an rms change of nearly 8% when comparing Figure 8d and 8e, whereas this is of 1.4% for Figure 12a and 12b. Thus, noise introduces some artifacts in the recovered models, which limits the performance of the current method, but the recovered P-wave models are slightly more accurate and well-resolved than those after conventional FWI of noisy elastic data.

Comparison with a full elastic FWI result

We have previously shown that the quality of the recovered P-wave velocity model after acoustic FWI of the elastic data can be improved by using the suggested method to mitigate elastic effects and an average guess of VP/VS. We now compare this result with that after elastic FWI of the elastic data. The P- and S-wave slownesses are simultaneously inverted using the same acquisition geometry and range of inverted frequencies as in the previous examples and a forward modeling code that solves the isotropic heterogeneous elastic wave equation, which is used to form the kernel of the FWI code (Guasch et al., 2012). The starting P-wave velocity model in Figure 8b and a starting S-wave velocity model generated from the latter and the average VP/VS guess in Figure 9c are used as starting models. Figure 13a and 13b shows the recovered P- and S-wave velocity models, respectively.

We observe that the recovered P-wave velocity model in Figure 13a is more resolved and accurate than that after acoustic inversion of the elastic data in Figure 8d because elastic effects are accounted for. However, the result is not as good as that after acoustic FWI of acoustic data (Figure 8c) because the recovered P-wave velocity model is affected by the quality of the recovered S-wave model in Figure 13b. This is due to crosstalk between the inverted parameters. Other multiparameter implementations might lead to improved models, but this is not discussed here. When compared with the recovered P-wave velocity model after applying the current method using an average VP/VS model (Figure 8f), the elastic inversion leads to a slightly better result, but the difference in terms of accuracy and resolution is not significant. Hence, the current method performs similarly to elastic FWI of elastic data in this case without the need of performing full elastic FWI, i.e., with a considerable reduction in computation time.

A land synthetic data set

We now apply the suggested workflow to mitigate elastic effects on a synthetic land data set based on the marine 2D SEG/EAGE overthrust model (Aminzadeh et al., 1997). We remove the top water layer from the latter, so it is effectively a land model. The true S-wave model is built by using a constant Poisson’s ratio of 0.24 (Brossier et al., 2009), i.e., a constant velocity ratio of VP/VS=1.71. We generate acoustic and elastic single-component pressure data from these models using a grid spacing of 10 m, which ensures that there is no dispersion of either P- or S-waves less than 17.5 Hz. A dense array of 961 receivers at 15 m depth with a receiver spacing of 10 m and a total number of 95 sources at 5 m depth uniformly distributed across the model with a source spacing of 100 m are used, in which the wavelet has a useful bandwidth that spans from 2 to 17 Hz. As in a realistic field data set, the sources are located in a layer with nonzero S-wave velocities and, hence, P- and S-waves are generated close to the sources. Absorbing boundaries are used to mitigate spurious reflections at the edges of the model.

Figure 14a and 14b shows the true acoustic and elastic data for a representative shot gather at the center of the model generated using the true P- and S-wave models. Their difference with respect to the acoustic modeled data is shown in Figure 14d and 14e, respectively. We observe a difference in amplitude of the P-wave arrivals, including P-wave refractions and reflections and the direct P-wave. This is due to S-waves also being generated close to the sources: The wavefield is already split into P- and S-waves already within the first layer. Additionally, we observe events in Figure 14b that correspond to S-waves and converted P- to S-waves, for example, there are refractions with different slopes that are not visible in the acoustic data.

Figure 15 shows the improvement in the recovered P-wave velocity when applying our workflow. Figure 15a is the true P-wave velocity model, and Figure 15b is the starting model used for the inversions, which has been obtained by smoothing the former (in slowness) 180 times with a Gaussian function of 30 m correlation length. Acoustic FWI is initially performed on the acoustic data, and the recovered P-wave velocity model is shown in Figure 15c. For this data set, all inversions are performed using an amplitude-balanced objective function during acoustic FWI, as described by Shen (2010). The different layers are well-recovered after inverting the data from 3.5 to 16.7 Hz with a total of 144 iterations. Figure 15d shows the result of performing acoustic FWI on the elastic data, which leads to a P-wave velocity model with poorer resolution. Note, for instance, the loss in resolution of the layers in the overthrust region when compared with Figure 15c, the higher average model rms difference within the dotted area, and how the low-velocity channel at approximately 800 m depth on the right side of the thrust is incorrectly pushed down due to the slow S-wave velocities at the top of the model. The use of an amplitude-balanced objective function to mitigate elastic effects is therefore not sufficient here because the elastic effects are too strong. We now combine the recovered result using this method and our suggested workflow.

We now generate acoustic and elastic modeled data using the recovered model in Figure 15d and an S-wave velocity model built by using the true VP/VS model (step 2). Matching filters that map the elastic to the acoustic modeled data are computed for nine traces simultaneously, and these filters are then convolved with the elastic true data (e.g., Figure 14b) to reduce the elastic effects from the data (step 3).

Acoustic FWI is then performed on the matched data (step 4), and the workflow is repeated once more to further mitigate the elastic effects (step 5). Figure 14c and 14e shows the resulting matched data and its difference with the true acoustic data, respectively. The reduced amplitude of all arrivals in Figure 14e compared with Figure 14d demonstrates that the elastic effects have been partially removed.

The result of applying acoustic FWI on the matched data set is shown in Figure 15e. The recovered P-wave velocity model shows the benefits of using our workflow to mitigate elastic effects, even when compared with the result of not addressing elastic effects and using an amplitude balanced objective function (Figure 15d): The different layers are significantly better resolved than those in the recovered model after acoustic FWI of elastic data in Figure 15d (as indicated by the rms values on the figure), and the low-velocity channel at the top of the model is now in the correct position.

To further assess this result, the rms error of the inverted model using the true model as the reference is computed as a function of depth and position within the dotted box in Figure 15 for the different recovered models and the starting model in Figure 15b, and it is presented in Figure 16. Again, this shows the benefit of using our approach because the smaller this value is, the closer it is to the true model. Note that the curve corresponding to the acoustic FWI of the matched elastic data (the gray line with squares) is closer to the true model than the curve corresponding to the acoustic FWI of the elastic data (the thick black line with circles). As in the previous example, the recovered model using the matched data is not better than the recovered model after acoustic FWI of acoustic data (the black line).

Finally, we compute the average model rms difference per iteration with respect to the true model within the dotted area in the models in Figure 15, and the result is presented in Figure 17. We make three observations: (1) The model rms decreases with increasing frequencies (iterations). (2) This decrease is less pronounced for the acoustic FWI of elastic data (the thick dark-gray line with circles), especially at high frequencies. (3) The curve corresponding to the acoustic FWI of the matched data (the light-gray line with squares) is closer to the true model than that of the acoustic FWI of elastic data (the dark-gray line with circles) and is nearly always very close and following the same trend as the acoustic FWI of the acoustic data (the black line).

These results suggest that the application of our workflow to a land synthetic data set has succeeded in mitigating the most important elastic effects, and it has provided a recovered P-wave velocity model that is closer to the true model and is better resolved than the model obtained using conventional acoustic FWI.

DISCUSSION

We have shown that the application of our workflow successfully addresses the elastic effects in acoustic FWI of marine and land synthetic data sets and it performs better than standard methods when elastic effects are strong. To reduce unwanted effects on the matched data due to the finite length of the filters and the fact that we match data with different spectra, we have computed Wiener filters that are optimum for more than one trace, thus smoothing the matched data laterally. We have chosen a different number of traces to be matched in each data set presented here because the grid samplings are different. However, we have observed that a reliable result can be obtained if we simultaneously match traces with spacing of at least 150 m between the first and last traces and that a more than 250 m separation degrades the result because the resulting matched data are too strongly smoothed. In terms of the length of the filters, we have observed that a length between 300 and 350 ms provides the best results in terms of continuity and quality of the matched data as well as an affordable computation cost of the matching filters.

We have shown with tests on 2D horizontally layered media that the current workflow does not provide accurate results when the starting model is poor at the sea bottom and elastic effects are strong because acoustic FWI of the observed elastic data provides inaccurate recovered P-wave velocity models to be used with the current workflow. We have suggested and implemented a modified workflow to successfully deal with these situations at a similar cost to a single acoustic inversion.

Our method has been applied to pressure data, but its application to multicomponent data is straightforward by modeling either pressure or the components of particle velocity during the acoustic and elastic modeling step (step 2). In terms of computational cost, the present method is two to three times more demanding than acoustic FWI per iteration of the workflow because each repetition requires two additional forward computations, a matching filter calculation and an acoustic inversion. In the tests shown here, one or two iterations of the workflow were enough to mitigate the most important elastic effects. Hence, it is an order of magnitude less demanding than elastic FWI when applied to 3D data sets (Guasch et al., 2010). To further reduce its computational cost, we suggest using a coarse grid for the acoustic inversion and the data matching steps and a fine spatial grid to model the acoustic and elastic data especially to avoid dispersion due to low S-wave velocities.

Future tests will focus on addressing other amplitude effects that are not typically accounted for in acoustic FWI. For instance, attenuation is generally disregarded when carrying out FWI. However, this can be addressed by estimating a model of attenuation and applying a procedure similar to the one outlined by Agudo et al. (2017). In fact, this strategy could be combined with the one used in this paper to mitigate viscoelastic effects all at once, thus taking into account attenuation and elasticity simultaneously.

Finally, we conclude that further research should also focus on exploring other approaches for matching the data and estimating VP/VS prior to the application of our workflow using global methods. The latter would be appropriate given that a smooth average guess of VP/VS is sufficient to mitigate most of the elastic effects and, hence, only a few parameters would be needed in a global inversion.

CONCLUSION

Acoustic FWI is typically performed to obtain P-wave velocity models of the subsurface, neglecting elastic effects, although these are generally present in the data. Here, we have suggested and implemented a method to mitigate elastic effects in acoustic FWI of pressure or multicomponent data based on matching filters and a guess of VP/VS, which is two to three times more expensive than acoustic FWI, but it is much less computationally demanding than elastic FWI. The application to marine and land synthetic data sets shows the improvement in terms of resolution and accuracy of the layers of the recovered P-wave velocity models with respect to those obtained after conventional acoustic FWI. These improvements are less significant in the marine environment because elasticity is not as strong in this case. However, we have shown that elasticity has a strong impact on the recovered P-wave velocity models of land data. The application of the current workflow to the land data set leads to significantly better P-wave velocity models, even at shallow depths. Additionally, tests on the accuracy of the VP/VS estimate show that in marine environments, only a smooth average guess of the ratio that contains a more accurate value at the vicinity of the seabed is required for the workflow to be successful. Future tests will focus on applying this method on different field data sets, which inherently contain elastic effects. The method can also be further extended to account for other amplitude effects that are typically not fully accounted for in acoustic FWI.

ACKNOWLEDGMENTS

We appreciate the comments from associate editor A. Guitton, reviewers G. Alves and J. Hobro, and three anonymous reviewers who have helped to improve the quality of this manuscript considerably. Imperial College London also wishes to thank the sponsors of the FULLWAVE consortium.

Freely available online through the SEG open-access option.