Localized time-lapse elastic waveform inversion using wavefield injection and extrapolation: 2-D parametric studies

We present a methodology to invert seismic data for a localized area by combining source-side waveﬁeld injection and receiver-side extrapolation method. Despite the high resolving power of seismic full waveform inversion, the computational cost for practical scale elastic or viscoelastic waveform inversion remains a heavy burden. This can be much more severe for time-lapse surveys, which require real-time seismic imaging on a daily or weekly basis. Besides, changes of the structure during time-lapse surveys are likely to occur in a small area rather than the whole region of seismic experiments, such as oil and gas reservoir or CO 2 injection wells. We thus propose an approach that allows to image effectively and quantitatively the localized structure changes far deep from both source and receiver arrays. In our method, we perform both forward and back propagation only inside the target region. First, we look for the equivalent source expression enclosing the region of interest by using the waveﬁeld injection method. Second, we extrapolate waveﬁeld from physical receivers located near the Earth’s surface or on the ocean bottom to an array of virtual receivers in the subsurface by using correlation-type representation theorem. In this study, we present various 2-D elastic numerical examples of the proposed method and quantitatively evaluate errors in obtained models, in comparison to those of conventional full-model inversions. The results show that the proposed localized waveform inversion is not only efﬁcient and robust but also accurate even under the existence of errors in both initial models and observed data


I N T RO D U C T I O N
Tomographic studies seek best models, of which synthetic data are in concordance with observed data.In general, seismic data are recorded on the Earth's surface.Amongst many proposed tomography methodologies in seismology, full waveform inversion has emerged as one of the most promising and powerful techniques since over three decades ago (Tarantola 1984;Woodhouse & Dziewonski 1984;Mora 1988;Pratt 1999;Virieux & Operto 2009;French & Romanowicz 2014).However, waveform inversion requires full wavefield modelling.It is thus computationally expensive to invert for a large area.
Recently, there are several attempts to focus wave-equation based seismic imaging methods on a small target area in the Earth rather than the entire region.Wang et al. (2016) iteratively inverted for a small region at the uppermost mantle in the vicinity of receiver arrays using teleseismic events and wavefield injection method for the source side.For the receiver side, wave extrapolation can mitigate the computational burden of waveform tomography.Dong et al. (2009) shows after wavefield extrapolation, the computational cost of reverse-time migration in a target region is significantly reduced.Yang et al. (2012) shows localized inversion after doubleside wavefield extrapolation is computationally more efficient, and converges faster compared to a standard full-model waveform inversion.Haffinger et al. (2013) indicated that due to the reduced computational domain after wavefield extrapolation, higher frequency contents could be introduced into the waveform inversion, which naturally improves the image quality.Baring the previous studies cited above for local seismic imaging in mind, in this study, we will (i) propose a methodology that allows to locally invert seismic waveform data for the region disconnected from both source and receiver arrays, combining source-side wavefield injection method and receiver-side residual wavefield extrapolation; and (ii) examine several 2-D synthetic tests to see the efficiency, robustness, accuracy and limitations of the proposed method.
During the oil and gas production in the subsurface, it is crucial to quantitatively estimate the model changes in the reservoir fluids.Reservoir monitoring requires a huge investment.In order to maximize the return on investment, one needs to take full advantage of a large volume of time-lapse seismic data.The method that we are proposing can be particularly appealing for time-lapse surveys.We have the 'baseline' data set that is acquired before the production and the 'monitor' data set that is acquired on a regular basis (daily or weekly) during or after production.We are then to interpret data difference for model difference.Unfortunately, full-model waveform inversion for every monitor data set would be still too expensive.Since the production does not change the overall substructures, there should be much redundancy in conventional full-model inversions.Furthermore, the structural changes will mainly occur in the region of fluid removal or replacement, which are far from source and receiver arrays.As the model alterations are small and deep, the corresponding changes in seismic data go only a few per cents above the background noise.For the reasons outlined above, it is indispensable to develop an efficient and accurate waveform-based method to image the local structural changes.
Fig. 1 shows how our methodology works.We have a source (array) and a receiver array.Instead of performing inversion inside the whole model region, we would like to focus on a subvolume in V + , where the model changes are expected to occur.All forward and backward simulations will be conducted within the reduced absorbing boundary (the dashed square shown in Fig. 1).Such a strategy of localized waveform inversion requires two levels of wavefield reconstructions: (i) representation of seismic sources in the subdomain (wavefield injection method) and (ii) extrapolation of residual wavefield between baseline and monitor surveys recorded at the physical receivers to virtual receivers in the subdomain (residual wavefield extrapolation).Borisov et al. (2015) proposed a localized waveform inversion technique for time-lapse surveys within a target region disconnected from source arrays.However, there are few studies that extended it to a region disconnected both from source and receiver arrays.As we will show in the later sections, receiverside wavefield extrapolation combining with source-side wavefield injection (Robertsson & Chapman 2000;Capdeville et al. 2003;Singh & Royle 2010;Monteiller et al. 2013;Masson et al. 2014;Borisov et al. 2015) can naturally filter out the non-physical extrapolated events caused by limited receiver aperture (the single-sided illumination problem).We formalize localized waveform inversion problem for time-lapse seismics and present the two levels of wavefield reconstruction methodology.To demonstrate the feasibility and robustness of our methodology, we perform a series of synthetic tests with noise and uncertainties in both model and data domain.

L O C A L I Z E D WAV E F O R M I N V E R S I O N F O R T I M E -L A P S E S U RV E Y
In this section, we formulate localized waveform inversion for timelapse surveys.We first start with the classical L 2 norm misfit function of waveforms under the framework of Tarantola (1984) and Mora (1987): with d and u(m) being the gathered waveform data points of observed and synthetic computed for a model m, respectively.By evaluating the first variation of the misfit function, we obtain a Newton method expression: and where which is partial derivatives (or Fréchet derivatives) of the wavefield with respect to model parameters and A T A is an approximated Hessian matrix.The right-hand side of eq. ( 2) is the gradient direction ∂ E(m)/∂m that can be calculated by either cross-correlation of partial derivatives and waveform difference between observed and synthetics or back propagation of waveform difference in the media m i .Multiple iterations (eq. 3) update model to overcome the nonlinearity of this inverse problem.The details of waveform inversion can be found in many literatures, some of which are cited herein.
We use above formulation to obtain the 'baseline' full model m base before the production mode.We then obtain 'monitor' data with some deformation in the subsurface media and we would like to invert for the model alteration from the 'baseline' model obtained during the first step.There are mainly three strategies to monitor the changes of target regions through inversions (Asnaashari et al. 2015): (i) 'parallel difference strategy' independently inverts time-lapse data sets, that is, the structure changes will be calculated through the subtraction of inverted baseline model from inverted monitor models; (ii) 'sequential difference strategy' will use the most recently acquired former monitor model as an initial model for inversions; (iii) 'double-difference strategy' inverts only time-lapse data difference for a model change instead of full-data inversions.Double-difference is widely applied in earthquake location and seismic imaging (Waldhauser & Ellsworth 2000;Zhang & Thurber 2003;Zheng et al. 2011;Zhang & Huang 2013;Asnaashari et al. 2015;Yang et al. 2015) and proves to be more robust and accurate in retrieving time-lapse variations and suppressing artefacts of inversions.Double-difference strategy requires strict repeatability of seismic surveys, which is generally satisfied by the state-of-the-art acquisitions, especially in time-lapse ocean-bottom seismics (Granger et al. 2005;Ronen et al. 2005).We thus seek to minimize the difference of observed and synthetic waveform Downloaded from https://academic.oup.com/gji/article/209/3/1699/3077207 by CNRS -ISTO user on 10 August 2022 residuals between baseline and monitor surveys: with superscripts 'base' and 'mon' denoting baseline and monitor, respectively.We define the waveform residual between baseline and monitor data as To minimize the new cost function defined above, we assume that the baseline model is close to the real structure and that the residuals between observed and synthetic waveforms are of the same order for baseline and monitor surveys.The Newton method for this new cost function is Note that A does not change from that of eq. ( 4).Since we are interested in a model perturbation δm mon happening inside the region far deep from the source and receiver arrays, we reconstruct d and u inside the small region of interest by wavefield injection method and residual wavefield extrapolation method that we discuss in the following sections.
Forward modelling is performed using 2-D staggered grid elastic isotropic finite-difference method with fourth-order accuracy in space and second-order accuracy in time (Virieux 1986;Levander 1988).We use the convolutional perfectly matched layers (Roden & Gedney 2000;Komatitsch & Martin 2007) to absorb undesirable seismic wave reflections from the model borders or outside the target regions.The gradient is obtained from adjointstate method (Liu & Tromp 2006;Plessix 2006) and the conjugate gradient method is used to approximate the Hessian matrix (Polak 1971).Compressional wave speed (V p ) and shear wave speed (V s ) are independently inverted, while V s (V s ≈ V p /1.67) and the density model is linearly linked to V p (Gardner et al. 1974).condition during wavefield injection.Borisov et al. (2015) applied it to waveform inversion for 3-D media with physical receivers inside the region.
Based on the injection boundary, the modelled region is divided into two subdomains, V + and V − (Fig. 1).The wavefield in either subdomain satisfies the same elastodynamic equation and can be expressed as where u denotes the displacement wavefield.The superscript + and − denote the wavefield inside and outside the injection boundary, respectively.The subscript c represents the complete wavefield and s represents the scattered wavefield caused by the model changes within injection boundary during time-lapse surveys.u b denotes the displacement wavefield excited by the stored boundary wavefield.The stress wavefield along the injection boundary satisfies the same condition as the displacement wavefield.
Even though u + in V + and u − in V − are achieved simultaneously, they are decoupled due to the boundary condition applied on injection boundary.Specifically, we add the stored boundary data for the wavefield inside the injection boundary, while we subtract the boundary data for that outside injection boundary.Since the physical sources are extrapolated from the surface to the source injection boundary during wavefield injection, we can naturally decrease the size of absorbing boundary (the dashed square in Fig. 1) in order to reduce the inversion area.This may prohibit the interactions between the scattered wavefield and the structures outside the absorbing boundary.However, the missing part of the wavefield is only second or higher-order scattered wavefield, whose amplitude is weaker and arrivals are later than that of the first interactions.Also, it is numerically expedient considering the associated computational cost by applying full-model absorbing boundary or the so-called exact boundary condition (van Manen et al. 2007;Willemsen et al. 2016).
In order to better illustrate the role that wavefield injection plays as the localized modelling engine during time-lapse surveys, here we perform a synthetic experiment and compare it with the fullmodel simulations (Fig. 2).The 2-D baseline (a perturbation-free model within the injection boundary) is 8 km in the x-direction and 4 km in the z-direction, while the reduced modelling region is 2.8 km in the x-direction and 1.4 km in the z-direction.To store the wavefield along injection boundary, we are supposed to perform the full-model forward modelling once at the beginning for each physical source.During time-lapse surveys, we assume that perturbations occur only inside the injection boundary (Figs 2e and i).h).However, at later times, the absorbing boundary in localized forward modelling will prevent the interactions between scattered wavefield and the medium outside the target region.Therefore, neglecting the second-or higher-order wavefield when applying wavefield injection method will cause small errors at later time steps, as observed in Fig. 2(l).However, we will show in Sections 5 and 6 that this imposes limited effect on the final inversion results through various synthetic experiments.

R E C E I V E R -S I D E R E S I D UA L WAV E F I E L D E X T R A P O L AT I O N
Extrapolation of wavefield recorded at a physical (dense) receiver array is needed in our methods in order to perform localized waveform inversion.In this study, it is identical to what is applied commonly in reverse time migration (Baysal et al. 1983;Schuster 2002).We extrapolate residuals between baseline and monitor with a truncated surface integral due to the receiver array geometry in frequency domain (Oristaglio 1989;Wapenaar & Fokkema 2006): where S denotes the surface covered by the physical receiver array, which cannot always be a closed boundary (cf.Fig. A1).u 1 n (r) is the ith component of perturbed displacement at r, n j denotes the normal vector, G 0 in (r, r A ) denotes the ith component of Green function from nth single force at r A to r for baseline model, whereas C 0 i jkl are elastic moduli for the baseline model.The derivation of the equation above can be found in Appendix A.
This leads to the so-called single-sided illumination problems which are commonly found in migration and seismic interferometry community (Snieder et al. 2006;Wapenaar & Fokkema 2006;Löer et al. 2014).The associated spurious events are particularly severe when the subsurface is strongly inhomogeneous and receiver apertures are very limited.In this study, we choose to apply the conventional extrapolation approach.Other extrapolation techniques, such as redatuming through multidimensional deconvolution (van der Neut et al. 2011;Wapenaar et al. 2011) 11).Panels (g) and (h) correspond to the 1-D velocity profile along the dashed lines in panels (e) and (f), respectively.Black, cyan and red curves denote the true, starting and inverted velocity respectively.For simplicity, we use the same notification in the following experiments and we will not mention them repeatedly.(Wapenaar et al. 2014;Meles et al. 2016), will be interesting to test against our methods in the future work.
The Green's functions G 0 * in (r, r A ) between physical and virtual receivers, serving as back-propagator in eq. ( 10) are calculated and stored just once using the baseline model obtained by the first full-model waveform inversion.We generate them by posing single source forces at r A .The fact that we truncate the surface integral to S means that we neglect the interactions between perturbed backpropagators (caused by time-lapse model variations) and the recorded wavefield (Ravasi & Curtis 2013).The evaluation of the integral in eq. ( 10) requires observed data of both displacement-(or particle velocity-) and stress-/strain-types.This is practically difficult since only displacement components (or velocity, acceleration) are recorded in most current land surveys but not stress nor strains.Even though replacing stress tensor by particle velocity can still guarantee the extrapolated wavefield being kinematically correct, it would result in additional spurious events.However, in marine surveys, stress tensor can be recorded by ocean-bottom 4C acquisitions (3C displacement/velocity plus 1C fluid pressure) since normal traction will be simplified to pressure due to fluid-solid boundary condition (Ravasi & Curtis 2013).Moreover, ocean-bottom seismics have additional advantages compared to conventional marine streamer acquisitions.Receivers fixed to the ocean bottom will not only improve repeatability for time-lapse surveys but also provide shear wave information, which leads to higher resolution due to shorter wavelength.Therefore, ocean-bottom 4-D acquisitions are particularly suitable for time-lapse monitoring and elastic wavefield extrapolation.

Wavefield extrapolation -2-D numerical examples
The proposed extrapolation method was tested on a 2-D elastic layered model (P-wave velocity shown in Fig. 3a).Density is linked to V p through Gardner's rule (Gardner et al. 1974).In Fig. 3(a), grey triangles and blue triangles represent physical receivers and extrapolated receivers, respectively.Extrapolated depth is 1 km.The time-lapse data difference is caused by the ellipsoid perturbation enclosed in the solid black square (Fig. 3a).Green's functions are obtained from the baseline model without the ellipsoid perturbation.Comparisons of observed (black curves) and extrapolated (red waveform) wavefield of selected two shots are shown in Figs 3(c) and (d).We can see that the extrapolated wavefield is generally accurate both in phase and amplitude.Some spurious events come out especially when shot position is at the edge of receiver plane due to loss of stationary phase points (Snieder et al. 2006).

L O C A L I Z E D WAV E F O R M I N V E R S I O N S -2 -D N U M E R I C A L E X A M P L E S W I T H L AY E R E D M O D E L S
In this section, we present numerical examples in order to examine the efficiency, robustness and accuracy of the proposed localized waveform inversion.Conventional full-model inversion is tested against our method.The grid spacing for wave propagation is 20 m   and d) is 8 km in width and 4 km in depth.The time-lapse model has a lens-shaped velocity perturbation inside V + (Figs 4a and b), which is around 6 per cent of the baseline velocity for both V p and V s .
In this section, we assume that the baseline model is identical to the 'true' baseline model except Sections 5.5 and 5.6, where we use the unknowingly perturbed and inverted baseline model, respectively.We define the model error associated with the inverted time-lapse perturbation within the target region as Model error (in per cent) we would like the reader to look the relative differences in Table 1 rather than to compare the face values with other inversion publications.The comparisons of each experiment are summarized in Table .1, where the test number corresponds to the number included in the title of following subsections.Note that all the full-model and localized inversion results are obtained after the same iterations.

Full-model inversion (#1)
In the first experiment, we show the conventional full-model waveform results starting from the true V p and V s baseline model (Figs 4c and d).We invert for the full time-lapse model shown in Figs 5(a) and (b), with the zoom-in of the target variation (Figs 5e and f).Besides the true perturbation (Figs 5c and d), we also find some inversion artefacts outside the target region as we can see in the 1-D velocity profile (Figs 5g and h).Due to the high ill-posedness and nonlinearity of the inverse problems, it can hardly ensure that the model update is focused on the target time-lapse perturbation.
The more unfavourable fact is that it is unavoidable to calculate the wavefield for the whole model every iteration, although the time-lapse perturbation occurs only in a certain small target region.Therefore, considering the computational cost, it is nearly impossible to conduct full-model waveform inversion on a daily or weekly basis for practical scale problems.

Localized inversion (#2)
The injection boundary data and the Green's functions are stored beforehand through a series of full-model forward modelling.With the source-and receiver-side extrapolation as we explained in Sections 3 and 4, the inversion area could be naturally reduced to a smaller region (the dashed square shown in Figs 4a and b) where time-lapse perturbation is expected to occur.Localized inversion provide us with better results for both V p and V s (Figs 6a and b) compared to the full-model inversion (Figs 5e and f).Since model update is solely allowed within the target region, the potential inversion artefacts can be mostly prevented.This is equivalent to a preconditioning of inverse problems which leads to a more accurate gradient direction.Also, localized inversion results look sharper compared to the full-model inversion due to the reduced offset and increased moveout, thus we obtain higher sensitivity to the smallscale perturbations.More importantly, the computing time in this test is reduced by a factor of 10 compared with the full-model inversion.The achieved computational efficiency is actually proportional to the size of the non-calculated area outside the target regions.This shows the feasibility of time-lapse elastic full waveform inversion in an efficient way for practical scale problems.

Single-shot inversion (#3)
In conventional waveform inversions, long-offset data sets are crucial to retrieve the long wavelength information and mitigate cycle skipping problems (Shipp & Singh 2002).However, during timelapse surveys, we expect that the long-wavelength structures (baseline model) remain almost unchanged and then the importance of long-offset can be reduced (Plessix et al. 2010).To invert smallscale velocity perturbations occurred locally in a target region, there is probably a redundancy when applying the same acquisitions for baseline surveys.To show the capability of the proposed localized inversion in reducing field cost (the number of shots and receivers), we perform a single shot inversion for both full-model and localized inversion.The results show that the full-model inversion (Figs 7a-d

Localized inversion with white noise (#4)
In this test, we would like to show the robustness and accuracy of the localized inversion in presence of random noise.Artificial Gaussian white noise is added to the recorded data set before wavefield extrapolation.The signal-to-noise ratio of extrapolated wavefield is around 6 dB (black curves in Fig. 9).Even with strong white noise, the localized inversion is still quite robust and accurate as it shows in Fig. 8.

Localized inversion with unknown baseline perturbation (#5)
In previous tests, there is an assumption that the baseline model (especially the parts outside the target regions) is unchanged during the time-lapse surveys.However, it is possible that small perturbation occurs during drilling or oil production.In order to test the robustness of the localized inversion in the case of slight model changes outside the target region, we introduce an unknown perturbation outside the target zone during time-lapse surveys (Figs 10a and b).Note that wavefield injection and extrapolation are still performed using the original baseline model, which does not include the shallow perturbation.The inverted perturbations are still quite accurate as is shown in Figs 10(c)-(f).
Due to the inaccurate extrapolation propagators which do not contain the effect of the unknown shallow perturbation, some nonphysical events will appear in the extrapolated wavefield (Fig. 11).The non-physical events result from the wrongly extrapolated shallow reflections caused by the unknown perturbations and come out before the first arrival at the virtual receivers.We can notice that the localized inversion does not fit the spurious events through comparing modelled wavefield and extrapolated wavefield (black curves in Fig. 11).This is because wavefield injection naturally serves as anti-causal filters and exclude the effect of the short wavelength extrapolation error existing before the physical first arrivals.In other words, the cross-correlations between the spurious extrapolated events and partial derivatives are quasi-orthogonal (eq.2).It is interesting to refer to Kawai et al. (2014) who performed 3-D localized waveform inversion for the base of the mantle using sources and receivers disconnected from the Earth's surface.They  also studied the risk of over-interpretation of data for the structure outside the inversion box but it turned out that the error is small due to quasi-orthogonality of partial derivatives for points inside and outside the box.
However, in case of severe perturbations, for instance, the associated first arrivals are apparently different from those of the unperturbed baseline model, we are supposed to do the full-model inversion once again to make sure the long wavelength information is accurate enough for wavefield injection and extrapolation.We might be able to control this by looking at the kinematics of timelapse data difference.We can qualitatively conclude that slightly perturbed baseline model will not seriously affect the accuracy of the proposed localized inversion.In fact, as long as the baseline perturbations do not dramatically affect the first arrivals, the localized inversion will be still reliable and stable.

Localized inversion using inverted baseline model (#6)
To acknowledge the fact that in practice, no precise information on the interfaces and medium parameters is available, the boundary data for wavefield injection and the Green's functions for wavefield extrapolation should be modelled in a more realistic baseline model.Therefore, we divide inversions into two parts in this test: (i) baseline model inversion (ii) time-lapse model inversion.Due to the lack of illumination considering the sensitivity, we need to extend the survey area to retrieve the desired baseline model (Figs 12a and  b).Starting from a smoothed model (Figs 12c and d), the baseline models of V p and V s are inverted (Figs 12e and f) with associated 1-D velocity profiles shown in Figs 13(a)-(d).
With inverted baseline models, we repeat the same procedures as in previous localized inversion experiments.The inverted local perturbation (Figs 14a-d) seems still quite accurate as compared to the localized inversion using true baseline model in Section 5.2.

L O C A L I Z E D WAV E F O R M I N V E R S I O N S -2 -D N U M E R I C A L E X A M P L E S W I T H M O D I F I E D M A R M O U S I M O D E L
We further validate our localized inversion strategy on the modified Marmousi model (Institut Franais du Pétrole (IFP), 1988) in comparison to full-model inversions.The V p model is 7.5 km × 2.5 km (Fig. 15a).V s (V s ≈ V p /1.5 + 300 m s −1 ) and density models are linearly linked to the V p model (Gardner et al. 1974).We generate 4 s of data recorded at 300 m depth with a regular spacing of 25 m.In total 271 receivers record two components of velocity wavefield plus stress wavefield.A Ricker wavelet with dominant frequency of 6.5 Hz is used as the source wavelet and added onto the vertical velocity component.The grid spacing for the wave propagation is 25 m and the time step is 2 ms.55 sources are deployed at 275 m depth with the interval of 125 m.The time-lapse velocity perturbations are about 5 per cent of the average baseline model (Fig. 15a).
The 2-D V p perturbations obtained from full-model and localized inversions under true baseline model are shown in Figs 15(c (d), respectively.Compared to the conventional full-model inversion, the localized inversion retrieves the velocity anomalies with a higher accuracy after the same 100 iterations.The corresponding 1-D velocity profiles are shown in Fig. 16.Fig. 17 shows the vertical velocity components of initial and final residuals from one shot gather.The residuals of the major events are better fitted after localized inversion (Fig. 17d) in comparison to the final full-model inversion residuals (Fig. 17b).We specially plot the theoretical initial residuals (Fig. 17e) in case that we directly record wavefield at the same positions as the extrapolated virtual receivers.The residual difference between Figs 17(c) and (e) is shown in Fig. 17(f).This difference essentially represents the spurious events after wavefield extrapolation.However, we notice that the localized inversion does not 'see' those spurious events when comparing it with final residuals of the localized inversion (Fig. 17d).
We finally apply the same strategy on a smoothed Marmousi model without changing the geometry and inversion parameters (Fig. 18).Besides the reduced computational cost, the localized in-version is also more accurate than the full-model inversion (Fig. 19) after the same 100 iterations.

D I S C U S S I O N
In order to develop a localized waveform inversion scheme, we combine source-side wavefield injection and receiver-side extrapolation methods, where we pose some simplifications and assumptions.The closed surface integral over the source-side wavefield injection can be conducted in a rigorous numerical way, while receiver-side residual wavefield extrapolation suffers from the truncation of surface integral and sparseness of receivers.Also, it requires a challenging technology that efficiently measures strains or stresses if we would like to rigorously complete all the terms in eq. ( 10).During the source-side wavefield injection method, we ignored higher interactions due to local model changes between target region and the model outside the reduced absorbing boundaries.However, it turns out that it does not seriously affect the seismic images since Downloaded from https://academic.oup.com/gji/article/209/3/1699/3077207 by CNRS -ISTO user on 10 August 2022    higher-order interactions are weaker and only influence later phases than that we use for inversions.Receiver-side residual wavefield extrapolation is, in addition to the difficulty described above, assuming no change (or error) in baseline model by definition.But it shows that long-wavelength structures of the outer region is sufficient to image the interior of target regions.
We show that wavefield injection is efficient and robust enough even though part of the higher order scatter wavefield is neglected due to the reduced modelling area.With the extrapolated wavefield, albeit some spurious events due to one-side illumination problem or slightly perturbed baseline model, the localized inversion is still able to provide sufficiently accurate time-lapse images.The localized inversion, which naturally focuses on the scattered wavefield due to perturbations within target regions, will not only contribute to the reduction of computational cost to a great extent, but also to the suppression of spurious events that cause wrong model update outside target zones.We further calculated the model errors as a function of iteration indice in Fig. 20, showing that localized inversions converge faster in both cases (layered and Marmousimodels).Hence, the localized inversions require fewer iterations to obtain models with the same accuracy as full-model inversions.This will further enhance the speed-up and efficiency.Although we fixed the number of iterations to be 100 for all the tests in this paper, the truncation criteria should be discussed in our future papers.The decreased computational areas will potentially improve the capability in inverting higher frequency seismic data, which will lead to a better resolution, even though it is not shown here.Based on the synthetic examples we show, our methodology is more accurate, efficient and robust for the same number of iterations.
In data domain, after wavefield extrapolation, without losing the field interaction, the overburden effect is mostly reduced so that the nonlinearity of waveform inversion and associated cycle-skipping problems are mitigated.This can partly explain why localized inversions converge faster and the final misfit is smaller compared with conventional full-model inversion.In model domain, the localized inversion is equivalent to applying a precondition on the gradient, which naturally prevents the spurious model update from the non-target regions.This strategy is specially effective if we can specify our target zones and assume that the baseline model has few perturbations.In case of large baseline changes, we are supposed to perform the full-model waveform inversion once more to achieve better baseline model for wavefield injection and extrapolation.Although we have focused on the residual waveform inversion, the method proposed here is equally validated for localized inversion of seismic data in general.

C O N C L U S I O N S A N D P E R S P E C T I V E S
We proposed and examined a methodology of inverting seismic waveforms for localized structures and show its feasibility through various synthetic experiments by adding errors in both model and data domain.We combine wavefield injection and extrapolation to perform the localized waveform inversion for time-lapse seismic surveys.It shows that the localized inversion can enhance the image quality of local time-lapse variation while reducing the computational cost to a large extent.The proposed strategy reveals that putting virtual sources and receivers in the vicinity of target region contribute to improving the accuracy and robustness of waveform inversion even under uncertainties of the outside area.It is as if we observed through a telescope to make the distant objects look closer and thus more clearly, even without knowing what is happening between them and us.
Due to the reduced modelling size, the computing time of 2-D localized inversion in the experiments above is decreased by a factor of 10 as compared to full-model inversions.We expect that the more severe computational cost of 3-D elastic or viscoelastic full waveform inversion during practical time-lapse surveys, could be potentially reduced to a more acceptable and economic level.Although the proposed methodology is examined on 2-D cases, the extension to general 3-D cases should be straightforward.

A P P E N D I X A : R E P R E S E N TAT I O N T H E O R E M F O R A S M A L L M O D E L C H A N G E
We review representation theorem for reconstruction of residual wavefield extrapolation.Let be V + the region where we have a model change.Note that this definition of V + is not always as same as that in Section 3. Receiver-side residual wavefield extrapolation is reconstructing the residual between baseline and monitor data ( d in eq. 6) at some point just outside V + from d recorded at the Earth's surface (Fig. A1).As the representation theorem writes in a compact form in frequency domain, we first recall the equation of motion for an (an)elastic medium in frequency domain where ω is angular frequency, ρ is the density, v i is the ith component of displacement, C ijkl is the elastic moduli tensor and has complex value when the medium is anelastic, g i is ith component of external body force distribution.Any internal or external boundary condition can be added in a form of surface integral to this equation but we omit the discussion on this topic in this paper.The comma A ,i denotes spatial partial derivative with respect to ith direction.The superscript A 0 denotes the seismic structure of unperturbed model.We can write the equation of motion for the perturbed model in the same manner: where ρ 1 and C 1 i jkl denotes the perturbation of seismic model with respect to the initial model, u 1 i is the perturbed wavefield due to the model change with a source distribution f i .A * denotes the Downloaded from https://academic.oup.com/gji/article/209/3/1699/3077207 by CNRS -ISTO user on 10 August 2022 Figure A1.Schematic illustration of the geometry of representation theorem discussed in Appendix A. We would like to extrapolate (residual) wavefield recorded at the surface S to the internal point r A by using only Green's functions inside V − .complex conjugate in frequency domain or time-reversal invariance in time-domain.Assuming the medium is lossless, we use complex conjugate type (or time reversal in time domain) equation of motion for eq.( A2) so that we can develop correlation type representation theorem for backward wavefield extrapolation.When needed, it is sufficient to take off all * to obtain the representation theorem of convolution type for forward wavefield extrapolation.
We then multiply eq. ( A1) by u * i = (u 0 * i + u 1 * i ) and eq.( A2) by v 0 i , respectively, taking the difference and taking a volume integral over an arbitrary V: This is the generalized reciprocity theorem that takes account for a model change.Unlike the classical reciprocity theorem, the volume integral over C ijkl does not have anymore symmetric form due to the model change.Now we set the source distribution as where δ in denotes Kronecker's delta and δ(•) denotes Dirac's delta function.Substituting eq.(A4) into eq.(A1), taking the volume integral over V, we obtain the equation of motion for Green's function:

Figure 1 .
Figure 1.A schematic illustration of the proposed localized inversion in a layered model.Red star denotes the physical source.Grey and blue triangles denote physical and virtual receivers, respectively.The dashed square represents the absorbing boundary and the solid square represents the source injection boundary.

Figure 2 .
Figure 2. Snapshots of vertical component velocity wavefield using full-model and localized forward simulation.Dashed square denotes absorbing boundary and solid square denotes wavefield injection boundary.Red star represents the physical source.(a-d) Simulations on an model with no perturbation inside the injection boundary.(e-l) A perturbed model.(b,f,j) Snapshots of full-model simulations.(c,g,k) Snapshots of localized forward simulations through wavefield injection.(d,h,l) Synthetic wavefield inside the red circle.Red curves and black curves denote the waveform obtained from full-model and localized simulations, respectively.

Figure 3 .
Figure 3. Wavefield extrapolation examples.(a,b) represent 2-D elastic layered V p model.Red stars denote source positions.Grey and blue triangles denote physical and extrapolated receivers, respectively.Panels (c) and (d) show the comparisons of observed (black curves) and extrapolated (red waveform) wavefield of selected two shots in (a) and (b), respectively.

Figure 4 .
Figure 4. V p and V s models for full-model and localized inversion tests.(a,b) True time-lapse V p and V s model.The dashed square encloses the target zone where we perform the localized inversion.(c,d) True baseline V p and V s models.

Table 1 .
Comparisons of V p and V s model errors for all the inversion experiments.At earlier stage, the scattered wavefield caused by the perturbations have not propagated long enough to interact with medium outside the target area.Thus the seismograms produced by full-model and localized simulations are exactly the same (Figs 2d and

Figure 5 .
Figure 5. Full-model time-lapse inversion results.(a,b) Inverted time-lapse V p and V s model.(c,d) True time-lapse perturbation in the target zone, corresponding to the difference between time-lapse and baseline model within the dashed square in panels (a) and (b).(e,f) Inverted time-lapse perturbation of V p and V s in the target zone.The error in the title is obtained though eq.(11).Panels (g) and (h) correspond to the 1-D velocity profile along the dashed lines in panels (e) and (f), respectively.Black, cyan and red curves denote the true, starting and inverted velocity respectively.For simplicity, we use the same notification in the following experiments and we will not mention them repeatedly.

Figure 6 .
Figure 6.Localized time-lapse inversion results.(a,b) Inverted time-lapse perturbation of V p and V s .Panels (c) and (d) correspond to the 1-D velocity profile along the dashed lines in panels (a) and (b), respectively.Panels (e) and (f) denote the initial and final residuals of localized inversion respectively.

Figure 7 .
Figure 7. Single-shot full-model and localized inversion results.(a-d) Inverted time-lapse perturbation and corresponding 1-D velocity profile of single-short full-model inversion.(e-h) Inverted time-lapse perturbation and corresponding 1-D velocity profile of single-shot localized inversion.

Figure 8 .
Figure 8. Localized inversion with white noise.(a-d) Inverted time-lapse model and corresponding 1-D velocity profile along the dashed lines in panels (a) and (b).

2 V×
|m mon true −m base true | + dV m mon true −m base true 2 100 per cent (11) with subscript 'inv' denotes inverted model, 'ini' the initial model and 'true' the true model.Note that this error definition is directly measuring the model perturbation and not model parameters themselves, resulting in the fact that the face values are deliberately large even for the best cases.There are still few localized waveform inversion studies where error is defined:

Figure 10 .
Figure 10.Localized inversion results with unknown baseline (a,b) Time-lapse V p and V s model with unknown perturbation at shallow depth.(c-f) Inverted time-lapse perturbation and corresponding 1-D velocity profile along the dashed lines in panels (c) and (d).

Figure 12 .
Figure 12.Extended baseline model inversion.(a,b) True baseline V p and V s model.(c,d) Starting V p and V s model obtained from the smoothed true baseline model.(e,f) Inverted baseline V p and V s model.

Figure 13 .
Figure 13.1-D V p and V s velocity profile along the dashed lines from inverted baseline model in Figs 12(e) and (f).

Figure 14 .
Figure 14.Localized inversion with inverted baseline model.(a,b) Inverted time-lapse perturbation of V p and V s .(c,d) corresponds to the 1-D velocity profile along the dashed lines in (a) and (b).Panels (e) and (f) denote the initial and final residuals of localized inversion respectively.

Figure 15 .
Figure 15.V p model and inverted results.(a) The baseline V p model selected from the true Marmousi model and the acquisition geometry.Red star and grey triangles denote physical source and receivers, respectively.White triangles denote the extrapolated virtual receivers.Two arrows point to the locations where time-lapse perturbation occurs.(b) The zooming true time-lapse V p perturbation within the dashed square in (a).(c) Full-model inversion result.(d) Localized inversion result.The two dashed lines (#1 and #2) in (c) and (d) indicate 1-D inverted velocity profiles and perturbations which are shown in Fig. 16.

Figure 16 .Figure 17 .
Figure 16.1-D inverted V p profiles and perturbations.The left two columns denote 1-D profiles from the full-model inversion result in Fig. 15(c) and the right two columns denote the 1-D profiles from the localized inversion result in Fig. 15(d).(a-d) 1-D profiles with true (black curves), starting (cyan curves) and inverted V p .(e-h) 1-D true (black curves) and inverted (red curves) V p perturbations.

Figure 18 .
Figure 18.(a) The smoothed baseline V p model with the same acquisition geometry and time-lapse perturbations shown in Figs 15(a) and (b).The two dashed lines (#1 and #2) indicate 1-D velocity profiles which are shown in (b) and (c).(b,c) profiles with true (black curves) and starting (cyan curves) velocities.

Figure 19 .
Figure 19.Full-model and localized inversion using the smoothed Marmousi baseline model.(a) The time-lapse V p perturbations from full-model inversion.The two dashed lines (#1 and #2) indicate 1-D inverted perturbations which are shown in (b) and (c).(b,c) 1-D true (black curves) and inverted (red curves) V p perturbations from full-model inversion.(d) The time-lapse V p perturbations from localized inversion.The two dashed lines (#1 and #2) indicate 1-D inverted perturbations which are shown in (e) and (f).(e,f) 1-D true (black curves) and inverted (red curves) V p perturbations from localized inversion.

Figure
Figure Comparisons of inverted errors as a function of iteration number between full-model and localized waveform inversions.(a,b) Using layered model (Sections 5.1 and 5.2).(c,d) Using modified Marmousi model (Section 6).