{'fraction_non_alphanumeric': 0.041744548286604365, 'fraction_numerical': 0.02367601246105919, 'mean_word_length': 4.4363533728478695, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Localized Wavefield Inversion (LWI): an Adaptation of Multi-Block ADMM for Localized FWI

Full-waveform inversion (FWI) is a high-resolution and computationally intensive imaging technique to reconstruct unknown parameters in the computational model in which the waves propagate; however, an accurate model of only part of this medium is required for some applications. To decrease the computational burden of such problems, target-oriented FWI was proposed where the redatumed data on the part of the medium or localized solvers for the wave equation are used. On the other hand, the classical formulation of FWI suffers from non-linearity and ill-posedness, which makes FWI sensitive to the initial model, the low-frequency content of the data, and limited illumination. In this study, we propose a localized version of the alternating direction method of multipliers (ADMM)-based FWI method, which was proposed to solve these problems in classical FWI. In our localized FWI or LWI, the medium is decomposed into a few subdomains, where some of them are updated, and the others are kept fixed based on an adaptation of multi-block ADMM, which is a powerful algorithm for solving inverse problems with decomposition and block separability. Numerical tests on the Marmousi model for a time-lapse application confirm the computational efficiency and robustness against background velocity model errors.


INTRODUCTION
Seismic full-waveform inversion (FWI) is an accurate imaging technique that can reconstruct the high-resolution elastic properties of the medium (Pratt et al., 1998). On the other hand, it suffers from non-linearity, ill-posedness, and a huge computational burden, all of which need to be addressed to make FWI a practical tool for imaging in different applications (Virieux and Operto, 2009). FWI generally estimates model parameters in the whole computational domain in which the seismic waves propagate, which is computationally expensive and not always necessary. For some specific target-oriented applications, often related to time-lapse imaging like the monitoring of reservoirs during hydrocarbon production or CO 2 injection (Willemsen, 2017), or other applications like salt boundary inversion (Lewis et al., 2012) or geothermal exploitation, one only seeks to update model parameters or their evolution over time within a localized region of interest. Such context is also encountered in earthquake seismology for lithospheric imaging from teleseismic events (Beller et al., 2018;Monteiller et al., 2013;Wang et al., 2022), or at the global scale for deep Earth investigations (Masson and Romanowicz, 2017). Target-oriented FWI Willemsen et al., 2016;Yuan et al., 2017) is a computationallyefficient method for such kinds of problems. Fig. 1 illustrates the concept: instead of updating the model parameters in the whole domain (m), we seek to update them only in a small subdomain, indicated by 2 in Fig. 1, when a good approximation of the surrounding medium, indicated by 1, is available.
The conventional methods for target-oriented FWI are based on data redatuming (van der Neut et al., 2015;Yang et al., 2012) or local wave-equation solvers (Borisov et al., 2015;Robertson and Chapman, 2000;van der Neut et al., 2015;van Manen et al., 2007;Willemsen et al., 2016). In the first method, the sources and receivers are moved to 2, hence, reducing the size of the problem and decreasing the computation time. The redatuming needs an integral over the physical source and receivers (Wiggins, 1984); however, the finite sampling of the integrals, limited acquisition coverage, and inaccurate velocity models introduce artifacts in the virtual redatumed data (Haffinger, 2013;Mulder, 2005). The second category of methods do not modify the data acquisition setup, but they try to reconstruct an accurate local wavefield in 2 (Broggini et al., 2017;Gillman et al., 2015;Masson and Romanowicz, 2017;Robertson and Chapman, 2000;Teng, 2003;van Manen et al., 2007). Robertson and Chapman (2000) propose the finite-difference injection method where initially the wavefields on the boundary of 2 are stored after a wavefield simulation in the entire medium. After alteration of 2, the wavefields are updated in 2 by using the stored wavefield, whereas the surrounding model is kept fixed. This method has a straightforward implementation, and it doesn't require any boundary conditions for 2. But, on the other hand, this method doesn't account for the scattering of the updated wavefield occuring in 1 and re-entering 2, which is called high-order long-range interaction by Robertson and Chapman (2000).
Inspired by van Manen et al. (2007); Vasmel et al. (2013), Willemsen (2017) proposed another localized solver method in the frequency domain by applying immersive boundary condition (IBC) to 2. Such IBCs take into account all of the wavefield interactions between altered subdomain 2 and unchanged subdomain 1, but, for each frequency, it requires as many wave-equation solutions as grid points discretizing the edges of the domain of interest for simulating the Green's functions in the entire medium. To keep the computational burden reasonable, Kumar et al. (2019) proposed a rank minimization-based framework to compute a low-rank factorized form of the Green's functions, which requires fewer wave-equation solutions. FWI is classically solved based on a reduced approach formulation with local optimization algorithms (Pratt et al., 1998), the only affordable approach to the problem considering the size of the data and model spaces, particularly in 3D. Therefore, FWI requires an initial subsurface model that is iteratively improved until convergence. The inversion is rapidly trapped in local minima when the initial model does not predict recorded traveltimes to less than half a period (Virieux and Operto, 2009). Overcoming this issue has attracted some attention over the last decade (Gholami et al., 2022, and the references therein). Among such methods, one category consists in extending the linear regime of FWI by enlarging the search space in the state-variable domain (Abubakar et al., 2009;Aghamiry et al., 2019;van Leeuwen and Herrmann, 2013) among others. Considering the bi-convexity of the FWI problem, Aghamiry et al. (2019) propose using the alternating direction method of multipliers (ADMM) (Boyd et al., 2010), which is an efficient method for solving bi-convex optimization problems. This method, which is called iteratively refined wavefield reconstruction inversion (IR-WRI), is an attempt to solve the FWI problem in full-space, i.e., when the wavefields, model parameters, and Lagrange multipliers are updated independently (Rieder, 2021), namely all-at-once methods. IR-WRI updates these three classes of variables in an alternating mode. The method relaxes the requirement to satisfy the wave equation at each iteration to improve the data fit. This amounts to reconstructing wavefields, namely data-assimilated (DA) wavefields, that best jointly fit the observations and satisfy the wave equation in the least-squares sense. Then, the subsurface parameters are updated by minimizing the wave-equation errors generated by the relaxation. Finally, the Lagrange multipliers (dual variables) are updated with a classical dual ascent approach, where they accumulate the remaining errors of the constraints.
In this study, we want to develop a localized version of IR-WRI (LWI) to update 2 when a good approximation of 1 is available. We first do a block decomposition on FWI to separate 1 and 2. By doing so, FWI becomes a constrained optimization problem with four unknowns, wavefields, and model parameters in 1 and 2. This problem can be solved with multi-block ADMM (Liu and Han, 2015), which updates all four variables in an alternating mode. This method will have a slower convergence rate, but it gets some computational gain from solving smaller systems compared to IR-WRI. In LWI, we propose an adaptation of multi-block ADMM, and we only update the model and wavefields in 2 at each iteration, while we update the model and wavefields in 1 only once. The proposed method is straightforward, and it is computationally cheap, but on the other hand, it estimates an approximate wavefield in the updated subdomain 2. Our numerical tests show that this effect is negligible in our formulation because we involve the data when we extract the wavefield in 1. We assess the performance of the method using a time-lapse example with the Marmousi model. The results show that the proposed method is robust to errors in the velocity model in 1. In addition, it can update the velocity model in 2 in a reasonable time with a quality close to the case when we apply conventional IR-WRI on the entire domain.

Theory
Frequency-domain FWI can be formulated as ( 2.1) where m ∈ C n×1 is the model parameter vector (squared slowness); n is the number of discretized points of the model; A(m) ∈ C n×n is the Helmholtz operator; U ∈ C n×ns , B ∈ C n×ns and D ∈ C nr×ns denote the wavefield, source term and the recorded data for n s sources and n r receivers, respectively; P ∈ R nr×n is the observation operator; R(m) is an appropriate regularization function on the model domain; and M is a convex set defined according to our prior knowledge of m (Aghamiry et al., 2020). Finally, A(m) = ∆ + ω 2 Diag(m), ω is the angular frequency, ∆ is the Laplacian, Diag(x) denotes a diagonal matrix with the entries of the vector x on its diagonal.
Suppose subvectors m 1 ∈ C n1×1 and m 2 ∈ C n2×1 are the partitions of vector m related to 1 and 2, respectively, the wave-equation can be decomposed as (Masson and Romanowicz, 2017) where A 1 (m 1 ) ∈ C n×n1 and A 2 (m 2 ) ∈ C n×n2 are over-determined matrices containing the columns of A(m) related to 1 and 2, respectively, and n 1 + n 2 = n. The same decomposition can be applied on P to extract the observation operator in 1 and 2 (P 1 , and P 2 , respectively). Accordingly, FWI can be written as minimize U1,U2,m1,m2∈M, In this study, we assume that the receivers are only in 1, .i.e. P 2 = ∅, which is the conventional setup for localized inversions in exploration seismology. But, generally, there isn't any restriction for the positions of receivers. For example, like teleseismic applications, they can be in 2, .i.e. P 1 = ∅.
LWI is based on multi-block ADMM (Boyd et al., 2010;Liu and Han, 2015) for solving (2.3). (Multi-block) ADMM finds the saddle point of a combination of classical Lagrangian and a penalty term, namely the augmented Lagrangian (AL), in an alternating mode for primal variables and dual variables as is the scaled form of AL, F is the Frobenius norm, • k is the value of • at iteration k, λ > 0 is the penalty parameter assigned to the wave equation constraint, andB k ∈ C n×ns ,D k ∈ C nr×ns are the scaled dual variables. Solving Eqs. (2.4) will converge to the solution of (2.1), but with a slower convergence rate compared to the conventional ADMM (IR-WRI). In LWI, we only want to update the model parameters in 2. To this end, we modify Eqs. (2.4) as follows: 1) Estimating the wavefield in 1: The wavefield in 1 is updated only once, at the first iteration of the inversion of each frequency, because we need it to update the wavefield and model parameters in 2 (Eqs. (2.4b) and (2.4d)). On the other hand, to extract the wavefield in 1 (Eqs. (2.4a)), we need the wavefield and model parameters in 2. At this stage, we can jointly estimate the wavefield in 1 and 2, which is nothing more than the DA wavefield of IR-WRI in the entire of the medium as (2.5) Then U 0 1 and U 0 2 are easily extracted from U 0 . The Lagrange multipliers (B, andD) are initialized with 0, thus these terms are absent in Eq. (2.5).
Having computed U 0 , we can update the model parameters at no extra cost (Aghamiry et al., 2019). This step is not mandatory for the algorithm.
2) Estimating model parameters in 2: Because of the bi-linearity of the wave-equation (Aghamiry et al., 2019), the medium parameters in 2 are updated as where U k+1 2,i refers to the wavefield in 2 related to i th source, L(U k+1 2,i and B i andB k i refer to the i th source and Lagrange multiplier, respectively and ∆ 2 is the Laplacian in 2. Eq. (2.6) can be solved using splitting methods in the presence of non-smooth regularization terms (Aghamiry et al., 2020).

3) Estimating wavefields in 2:
We update the wavefields in 2 by finding the minimizer of Eq. (2.4b), which has a closed form solution (2.7) The term B +B k − A 1 (m 0 1 )U 0 1 is called redatumed data by Masson and Romanowicz (2017), because the sources and receivers virtually are moved onto the boundaries of 2. The important point about LWI is that U 0 1 in Eq. (2.7) is reconstructed using the DA wavefield in Eq. (2.5), which considers the data D for extracting the wavefield. This helps LWI to minimize the interaction of the wavefield between 1 and 2, when 2 changes. The inclusion test (Fig. 2) is designed to show this effect.

4) Updating the dual variables:
At the final step, like IR-WRI, we update dual variables by the running sum of the constraint violations (source and data residuals), Eqs. (2.4e)-(2.4f). For a frequency batch with n b iterations, the proposed method requires n s PDE solutions of size n × n for step 2 and n b * n s PDE solutions of size n 2 × n 2 for step 5, while it is n b * n s PDE solutions of size n × n in IR-WRI. The computational burden of other steps are negligible in LWI and IR-WRI.

Numerical results
Inclusion test: When the model parameters change in 2, the interaction of the wavefield between 2 and 1 are missed if we don't update the wavefield in 1. This simple test is devoted to showing the effects of the DA wavefield U 0 1 of 1 in Eq. (2.7) instead of using the wavefield extracted by solving the wave equation in the full domain. The test involves a 4.2 km × 4.2 km homogeneous velocity model with v = 2 km/s and a circular inclusion with v = 2.8 km/s in the center depicted in Fig. 2a, where the boundary between 1 and 2 is shown with dashed lines. The acquisition consists of 120 receivers along a ring with a radius of 1.9 km and a source shown by a red asterisk in Fig 2a. The 5 Hz wavefield related to this model is shown in Fig. 2b. Suppose we want to extract the wavefield in 2 when the wavefield in 1 (and the values on the boundaries of 2) is calculated in the homogeneous medium with v = 2 km/s, without the  Because the DA wavefield involves the data term when solving the wave equation, it is a better approximation of the true wavefield in 1 (compare wavefields in 1 in Figs. 2b and 2d). Finally, the extracted wavefield in 2 using (2.7) with Figs. 2c-2d are shown in Figs. 2e and 2g, respectively. The difference between these wavefields and the true wavefield in 2 are shown in Figs. 2f and 2h, respectively. The extracted wavefield in 2 using DA wavefield, Figs. 2g-2h, can take into account the interactions between 1 and 2 and hence gives a better approximation of the wavefield in 2.
4D FWI example: Here the goal is to rapidly estimate the local changes that happen because of injected fluids or gas in the subsurface between a baseline and monitor data set (Asnaashari et al., 2012). In the past decade, different methods have been proposed for this problem. Among them, the most straightforward one is to apply conventional FWI on the monitor data when the recovered model from the baseline data is used as the initial model and then subtract the estimated models, which is expensive. To bypass this issue, one can use the velocity extracted from the baseline data as the initial model for applying localized FWI on the monitor data, which can decrease the computational burden significantly compared to the former case . For this test, we use the Marmousi II benchmark (Fig. 3a) as the baseline model. The fixed-spread acquisition contains 57 point sources spaced 300 m apart and the source signature is a 10 Hz Ricker. A line of receivers spaced 50 m apart at the surface is used to record data. We apply the inversion in the 3 Hz -13 Hz frequency band with a frequency interval of 1 Hz. We perform successive mono-frequency inversions proceeding from the low frequencies to the higher ones. We perform three passes through the frequencies, using the final model of one path as the initial model of the next one. The starting and final frequencies of the three paths are [3,5], [3,9], [5,13] Hz, respectively. The initial model and the final estimated model using IR-WRI without regularization are shown in Figs. 3b-3c, respectively. The difference between the true model and the final model of IR-WRI is shown in Fig. 3d. We create the monitor model by adding a small perturbation (Fig. 4a) to the true model. Also, the difference between monitor and baseline model is shown in Fig. 4b. Only the part of the Marmousi model, which is indicated by white rectangular in Fig. 3a is shown in Fig. 4. In this test 2 contains three different parts as indicated by black rectangles in Figs. 3a and 4. We use three frequencies for the monitor data inversion, [5, 10, 15]Hz, and a successive mono-frequency inversion is applied with a stopping criteria of five iterations for each frequency. In addition, we perform two paths through the frequencies, which means that the total number of iterations is 30 for all the inversion results in Fig. 4. The extracted LWI model when the true baseline model (Fig. 3a) is used for the initial model of LWI is shown in Fig. 4c. We continue by applying LWI when the extracted model from IR-WRI of baseline data (Fig. 3c) is used as the initial model. These LWI results demonstrate the impact of the initial models used for the inversion with obvious differences in the largest, deepest rectangle. We repeat the LWI test in Fig. 4d, but we allow for updating the velocity model in 1 one time without extra cost as discussed above (Fig. 4e). Comparing Figs 4d-4e show significant improvements in the reconstructed model. Also, we apply IR-WRI (inversion in the entire of the medium) with the same setup and configuration of the LWI test on the monitor data (Fig. 4f). We observe that LWI and IR-WRI in 2 (black rectangles in Figs. 4e-4f) agree well.

Conclusion
We implemented a localized version of the recently proposed IR-WRI, called LWI, to reduce computational cost for target-oriented applications. The proposed method does a block decomposition on the original FWI problem and decomposes the computational domain to a zone of interest 2 and the rest (1). The new optimization problem is solved using an adaptation of multi-block ADMM when the subproblems related to 1 are solved only once, but the others are solved normally. The proposed method is fast and straightforward, but it only approximately solves the subproblems related to 2 when its velocities change. We show that this error is small because of the ability of the DA wavefield to approximate the true wavefield in 1. Also, we show that the algorithm can update the background medium without additional computational cost when the wavefield of 1 is updated. This update gives a better reconstruction in the targeted area when the velocity in 1 is inaccurate. Numerical tests show promising results with this new algorithm.