Fast Computation of Green Function for Layered Seismic Field via Discrete Complex Image Method and Double Exponential Rules

: A novel computational method to evaluate the Sommerfeld integral (SI) efﬁciently and accurately is presented. The method rewrites the SI into two parts, applying discrete complex image method (DCIM) to evaluate the inﬁnite integral while using double exponential quadrature rules (DE rules) for the computation of the ﬁnite part. Estimation of signal parameters via rotational invariance techniques (ESPRIT) is used to improve the accuracy and efﬁciency of extracting DCIM compared to the generalized pencil of function (GPOF). Due to the symmetry of the horizontal layered media, the Green function, representing the seismic ﬁelds due to a point source, can be written in the form of Sommerfeld integral in cylindrical coordinate system and be calculated by the proposed method. The performance of the method is then compared to the DE rules with weighted average partition extrapolation (WA), which shows a good agreement, with computational time reduced by about 40%.


Introduction
Green function for the horizontal layered seismic field is usually derived by reflectivity method, which was proposed by Fuchs and Müller [1] and extended to many other kinds [2], like reflection and transmission coefficient matrix method [3], discrete wavenumber method [4], discrete wavenumber finite element method [5], and generalized reflection transmission coefficient matrix method [6]. In the frequency domain, the Green function, derived by the reflectivity method, can be written in the Sommerfeld integral form in cylindrical coordinate system for symmetrical media.
It is well known that the numerical evaluation of Sommerfeld integral (SI) is computationally expensive due to the oscillatory and slow convergence of the integrands. To overcome this problem, several approaches have been proposed, which can be divided into two main categories: one is the approximation of the spatial domain Green functions in a closed form where no numerical integration is needed, and the other is the numerical integration of SI in conjunction with some acceleration techniques [7]. Within the first category, discrete complex image method (DCIM), which approximates the integrand of Sommerfeld integral by a series of complex exponential functions, is commonly used for the advantages of high computational efficiency, but it needs to handle the surface wave poles contributions, which not only makes the computation complicated but also brings singularity to the near region [8], and also the calculation accuracy and effective range are difficult to be accurately estimated. For the latter category, the common practice is dividing the whole Sommerfeld integral into two parts: the first part is the path to bypass the singularity; the second part, the path to infinity, is the Sommerfeld tail integral. The finite-range integrals may readily be evaluated by the Gauss-Jacobi quadrature [9] or by the double-exponential (DE) rule [10,11]. Mosig first used DE rules to calculate Sommerfeld integrals in [11,12] which indicated its validity of suppressing endpoint singularities. The calculation of tail integral is difficult to converge due to Bessel's oscillation and slow attenuation characteristics, so this kind of method generally requires extrapolation to accelerate convergence. The WA method has shown higher levels of convergence among various extrapolation methods [13][14][15]. This kind of method does not need to strictly locate the position of singularity in Sommerfeld integral but only needs to ensure that the first integral path avoids all the singularities. It has good adaptability and controllable numerical accuracy; however, this depends on the number of intervals (n) that are chosen to evaluate the tail region. The computational time also rapidly increases as the value of (n) increases [16]. Given the advantages and disadvantages of these two methods, we propose a method by combining DE rules and DCIM to calculate Sommerfeld integral.
This paper first presents the Green function of a point source in a multilayer half-space in Section 3, explaining the mathematical manipulation required to obtain the solution as a Sommerfeld integral form in the frequency domain; it then describes the principle of the proposed method with DE quadrature rules and DCIM in Section 4; finally, it corroborates the correctness of the algorithm by the frequency responses obtained from the proposed approach with those where DE rules and WA partition-extrapolation are used for half-space model, and the finite element method is used for three layers model.

Seismic Wave Equation
The propagation of seismic wavefield in the time domain can be simplified by the following three-dimensional acoustic wave equation: where u (x, y, z, t), v (x, y, z, t), and f (x, y, z, t) represent displacement, velocity, and source term, respectively. f (x, y, z, t) = −δ(x − x s , y − y s , z − z s )s(t), s(t) is the wavelet, and Ricker wavelet is used in this paper; and δ(x − x s , y − y s , z − z s ) is the Dirac function at the source point (x s , y s , z s ). By Fourier-transform of Equation (1), the two-dimensional acoustic wave equation in frequency domain is obtained, and therefore, Green function for the problem is defined by the following equation: where, G denotes the Green function, F(x, y, z, ω) = −δ(x − x s , y − y s , z − z s )S(ω) is the source term in the frequency domain, k(x, y, z) = ω/v is wave number, S(ω) is Ricker wavelet in the frequency domain, and ω is the angular frequency. However, the underground medium is always viscous, which leads to wave energy loss and phase change in the process of propagation. The visco-acoustic wave equation is established to better describe the propagation of the seismic waves in this viscous medium, which is the same form as Equation (2), but the complex velocity is introduced to simulate the viscous effect [17][18][19]. The reciprocal of complex velocity is defined as [18], where Q is quality factor, so the complex wavenumber is set to be In this paper, the value Q generated by Li Qingzhong's empirical formula is used for numerical simulation [20]

Green Function in Full-Space
The Green function of Equation (2) in homogeneous full-space can be expressed as is Ricker wavelet in the frequency domain, ω is the angular frequency, and k 1 is the wavenumber of the medium. The above formula can be rewritten in the Sommerfeld integral form in cylindrical coordinate system:

Green Function in Layered Half-Space
Consider n layers symmetric structure of homogeneous medium defined by interfaces located at z 1 , z 2 , · · · , z n−1 , as shown in Figure 1. The density of each layer is ρ 1 , ρ 2 , · · · , ρ n ; and the velocity is v 1 , v 2 , · · · , v n . In this paper, the source is placed in the second layer.
In this paper, the value Q generated by Li Qingzhong's empirical formula is used for numerical simulation [20]

Green Function in Full-Space
The Green function of Equation (2) in homogeneous full-space can be expressed as S ω is Ricker wavelet in the frequency domain, ω is the angular frequency, and 1 k is the wavenumber of the medium. The above formula can be rewritten in the Sommerfeld integral form in cylindrical coordinate system:

Green Function in Layered Half-Space
Consider n layers symmetric structure of homogeneous medium defined by interfaces located at 1 2 1 , , , n z z z −  , as shown in Figure 1. The density of each layer is 1 2 , , , n ρ ρ ρ  ; and the velocity is 1 2 , , , n v v v  . In this paper, the source is placed in the second layer. Each layer satisfies the acoustic Equation (2) with parameters of the Green function G , wavenumber k , velocity v , density ρ , layer thickness h , and quality factor Q , respectively. We obtain the equations as follows: Each layer satisfies the acoustic Equation (2) with parameters of the Green function G, wavenumber k, velocity v, density ρ, layer thickness h, and quality factor Q, respectively. We obtain the equations as follows: where At the interface, the pressure P i = ρ i G i as well as the gradient of the potential for the vertical direction ∂G i ∂z i are continuous [21]. Therefore, the following boundary conditions can be imposed on the Green function The solutions of Equations (6) and (7) can be regarded as the summation of separate up-going and down-going waves, and therefore, it can be written in the form of Sommerfeld integral in cylindrical coordinate system as follows: where Equations (9) and (10), e m i z and e −m i z may be infinity. To maintain numerical stability, rewrite Equations (9) and (10) into By using the boundary conditions (8), the unknowns C 1 , C 2 , D 2 , · · · , C i , D i , · · · , C n−1 , D n−1 , D n in the above formula are solved. The coefficients of the source layer are derived firstly, and other coefficients can be obtained by recursion; then, the expression of the Green function of the layered medium is obtained. where

Methods
The expressions (11) and (12) contain Sommerfeld integrals, which is an infinite integral with the highly oscillatory and slow-decaying kernel. In this paper, the partial closed form of Sommerfeld integral is derived, and ESPRT is applied to extract DCIM, while DE rules are used for the computation of the finite integration.

Partial Closed Form Expression
The Sommerfeld integration in Equations (11) and (12) can be written in the following form: Symmetry 2021, 13, 1969

of 12
The integration is then divided into two parts, (18) where p is a reasonably selected integral breakpoint. The second integral on the right side of the Equation (18) can be asymptotically approximated and then be written as On substituting (19) in (18), we get The kernel function g(m i ) in (17) can be approximated by an exponential function, where pb is the number of exponentials used for approximation. In this paper, the coefficients a (l), b (l) are solved by the ESPRIT algorithm, which will be discussed in the next section. From the Sommerfeld identity expressed by formula (5), the closed form of the second integral can be obtained, where The first part of the g(m i ) function usually contains singularity, and the tail is smooth and decays fast. Therefore, if the appropriate p value is selected, the approximate fitting will be very accurate by avoiding the singular value in the front part. The finite integral with singularity in (20) can be evaluated directly by the numerical integration method. We choose DE quadrature rules here for integration with the advantage of dealing with singular points and high precision.

ESPRIT Algorithm
The signal g(m i ) can be sampled as [8] The value of T 01 can be selected from 100~200, T 02 usually set between 1~3, and in this paper, T 01 = 200, T 02 = 2, t is an integer. According to the relationship m 2 = m 2 i + k 2 i , the first sampling in m-plane is k i 1 + T 2 02 . Therefore, the approximation of g(m i ) starts from m = k i 1 + T 2 02 , and the parameter p in formula (20) should be set not less than k i 1 + T 2 02 to insure the integration accuracy. Then the sampling sequence can be expressed as The relation of A(l), B(l), and a (l), b (l) can be obtained from Equations (23) and (24).
Then, the unknown coefficients a (l), b (l) in (24) are obtained.
For the sampling sequence y(0), y(1), · · · , y(N − 1), a data matrix can be constructed: where N is the sample number, L is called the pencil parameter, and its value should be between N/3 and N/2 [22].The data matrix Y can be decomposed by SVD, Find the eigenvalues λ l of the matrix Ψ, For N sampled signals, and According to the least square method, we can obtain So far, the coefficients A (l) and B (l) are obtained. According to (22), (26), and (27), the second part of Equation (20) is solved, and the integrand of the first part is gained, and then DE rules are applied to compute the finite integral.

DE Rules
The double exponential transformation was first proposed by Takahasi and Mori in 1974 [16]. It can be seen from Equations (11)-(16) that the integration kernel is singular at m i = k i . DE rules are insensitive to endpoint singularity and simple to program since the weights and nodes are easily generated [12].
Consider the following form of integral: Let a variable transform: The DE rules are transformed by the tanh-sinh formula, The standard trapezoidal rule for numerical integration is applied with h as grid interval when the integral is defined on the interval [−∞, +∞], and n is the sample point, which is truncated at ±N. Then we can approximate the definite integral via with the nodes ξ k and weights ω k defined as where For arbitrary integral interval [a, b] may be mapped onto [−1, 1] by the linear transfor- Hence, for an arbitrary interval [a, b], the nodes and weights become σξ k + γ and σω k , respectively, and (43) is transformed to As with any other quadrature rule, singularities of the integrand near the integration path adversely affect the convergence. However, any singularities on the integration path Symmetry 2021, 13, 1969 8 of 12 are easily treated by splitting the integration range so that the singularities are placed at the endpoints [12]. Therefore, the integration path in the first part of Equation (20) should be separated into (45) to ensure the convergence of the integration.
where breakpoint bk is set to the real part of wavenumber k i in this paper.

Results
In this section, the half-space model is designed to test the correctness of the proposed method by comparing it with DE rules and partition-extrapolation WA algorithm [12], which was well accepted and known as one of the most accurate and efficient ones. Finally, the DE_DCIM method is utilized for the calculation of the Green function of the three-layer model in comparison with the finite element method [24].

Half-Space
It is assumed that the size of the study area is 1410 m * 1410 m * 710 m; the sampling interval in the horizontal and vertical directions is 10 m, taking the main frequency of 20 Hz Ricker wavelet as the source; and the simulation frequency is 10 Hz (all the study areas in the following are consistent with this study area). Consider a half-space defined by interface z 1 = 350 m, with a point source located at (x s , y s , z s ) = (710 m, 710 m, 360 m). Assume that the velocity and density parameters are set as v 1 = 340 m/s, ρ 1 = 0.00129 g/cm 3 , v 2 = 2000 m/s, ρ 2 = 1.5 g/cm 3 . The integration path is divided into three segments, the breakpoint bk= real(ω/ v i ), and the breakpoint p is set as k i 1 + (198 * 5 * T 02 /200) 2 . Figure 2 shows the comparison of ESPRIT and GPOF methods to extract DCIM of the lower half-space. g(m 2 ) is defined by (13) and (14); sampling point t is defined in (23). From Figure 2a, both methods gain a good fit with the original data, but the ESPRIT method is slightly more accurate, with a fitting error less than 0.008, while GPOF [25] is less than 0.014. Further, the time cost in computing half-space DCIM is also presented in Table 1, which shows ESPRIT also reduces the calculation time. The more layers there are, the more obvious time saving will be seen.
[ ] As with any other quadrature rule, singularities of the integrand near the integration path adversely affect the convergence. However, any singularities on the integration path are easily treated by splitting the integration range so that the singularities are placed at the endpoints [12]. Therefore, the integration path in the first part of Equation (20) should be separated into (45) to ensure the convergence of the integration.
where breakpoint bk is set to the real part of wavenumber i k in this paper.

Results
In this section, the half-space model is designed to test the correctness of the proposed method by comparing it with DE rules and partition-extrapolation WA algorithm [12], which was well accepted and known as one of the most accurate and efficient ones. Finally, the DE_DCIM method is utilized for the calculation of the Green function of the three-layer model in comparison with the finite element method [24].

Half-Space
It is assumed that the size of the study area is 1410 m * 1410 m * 710 m; the sampling interval in the horizontal and vertical directions is 10 m, taking the main frequency of 20 Hz Ricker wavelet as the source; and the simulation frequency is 10 Hz (all the study areas in the following are consistent with this study area). Consider a half-space defined by  Figure 2 shows the comparison of ESPRIT and GPOF methods to extract DCIM of the lower half-space. 2 ( ) g m is defined by (13) and (14); sampling point t is defined in (23).
From Figure 2a, both methods gain a good fit with the original data, but the ESPRIT method is slightly more accurate, with a fitting error less than 0.008, while GPOF [25] is less than 0.014. Further, the time cost in computing half-space DCIM is also presented in Table 1, which shows ESPRIT also reduces the calculation time. The more layers there are, the more obvious time saving will be seen.   Figure 3 shows the symmetrical wavefield of (x, y s , z) plane. Comparing the solution of DE_DCIM and the numerical integration with DE_WA, the relative errors are shown in (c) and (f). Excellent agreement is obtained with a relative error of real part less than 2.5 × 10 −3 and of image part less than 5.8 × 10 −4 , which assesses the validity of the present method. The calculation time of the half-space with different parameters is shown in Table 2. It can be seen from Table 2 that the proposed DE_DCIM method reduced the computational time by about 40% when compared to the DE_WA method, with ρ 1 = ρ 2 = 1.0 g/cm 3 .

Method
GPOF ESPRIT Computation time 0.135 s 0.025 s Figure 3 shows the symmetrical wavefield of ( , , ) s x y z plane. Comparing the solution of DE_DCIM and the numerical integration with DE_WA, the relative errors are shown in (c) and (f). Excellent agreement is obtained with a relative error of real part less than 3 2.5 10 − × and of image part less than 4 5.8 10 − × , which assesses the validity of the present method. The calculation time of the half-space with different parameters is shown in Table 2. It can be seen from Table 2 that the proposed DE_DCIM method reduced the computational time by about 40% when compared to the DE_WA method, with 1

Three-Layer Model
Consider a three-layer structure defined by two interfaces placed at

Three-Layer Model
Consider a three-layer structure defined by two interfaces placed at z 1 = 200 m, z 2 = 500 m, with a point source located at (x s , y s , z s ) = (710 m, 710 m, 200 m). Assume that the velocity and density parameters are set as v 1 = 1000 m/s, ρ 1 = 1.5 g/cm 3 , v 2 = 2000 m/s, ρ 2 = 2.0 g/cm 3 , v 3 = 3000 m/s, ρ 3 = 3.0 g/cm 3 . The method in this paper is used to solve the layered model, and the symmetrical wavefield, as shown in Figure 4, is obtained. It can be seen from (a) and (d) that in the first layer, there are only up-going wave fields; in the third layer, only down-going wave fields; and in the second layer, there are upward and downward wave fields, which are mixed. Comparing with the FEM [24], the results are consistent in shape, and there are some numerical differences. The relative error of the real part is less than 0.09, and the imaginary part is less than 0.04. Figure 4 also indicates that it is correct to calculate the layered space wave field according to Formulas (11) and (12).
in the third layer, only down-going wave fields; and in the second layer, there are upward and downward wave fields, which are mixed. Comparing with the FEM [24], the results are consistent in shape, and there are some numerical differences. The relative error of the real part is less than 0.09, and the imaginary part is less than 0.04. Figure 4 also indicates that it is correct to calculate the layered space wave field according to Formulas (11) and (12). . The three-level DCIM with surface wave extraction [26] is adopted. It can be observed in the figure that the DE_DCIM result is more accurate than the three-level DCIM for about two orders of magnitude when both compared to DE_WA, especially when the distance between source and field point is large. As discussed in reference [27], for multilayer media, it is very difficult to find surface wave poles, and the inaccurate extraction of the surface wave will bring unpredictable errors to the results.  To show the advantage of DE_DCIM over DCIM for accuracy, consider the three-layer structure defined above with a point source located at (x s , y s , z s ) = (0 m, 710 m, 200 m). Figure 5 compares the solutions of DE_WA, DE_DCIM, and DCIM in line (x, y s , z s ). The three-level DCIM with surface wave extraction [26] is adopted. It can be observed in the figure that the DE_DCIM result is more accurate than the three-level DCIM for about two orders of magnitude when both compared to DE_WA, especially when the distance between source and field point is large. As discussed in reference [27], for multilayer media, it is very difficult to find surface wave poles, and the inaccurate extraction of the surface wave will bring unpredictable errors to the results.
in the third layer, only down-going wave fields; and in the second layer, there are upward and downward wave fields, which are mixed. Comparing with the FEM [24], the results are consistent in shape, and there are some numerical differences. The relative error of the real part is less than 0.09, and the imaginary part is less than 0.04. Figure 4 also indicates that it is correct to calculate the layered space wave field according to Formulas (11) and (12).  . The three-level DCIM with surface wave extraction [26] is adopted. It can be observed in the figure that the DE_DCIM result is more accurate than the three-level DCIM for about two orders of magnitude when both compared to DE_WA, especially when the distance between source and field point is large. As discussed in reference [27], for multilayer media, it is very difficult to find surface wave poles, and the inaccurate extraction of the surface wave will bring unpredictable errors to the results. s s x y z of 3-layer model. Figure 5. The results of log 10 |G| and relative errors in line (x, y s , z s ) of 3-layer model.

Conclusions
Sommerfeld integral is included in the Green function for the seismic field in horizontal layered half-space. The numerical technique is used to compute the Sommerfeld integrals by deriving the integral into two parts, the infinite integral part and the finite integral part, and by applying DE quadrature rules to evaluate the finite part and DCIM to calculate the infinite part. Compared with the DE_WA method, the new method can get an accurate result with a relative error less than 2.5 × 10 −3 and increase time saving by about 40%. The ESPRIT method is introduced to extract DCIM for better accuracy and efficiency. Finally, the fast numerical method is applied to the calculation of the seismic field in horizontal layered half-space. The method in this paper takes both efficiency and accuracy into account, and theoretically, higher accuracy can be achieved by controlling the parameters.