On the applicability of a renormalized Born series for seismic wavefield modelling in strongly scattering media

Scattering theory is the basis for various seismic modeling and inversion methods. Conventionally, the Born series suffers from an assumption of a weak scattering and may face a convergence problem. We present an application of a modified Born series, referred to as the convergent Born series (CBS), to frequency-domain seismic wave modeling. The renormalization interpretation of the CBS from the renormalization group prospective is described. Further, we present comparisons of frequency-domain wavefields using the reference full integral equation method with that using the convergent Born series, proving that both of the convergent Born series can converge absolutely in strongly scattering media. Another attractive feature is that the Fast Fourier Transform is employed for efficient implementations of matrix–vector multiplication, which is practical for large-scale seismic problems. By comparing it with the full integral equation method, we have verified that the CBS can provide reliable and accurate results in strongly scattering media.


Introduction
The integral equation (IE) method based scattering theory (Aki & Richards 1980;Zhdanov 2002) is a powerful tool in the modeling of wave propagation, which has a wide application in data processing (Weglein et al. 1997(Weglein et al. , 2003, modeling (Innanen 2009), and seismic inversion (Wu & Zheng 2014;Snieder 1990;Berkhout 2012;Zhang & Weglein 2009;Alkhalifah & Wu 2016;Alkhalifah 2016;Wu & Alkhalifah 2017;). An attractive features of the IE method is that only the anomalous volume (scattering volume) needs to be discretized, which leads to more efficient computation (Malovichko et al. 2017). The implementation of the IE method involves dividing the medium into background and anomalous parts (Zuberi & Alkhalifah 2014).
The Born series has an assumption of weak scattering (Wu & Toksoz 1987;Kouri & Vijay 2003). Convergence issues may occur in strongly scattering areas, such as salt structures. It is important for seismic imaging in such strong-contrast regions to address the weak-scattering assumption. One important approach that addresses the divergence problem is to renormalize the Born series using various renormalization approaches (Eftekhar et al. 2018).
There are several approaches to develop a convergent scattering series. There have been successful attempts to introduce the De Wolf approximation (De Wolf 1971, 1985 into seismic scattering series (Wu & Huang 1995). The renormalized scattering series is derived by Jakobsen & Wu (2016) using the T-matrix and De Wolf series. The T-matrix is a transition operator that includes all the effects of multiple scattering. Renormalization group method has been applied to seismic waveform inversion (Wu et al. 2015) and envelope inversion . Signicant progress has been made by Yao et al. (2015) by dividing the renormalized Lippmann-Schwinger equation into two sub-Volterra type integral equations and introducing wavefield separation.
Recently, by employing the renormalization group (RG) theory, we developed a renormalized version of the Born series. Numerical tests showed that this solution can be convergent for large-contrast media ( Jakobsen et al. 2018). Our renormalization group approach is based on the use of an auxillary set of scale-dependent scattering potentials, which gradually evolves toward the real physical scattering potential.
Another interesting convergent Born series (CBS) was proposed by Osnabrugge et al. (2016) to solve the Helmholtz equation. The convergent Born series can guarantee convergence by localizing the wavefields, in which the contracted preconditioner must be specified. Actually, the CBS can be understood as a kind of renormalized Born series based on the RG theory. In the early 1970s, the renormalization procedure was proposed by Gell-Mann & Low (1954) for problems of infinity and divergence. From the early 1970s (Wilson 1971), RG theory has been widely used to remove divergence in quantum physics, critical phenomena, dynamical systems and statistical mechanics, etc. The major purpose of the RG theory is to obtain stable properties of physical systems (Goldenfeld 1992). Based on the above fact, Chen et al. (1994Chen et al. ( , 1996 applied the RG approach to deriving global asymptotic solutions of differential equations. Since then, the RG theory has been well developed and significant progress in renormalizing perturbation series (Yakhot & Orszag 1986;Pelissetto & Vicari 2002;Delamotte 2004;Kirkinis 2008Kirkinis , 2012 has been made. The purpose of this paper is to extend the CBS of Osnabrugge et al. (2016) to seismic scattering problems for strongly scattering media and compare the CBS with the full integral equation method. The convergent Born series is obtained by localizing the Green's function with a dampling factor. From the technical point of view, the CBS removes the divergence by localizing the wavefields and controlling convergence using a preconditioner. Thus, the CBS can be understood as a kind of renormalized Born series. After presenting the convergent Born series, we analyze the theoretical background of the convergent Born series from the renormalization group theory prospective and its nature of localization. Then, we give numerical results of frequencydomain wavefields. To quantatively compare the results from the CBS and full integral equation methods, we present numerical results for results of pressure wavefields in strongly scattering media.

The Lippmann-Schwinger equation
The Helmholtz equation can be written as (Morse & Feshback 1953): with the wavenumber k and source signal s(r), where ∇ is the second-order differentiation. Defining b (r) as the background field, where k b is the background wavenumber. The total wavefield (r) can expressed as where k s is the scattered wavenumber and k s = k − k b . Substituting equation (3) into (1), we have the following equation for the scattered wavefields s (r): From equations (2) and (4), we get the scattered wavefields s (r) where the background Green's function G b can also be calculated with analytical expressions for homogeneous media, and the ray theory (Cerveny 2005; Huang & Greenhalgh 2019), Gaussian beam (Huang et al. 2016a(Huang et al. , b, 2018Huang 2018) or finite difference method (Carcione 2007) for inhomogeneous media. Finally, we can get the Lippmann-Schwinger equation:

The conventional Born series
Equation (7) has the formal solution as This equation can be solved iteratively: (r g , r s , k) = b (r g , r s , k b ) After taking the first-order term of the conventional Born series, we have the Born approximation: (r g , r s , k) = b (r g , r s , k) The Born series is the basis for seismic forward and inverse problems. However, because the Born series assumes weak scattering, it can only converge when the scattering potential is weak (Kirkinis 2008;Wu et al. 2007). For real applications, the strength of the contrast in the medium is relatively strong. Divergence may occur in media with strong contrasts.

Convergent Born series
This section revisits the equations for the CBS (Osnabrugge et al. 2016). Equations are presented for an arbitrary strong medium. The convergent Born series refers to the situation in which the value of the coefficient of each iteration term is less than unity. The key point for modification of the conventional Born series is to introduce a damping parameter and a preconditioner . Here, we review the derivation of the CBS. To this end, we start with the conventional Born series and apply the preconditioner to both sides of the conventional Born series. It should be noted that the damping parameter , which is related to the attenuation of the wavefields in the background medium, is also important for the convergence.
In operator form, the Born series (9) can be written as follows: where S is the source term, which represents the source wavelet in frequency-space domain, and G is the Green's operator, which represents the Green's function G b in equation (9). Note that S and are vectors, and G is dense operator filled with various Green's functions. is a diagonal matrix. By applying a preconditioner to both sides of equation (11), Osnabrugge et al. (2016) Then, we reformulate equation (12) as with where The combination of the parameters and can ensure that the largest eigenvalue of Λ is smaller than unity. The details of choice of parameter will be discussed in the section Implementation. The modified Born series is explained by the following renormalized Born series: The iteration form solution is = Λ + GS with the initial solution b = GS. The background Green's function can be written as (Osnabrugge et al. 2016) and in the real-space domain (Osnabrugge et al. 2016) where k is the wavenumber vector in the real medium.
Because of the introduction of the coefficient of each iteration term Λ and preconditioner into the Born series, the modified Born series is convergent. Mathematically, Osnabrugge et al. (2016) demonstrated that the modified Born series can converge by combining the two parameters and given the suggestions of choice of the parameters in the optimum scale.
Here, we aim to apply the CBS method to seismic wave modeling problems and provide the renormalization interpretation of the convergent Born series. Due the strong contrast in the seismic problems, it is more challenging. We perform numerical tests and investigate how to choose the coefficient Λ and preconditioner for each iteration term, and the dependence on the parameter . We will look at how the parameter changes for different models with different strong scattering cases.

Interpretation of the convergent Born series
In this section, we discuss the convergent Born series from the renormalization perspective.

Step-by-step local interaction
The modified version of the Born series developed by Osnabrugge et al. (2016) is called the convergent Born series. In this modified Born series, a preconditioner is introduced. By combining the preconditioner and parameter , the iteration computation satisfies the convergence condition of the Born series.
Actually, the concepts of the locality of wavefields explain how step-by-step propagators prevent the CBS from divergence. From equation (13), we can see that for the conventional Born series, each term involves integrations over the whole volume, which leads to the divergence problems of strongly scattering medium. The CBS makes the total energy in the background medium localized and finite so the volume integral in each term will not blow up. It compensates the damped wavefield by introducing an imaginary part with an opposite sign into the scattering potential V. This means that in the latter procedure the wavefield will grow when interacting with the scattering potential, and therefore compensate the energy loss during propagation in the background medium. In this way, those interactions always act locally, and thus can be regarded as short-range interactions. The iterations will continue until the wavefields cover the whole region with accepted accuracy. Physically, this can be explained as the renormalization process. According to Wilson's RG theory, one can first integrate out the local interactions and then derive the effective action operator, and then go to the next level to calculate the local interactions based on the effective interaction operator. The RG procedure in CBS is more like the mathematician's renormalization procedure  (Chen et al. 1996) using RG theory as a floating initial condition. At each step, the calculated field is treated as a new initial wavefield for further propagation. This is different from Wilson's multi-scale RG procedure.

Implementation
In this section, we give the coordinate representation for the CBS. The advantage of this representation is that it is easy to relate two adjacent scattering potentials. Using the coordinate representation, equation (13) can be rewritten as (r) = (r) G (r) S + Λ (r, r ′ ) (r).
The entire model is discretized into N x × N z in the twodimensional case. Then, we have where i = 1,…, N x and j = 1,…, N z .
It should be noted that, to compute the wavefields at the receivers along the surface, we need to compute the wavefields from the sources to any subsurface point and the background Green's functions from receivers to any subsurface point. For computing the wavefields at any subsurface point, we need compute the wavefields in the background medium and the Green's functions from any subsurface point to any subsurface point (the Green's function from volume to volume G VV ). Then we need compute Λ, which is used for the high-order terms of the CBS. The workflow for implementing the renormalized Born series can be found in Algorithm 1.

Choice of parameter in the scale of seismic modeling
An important issue for the CBS is to choose the parameter . From the analysis in the above section, it can be found that the stronger the scatters (large-contrast), the higher the required parameter . This is because to eliminate the divergence the wavefields should be strongly localized. From equation (19), one can see that the higher the parameter , the stronger the attenuation of the background Green's functions. This means that there is a compromise between the demand on the convergence of the CBS and the computational cost. After conducting numerical tests, we find that the parameter should be chosen as follows: After investigating the convergence of the convergent Born series, Fast-Fourier Transform (FFT) is used, which can accelerate the computation in the implementation. Following Osnabrugge et al. (2016), we employ the FFT technique for an efficient matrix-free implementation. The FFT method has been used for integral equation modeling (Liu et al. 2001;Gao & Torres-Verdin 2006). The matrix-vector multiplications can be expressed as where  is the forward 2D FFT operator and  −1 is the inverse 2D FFT operator. It should be noted that the product is performed in the size of 2N x × 2N z . The computational complexity is O[N x N z log(N x N z )] and the memory complexity is O[N x N z ].

Comparison of frequency-domain wavefields
In this section, we share the frequency-domain wavefields for different models, including two-layers and the SEG/EAGE salt models as well as, compare the convergence property of the CBS by calculating the normalized errors and share the pressure response along receiver line for different iterations. To demonstrate the accuracy of the CBS, we compare it with the full integral equation method ( Jakobsen & Wu 2016). We have used homogeneous background media in the tests and the velocity is 2000 m s −1 . We first compare the renormalized Born series against the T-matrix method in an acoustic two-layer model (figure 1). The model measures 1280 m × 1280 m with grid intervals for the simulations of 10 m × 10 m.
Snapshots of frequency-domain wavefields computed with the full integral equation method and the convergent Born series are shown in figures 2 and 3, respectively. In each figure, we show the wavefields of two frequencies, 10 Hz and 15 Hz. From figures 2 and 3, one can make the following observations: (1) the wavefields using all the methods display an obvious change around the boundary; (2) the wavefields using the CBS match well with those from the reference integral equation method. Figure 4 shows the resampled SEG/EAGE salt model for this example. The model grid is 10 m × 10 m. The model represents a uniform mesh of 128 × 128 nodes. We have performed simulations of frequency-domain wavefields in which the frequencies of 10 Hz and 15 Hz are used. Figure 5 shows the wavefields for 10 Hz and 15 Hz obtained by the reference integral equation method with 100 iterations. Figure 6 shows the wavefield snapshots for 10 Hz obtained by the conventional Born series (BS) method with 20, 50, 80 and 100 iterations. Figures 7 and 8 show the wavefields at frequencies 10 Hz and 15 Hz, respectively, obtained by the CBS method with 20, 50 ,80 and 100 iterations.

Convergence property of the CBS
To investigate the convergence property of the CBS, we calculate the normalized error. Figures 9 and 10 show the results for the two-layered and the SEG/EAGE salt models, respectively . The figures show that the CBS has a similar convergence property in different models and frequencies, but the error decreases in a different way. With the same iterations, the error in the two-layered model is smaller than those of the SEG/EAGE salt model. From the figures, one can observe that after around 100 iterations the error of the CBS is very small. This suggests that the CBS can give a good match with the reference solution.

Frequency-domain wavefields with FFT
Because we use FFT in the implementation, some periodic boundary condition problems may occur (Osnabrugge et al. 2016). To prevent the reflection from the boundaries, we use an absorbing boundary condition in the implementation 12 RG theory: real(G) of the CBS. The absorbing boundary condition has been used in the context of wave modeling by different authors. More specifically, for the CBS, we use a width of absorbing boundary layer with grids of 40 × 40. Figure 11 shows the frequency-domain wavefields of 10 Hz using the CBS with FFT. We have estimated the computational cost. The computational times of the reference integral equation method is 492 s. The computational time of the CBS with FFT is 100 s.

Anomalous pressure response along receiver line
We now consider the same simulation as the SEG/EAGE salt models but for the synthetic pressure response along the receiver line. This example is designed to test the accuracy of the numerical scheme. Figure 12 shows the results for a twolayered model, in which the frequencies of 10 Hz and 15 Hz are used. Figures 13 and 14 show the pressure response in a salt model using the BS method with the frequency of 10 Hz for different iterations. Figures 15 and 16 show the pressure using the CBS method with the frequency of 10 Hz for the salt model. Figures 17 and 18 show the pressure using the CBS method with the frequency of 15 Hz for the salt model. For all the tests, the point source is placed at the same position. A receiver line is located at the surface. From figures 13 and 14, one can observe that the results from the BS method do not agree with the pressure wavefields from the full integral equation method. From figures 15-18, we observe that the pressure response using the CBS works very well compared with the result using the reference integral equation method.

Discussion
Before we discuss the convergence, computational complexity and potential application of the CBS, we would like to clarify that we have presented the theory and performed the numerical tests in the frequency domain because scattering theory is naturally formulated in the frequency domain and we do not have to generate time-domain waveforms to perform inversion in the frequency domain. The main reason for using a homogeneous reference medium is that we want to use FFT, which depends on the fact that the Greens function for a homogeneous medium is directly related to the difference between x and x ′ . Another point is that the contrast is frequency-dependent. In our tests, we used different frequencies and investigated different choices of parameter .
The application of the BS to seismic forward modeling requires small contrasts to achieve convergence. Here, by applying a preconditioner to the both sides of the BS and introducing the parameter to the background Green's function, the convergence of the BS is guaranteed. Figures 19  and 20 show the difference of frequency-domain wavefields for SEG/EAGE salt model using full integral equation and CBS methods. Osnabrugge et al. (2016) have already provided a general proof of convergence, we have used numerical examples to verify that the general proof holds for our specific models. However, for the case where the contrast is very large, e.g. salt areas, more iterations are needed to achieve convergence. This is also related to the choice of the parameter . The stronger the scatters (large-contrast), the higher the needed parameter .
Compared to the conventional BS there is no additional computational cost for each term in the real-space implementation. It should be noted that the accuracy of the wavefields depend on the number of iteration. This is different from the full IE method. One important thing we would like to mention is that, due to the FFT implementation for the CBS, the computational cost is reduced significantly. The computational complexity for such an implementation stone to developing modifications for one-way propagators. Also, it can be used to establish the Frechet derivatives for multi-scattering full waveform inversion (Alkhalifah & Wu 2016).

Conclusions
Seismic scattering theory is an effective method for seismic wave modeling and is the basis of seismic inversion. However, the Born series assumes weak scattering, which