Next Article in Journal
Improving the Convergence of Interval Single-Step Method for Simultaneous Approximation of Polynomial Zeros
Previous Article in Journal
Matching of Manufacturing Resources in Cloud Manufacturing Environment
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
School of Geosciences and Info-Physics, Central South University, Changsha 410083, China
2
Key Laboratory of Metallogenic Prediction of Nonferrous Metals and Geological Environment Monitoring, Ministry of Education, Changsha 410083, China
3
College of Earth Sciences, Guilin University of Technology, Guilin 541004, China
*
Author to whom correspondence should be addressed.
Symmetry 2021, 13(10), 1969; https://doi.org/10.3390/sym13101969
Submission received: 23 September 2021 / Revised: 13 October 2021 / Accepted: 17 October 2021 / Published: 19 October 2021
(This article belongs to the Topic Applied Metaheuristic Computing)

Abstract

:
A novel computational method to evaluate the Sommerfeld integral (SI) efficiently and accurately is presented. The method rewrites the SI into two parts, applying discrete complex image method (DCIM) to evaluate the infinite integral while using double exponential quadrature rules (DE rules) for the computation of the finite part. Estimation of signal parameters via rotational invariance techniques (ESPRIT) is used to improve the accuracy and efficiency 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 fields 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%.

Graphical Abstract

1. 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.

2. Seismic Wave Equation and Green Function

2.1. Seismic Wave Equation

The propagation of seismic wavefield in the time domain can be simplified by the following three-dimensional acoustic wave equation:
2 u ( x , y , z , t ) = 1 v ( x , y , z ) 2 2 u ( x , y , z , t ) t 2 + f ( x , y , z , t )
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:
2 G ( x , y , z , ω ) + k 2 G ( x , y , z , ω ) = F ( x , y , z , ω )
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 1 v ˜ = 1 v 1 j 2 Q [18], where Q is quality factor, so the complex wavenumber is set to be k = ω v 1 j 2 Q , j = 1 . In this paper, the value Q generated by Li Qingzhong’s empirical formula is used for numerical simulation [20]
Q = 14 × ( v / 1000.0 ) 2.2

2.2. Green Function in Full-Space

The Green function of Equation (2) in homogeneous full-space can be expressed as
G ( x , y , z , ω ) = S ( ω ) e i k 1 R 4 π R
where R = ( x x s ) 2 + ( y y s ) 2 + ( z z s ) 2 , S ( ω ) 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:
S ( ω ) e i k 1 R 4 π R = S ( ω ) 4 π 0 m m 1 e m 1 z z s J 0 ( m r ) d m
where r = ( x x s ) 2 + ( y y s ) 2 , m 1 = m 2 k 1 2 .

2.3. 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.
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:
2 G i + k i 2 G i = 0 ,           i = 1 , 3 , , n
2 G i + k i 2 G i = S ( ω ) δ ( R R 0 )           i = 2
where k i = ω v i 1 j 2 Q , j = 1 .
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
G i z = G i + 1 z , ρ i G i = ρ i + 1 G i + 1 , ( i = 1 , 2 , , n 1 )
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:
G i = S ( ω ) 4 π 0 C i e m i z + D i e m i z J 0 ( m r ) d m ,           i = 1 , 3 , , n
G i = S ( ω ) 4 π e i k i R R + 0 C i e m i z + D i e m i z J 0 ( m r ) d m ,             i = 2
where m i = m 2 k i 2   ,   k i 2 = ω 2 v ˜ i 2 , v ˜ i = v i 1 j 2 Q n , i = 1 , 2 , , n , j = 1 . In 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
G i = S ( ω ) 4 π 0 C i e m i ( z z i ) + D i e m i ( z z i 1 ) J 0 ( m r ) d m ,           i = 1 , 3 , , n
G i = S ( ω ) 4 π e i k i R R + 0 C i e m i ( z z i ) + D i e m i ( z z i 1 ) J 0 ( m r ) d m ,       i = 2
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.
D 2 = m m 2 H 2 d * H 3 u e m 2 z 2 z s e m 2 h 2 e m 2 z 1 z s 1 H 3 u H 2 d * e 2 m 2 h 2
C 2 = H 3 u m m 2 e m 2 z 2 z s + D 2 e m 2 h 2
D i = D i 1 m i 1 e m i 1 h i 1 1 + H i u m i 1 + H i + 1 u e 2 m i h i
C i = H i + 1 u D i e m i h i
where h i = z i + 1 z i , H 2 d = ρ 1 m 2 ρ 2 m 1 , H 2 d * = 1 H 2 d 1 + H 2 d , H n u = ρ n m n 1 ρ n 1 m n , H n u = 1 H n u 1 + H n u , H i + 2 u = 1 H i + 2 u 1 + H i + 2 u , H i + 2 u = ρ i + 2 m i + 1 m i + 2 ρ i + 1 1 H i + 3 u e 2 m i + 2 h i + 2 1 + H i + 3 u e 2 m i + 2 h i + 2 , i = 3 , , n 1 .

3. 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.

3.1. Partial Closed Form Expression

The Sommerfeld integration in Equations (11) and (12) can be written in the following form:
I = 0 m m i g ( m i ) e m i z z s J 0 ( m r ) d m = 0 Γ m J 0 ( m r ) d m
The integration is then divided into two parts,
0 Γ m J 0 ( m r ) d m = 0 p Γ m J 0 ( m r ) d m + p Γ m J 0 ( m r ) d m
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
p Γ m J 0 ( m r ) d m 0 Γ m J 0 ( m r ) d m 0 p Γ m J 0 ( m r ) d m
On substituting (19) in (18), we get
I 0 p Γ m Γ m J 0 ( m r ) d m + 0 Γ m J 0 ( m r ) d m
The kernel function g ( m i ) in (17) can be approximated by an exponential function,
Γ m = g ( m i ) m m i e m i z z = l = 1 p b a ( l ) exp b ( l ) m i m m i e m i z z
where p b 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,
0 Γ m J 0 ( m r ) d m = l = 1 p b a   ( l ) exp ( j k i R l ) R l
where R l = r 2 + b   ( l ) z z s 2 ,
r = ( x x s ) 2 + ( y y s ) 2 , ( x s , y s , z s ) is the location of the source.
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.

3.2. ESPRIT Algorithm

The signal g ( m i ) can be sampled as [8]
m i = k i T 02 + T 01 T 02 T 01 t , 0 t T 01
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 i 2 + k i 2 , the first sampling in m-plane is k i 1 + T 02 2 . Therefore, the approximation of g ( m i ) starts from m = k i 1 + T 02 2 , and the parameter p in formula (20) should be set not less than k i 1 + T 02 2 to insure the integration accuracy.
Then the sampling sequence can be expressed as
g ( m i ) = y ( t ) l = 1 p b A   ( l ) exp B   ( l ) t
The relation of A ( l ) , B ( l ) , and a   ( l ) , b   ( l ) can be obtained from Equations (23) and (24).
l = 1 p b a   ( l ) exp b   ( l ) m i = l = 1 p b A ( l ) exp B ( l ) t
Then, the unknown coefficients a   ( l ) , b   ( l ) in (24) are obtained.
b   ( l ) = B   ( l ) T 01 k i ( T 01 T 02 )
a   ( l ) = A   ( l ) exp b   ( l ) k i T 02
For the sampling sequence y ( 0 ) , y ( 1 ) , , y ( N 1 ) , a data matrix can be constructed:
Y = y ( 0 ) y ( 1 ) y ( L ) y ( 1 ) y ( 2 ) y ( L + 1 ) y ( N L 1 ) y ( N L ) y ( N 1 ) ( N L ) × ( L + 1 )
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,
Y = U Σ V H = U s U n Σ s 0 0 Σ n V s H V n H
where U is ( N L ) × ( N L ) orthogonal matrix, and V is ( L + 1 ) × ( L + 1 ) orthogonal matrix, Σ is ( N L ) × ( L + 1 ) diagonal matrix with main diagonal element σ l , which is the singular value of matrix Y . For signals without noise, Y has pb non-zero singular values σ l (l= 1,2,..., pb), and pb represents the highest order of the exponential signal of formula (24). If the signal contains noise, mode number pb can be recorded by setting a minimum threshold for σ l .
Take out the first pb dominant right singular vectors in V s matrix to form ( L + 1 ) × p b matrix V s p b .The last line of V s p b is deleted to obtain L × M matrix V 1 ; the first line of V s p b is deleted to get L × M matrix V 2 . Construct a matrix Ψ [23]
Ψ = V 1 H V 1 1 V 1 H V 2
Find the eigenvalues λ l of the matrix Ψ ,
B   ( l ) = log λ   l ,   l = 1 , 2 , , p b
For N sampled signals,
λ = 1 1 1 λ 1 λ 2 λ p b λ 1 N 1 λ 2 N 1 λ p b N 1 , Y = y ( 0 ) y ( 1 ) y ( N 1 ) , A = A ( 1 ) A ( 2 ) A ( p b )
and
Y = λ A
According to the least square method, we can obtain
A = λ H λ 1 λ H Y
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.

3.3. 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:
I = 1 1 f f ( ξ ) d ξ
Let a variable transform:
ξ = φ ( t )   and   φ ( ) = 1 , φ ( + ) = 1
be applied into (35) so as to change the interval [ 1 , 1 ] into the infinite interval [ , + ]
I = + f f ( φ ( t ) ) φ ( t ) d t
The DE rules are transformed by the tanh–sinh formula,
φ t = tanh g s t = tanh sinh t
φ ( t ) = g s ( t ) sech 2 g s ( t ) = cosh ( t ) cosh 2 ( sinh ( t ) )
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
I = h n = f f ( φ ( n h ) ) φ ( n h ) h n = N N ω n f f ( ξ n )
with the nodes ξ k and weights ω k defined as
ξ n = 1 δ n , ω n = 2 g ( n h ) δ n ( 1 q n ) 1
where
δ n = 2 q n ( 1 + q n ) 1 , q n = e 2 g s ( n h )
For arbitrary integral interval a , b may be mapped onto 1 , 1 by the linear transformation ξ = σ x + γ with σ = ( a b ) / 2 , γ = ( a + b ) / 2 , which leads to
a b f f ( ξ ) d ξ = σ 1 1 f f ( σ x + γ ) d x
Hence, for an arbitrary interval a , b , the nodes and weights become σ ξ k + γ and σ ω k , respectively, and (43) is transformed to
a b f f ( ξ ) d ξ σ h g s ( 0 ) f f ( γ ) + n = 1 N ω n f f ( a + σ δ n ) + f f ( b σ δ n )
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.
0 p Γ m Γ m J 0 ( m r ) d m = 0 b k + b k p Γ m Γ m J 0 ( m r ) d m
where breakpoint bk is set to the real part of wavenumber k i in this paper.

4. 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].

4.1. 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/cm3, v 2 = 2000 m/s, ρ 2 = 1.5 g/cm3. The integration path is divided into three segments, the breakpoint b k = 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.
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/cm3.

4.2. 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/cm3, v 2 = 2000 m/s, ρ 2 = 2.0 g/cm3, v 3 = 3000 m/s, ρ 3 = 3.0 g/cm3. 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).
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.

5. 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.

Author Contributions

Conceptualization, Z.Z. and S.D.; methodology, S.L.; software, S.L.; validation, S.L., Y.Y. and Z.Z.; formal analysis, S.L.; writing—original draft preparation, S.L.; writing—review and editing, I.I. and Y.Y.; supervision, Z.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Fundamental Research Funds for the Central Universities of Central South University 2019ZZTS300.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author. The data are not publicly available due to some institutional reasons.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Fuchs, K.; Müller, G. Computation of Synthetic Seismograms with the Reflectivity Method and Comparison with Observations. Geophys. J. R. Astron. Soc. 1971, 23, 417–433. [Google Scholar] [CrossRef] [Green Version]
  2. Yao, Z.; Zheng, T. A Generalized Reflection-transmission Coefficient Matrix and Discrete Wavenumber Method for Synthetic Seismograms (ii)- Multiple sources at Different Depths. Chin. J. Geophys. 1984, 27, 338–348. [Google Scholar]
  3. Apsel, R.J. Dynamic Green’s functions for layered media and applications to boundary-value problems. Aviat. Week Space Technol. 1979, 125, 209–211. [Google Scholar] [CrossRef] [Green Version]
  4. Bouchon, M. A simple method to calculate Green’s functions for elastic layered media. Bull. Seism. Soc. Am. 1981, 71, 959–971. [Google Scholar] [CrossRef]
  5. Olson, A.H. Forward Simulation and Linear Inversion of Earthquake Ground Motions; University of California: San Diego, CA, USA, 1982. [Google Scholar]
  6. Yao, Z. A generalized reflection-transmission coefficient matrix and discrete wavenumber method for synthetic seismograms. Bull. Seism. Soc. Am. 1983, 73, 647–654. [Google Scholar]
  7. Karabulut, E.P.; Aksun, M.I. A Novel approach for the Efficient and Accurate Computation of Sommerfeld Integral Tails. In Proceedings of the 2015 International Conference on Electromagnetics in Advanced Applications (ICEAA), Torino, Italy, 7–11 September 2015. [Google Scholar]
  8. Guo, L. Research and Application of Fast Integral Equation Method Based on Layered Medium Green’s Function. Ph.D. Thesis, School of Electronic Engineering of UESTC, Chengdu, China, 2016. [Google Scholar]
  9. Press, W.H.; Flannery, B.P.; Teukolsky, S.A.; Vetterling, W.T.; Chipperfield, J.R. Numerical recipes: The art of scientific computing: Cambridge University Press, Cambridge, 1986 (ISBN 0-521-30811-9). xx + 818 pp. Price 25.00. Anal. Chim. Acta. 1987, 199, 293–294. [Google Scholar] [CrossRef]
  10. Golubovic, R.; Polimeridis, A.G.; Mosig, J.R. Efficient Algorithms for Computing Sommerfeld Integral Tails. IEEE Trans. Antennas Propag. 2012, 60, 2409–2417. [Google Scholar] [CrossRef] [Green Version]
  11. Niciforovic, R.G.; Polimeridis, A.; Mosig, J.R. Fast Computation of Sommerfeld Integral Tails via Direct Integration Based on Double Exponential-Type Quadrature Formulas. IEEE Trans. Antennas Propag. 2011, 59, 694–699. [Google Scholar] [CrossRef] [Green Version]
  12. Michalski, K.A.; Mosig, J.R. Efficient computation of Sommerfeld integral tails—Methods and algorithms. J. Electromagn. Waves Appl. 2016, 30, 281–317. [Google Scholar] [CrossRef]
  13. Mosig, J. The Weighted Averages Algorithm Revisited. IEEE Trans. Antennas Propag. 2012, 60, 2011–2018. [Google Scholar] [CrossRef]
  14. Michalski, K.A. Extrapolation methods for Sommerfeld integral tails. IEEE Trans. Antennas Propag. 1998, 46, 1405–1418. [Google Scholar] [CrossRef]
  15. Mosig, J.R.; Gardiol, F.E. A Dynamical Radiation Model for Microstrip Structures. Adv. Electron. Electron Phys. 1982, 59, 139–237. [Google Scholar]
  16. Takahasi, H.; Mori, M. Double Exponential Formulas for Numerical Integration. Publ. Res. Inst. Math. Sci. 1973, 9, 721–741. [Google Scholar] [CrossRef] [Green Version]
  17. Rao, Y.; Wang, Y. Fracture effects in seismic attenuation images reconstructed by waveform tomography. Geophysics. 2009, 74, 25. [Google Scholar] [CrossRef]
  18. Malinowski, M.; Operto, S.; Ribodetti, A. High-resolution seismic attenuation imaging from wide-aperture onshore data by visco-acoustic frequency-domain full-waveform inversion. Geophys. J. Int. 2011, 186, 1179–1204. [Google Scholar] [CrossRef] [Green Version]
  19. Carpio, A.; Rapún, M.L. Multifrequency Topological Derivative Approach to Inverse Scattering Problems in Attenuating Media. Symmetry 2021, 13, 1702. [Google Scholar] [CrossRef]
  20. Li, Q. The Way to Precise Exploration: Systematic Engineering Analysis of High Resolution Seismic Exploration; Petroleum Industry Press: Beijing, China, 1994. [Google Scholar]
  21. Touhei, T. A scattering problem by means of the spectral representation of Green’s function for a layered acoustic half-space. Comput. Mech. 2000, 25, 477–488. [Google Scholar] [CrossRef]
  22. Faisal, K.M. Improved Matrix Pencil Methods for Parameters Estimation of Plane Wave Signals. Ph.D. Thesis, Department of Electrical Engineering, Pakistan Institute of Engineering & Applied Sciences (PIEAS), Nilore, Islamabad, Pakistan, 2011. [Google Scholar]
  23. Zhang, X. Modern Signal Processing, 2nd ed.; Tsinghua University Press: Beijing, China, 2002. [Google Scholar]
  24. Xu, K.; Wang, M. Parameter inversion of acoustic wave equation in frequency domain with Finite element method. Chin. J. Geophys. 2001, 44, 852–864. [Google Scholar] [CrossRef]
  25. Hua, Y.; Sarkar, T.K. Generalized pencil-of-function method for extracting poles of an EM system from its transient response. IEEE Trans Antennas Propagation. 1989, 37, 229–234. [Google Scholar] [CrossRef] [Green Version]
  26. Alparslan, A.; Aksun, M.I.; Michalski, K.A. Closed-Form Green’s Functions in Planar Layered Media for All Ranges and Materials. IEEE Trans. Microw. Theory Tech. 2010, 58, 602–613. [Google Scholar] [CrossRef]
  27. Boix, R.R.; Fructos, A.L.; Mesa, F. Closed-Form Uniform Asymptotic Expansions of Green’s Functions in Layered Media. IEEE Trans. Antennas Propag. 2010, 58, 2934–2945. [Google Scholar] [CrossRef]
Figure 1. n layers structure of homogeneous medium.
Figure 1. n layers structure of homogeneous medium.
Symmetry 13 01969 g001
Figure 2. Comparison of ESPRIT and GPOF results with original data. (a) Fitting results; (b) fitting error.
Figure 2. Comparison of ESPRIT and GPOF results with original data. (a) Fitting results; (b) fitting error.
Symmetry 13 01969 g002
Figure 3. Comparison of DE_DCIM and DE_WA solution in half-space. (a) real part of DE_DCIM solution; (b) real part of DE_ WA solution; (c) real part of relative error between (a,b). (d) Imaginary part of DE_DCIM solution; (e) imaginary part of DE_ WA solution; (f) imaginary part of relative error between (d,e).
Figure 3. Comparison of DE_DCIM and DE_WA solution in half-space. (a) real part of DE_DCIM solution; (b) real part of DE_ WA solution; (c) real part of relative error between (a,b). (d) Imaginary part of DE_DCIM solution; (e) imaginary part of DE_ WA solution; (f) imaginary part of relative error between (d,e).
Symmetry 13 01969 g003
Figure 4. Comparison of the wavefield of 3-layer model (DE_DCIM method and FEM method). (a) Real part of DE_DCIM solution; (b) real part of FEM solution; (c) real part of relative error between (a,b). (d) Imaginary part of DE_DCIM solution; (e) imaginary part of FEM solution; (f) imaginary part of relative error between (d,e).
Figure 4. Comparison of the wavefield of 3-layer model (DE_DCIM method and FEM method). (a) Real part of DE_DCIM solution; (b) real part of FEM solution; (c) real part of relative error between (a,b). (d) Imaginary part of DE_DCIM solution; (e) imaginary part of FEM solution; (f) imaginary part of relative error between (d,e).
Symmetry 13 01969 g004
Figure 5. The results of log 10 G and relative errors in line ( x , y s , z s ) 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.
Symmetry 13 01969 g005
Table 1. Computation time comparison with different methods.
Table 1. Computation time comparison with different methods.
MethodGPOFESPRIT
Computation time0.135 s0.025 s
Table 2. Computation time comparison with different parameters (computation of ( x , y s , z ) plane).
Table 2. Computation time comparison with different parameters (computation of ( x , y s , z ) plane).
Method\Parameters v 1   =   340   m / s
v 2   =   2000   m / s
v 1   =   1000   m / s
v 2   =   2000   m / s
v 1   =   2000   m / s
v 2   =   3000   m / s
DE_WA91.2 s67.8 s57.6 s
DE_DCIM50.8 s38.9 s32.5 s
Time saving44.3%42.6%43.6%
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Liu, S.; Zhou, Z.; Dai, S.; Iqbal, I.; Yang, Y. Fast Computation of Green Function for Layered Seismic Field via Discrete Complex Image Method and Double Exponential Rules. Symmetry 2021, 13, 1969. https://doi.org/10.3390/sym13101969

AMA Style

Liu S, Zhou Z, Dai S, Iqbal I, Yang Y. Fast Computation of Green Function for Layered Seismic Field via Discrete Complex Image Method and Double Exponential Rules. Symmetry. 2021; 13(10):1969. https://doi.org/10.3390/sym13101969

Chicago/Turabian Style

Liu, Siqin, Zhusheng Zhou, Shikun Dai, Ibrar Iqbal, and Yang Yang. 2021. "Fast Computation of Green Function for Layered Seismic Field via Discrete Complex Image Method and Double Exponential Rules" Symmetry 13, no. 10: 1969. https://doi.org/10.3390/sym13101969

APA Style

Liu, S., Zhou, Z., Dai, S., Iqbal, I., & Yang, Y. (2021). Fast Computation of Green Function for Layered Seismic Field via Discrete Complex Image Method and Double Exponential Rules. Symmetry, 13(10), 1969. https://doi.org/10.3390/sym13101969

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop