(Department of Earth Sciences, Petroleum University, Beijing 100083)
Teng Zhang Ji Wen intermediary
(Institute of Geophysics, Chinese Academy of Sciences, Beijing10010/)
In exploration seismology, it is of great significance to simulate the propagation of seismic waves in anisotropic media quickly and effectively. The stability analysis of the algorithm is very necessary to speed up the calculation. In this paper, firstly, the wave equation is described by matrix and vector, and a fast finite difference method with less memory consumption is proposed for the elastic wave equation in two-dimensional and three-dimensional general anisotropic media. Then, the stability conditions of finite difference schemes for two-dimensional homogeneous and inhomogeneous anisotropic wave equations are systematically studied. In addition, the stability formula of finite difference method under some special anisotropic conditions is given. Finally, the stability of the three-dimensional finite difference scheme is also studied.
Elastic wave equation; Finite difference; Stable anisotropic medium
1 Introduction
Numerical simulation of seismic wave propagation is of great significance in earth science. Among various methods of anisotropic earthquake simulation, based on Kennett's research work. Kosloff et al [13] used Fourier method to simulate seismic anisotropy. Chen [7] uses finite element method to simulate non-uniform anisotropy. Although the finite difference method has been widely used to simulate elastic waves in isotropic media, it is not common to simulate seismic anisotropy by this method. Tsingas et al. [16] developed a simulation algorithm for solving partial differential equations in transversely isotropic media by using finite difference operator, which is based on MacCormack-type separation scheme [2]. Faria et al. [8] based on staggered grid scheme [17], used finite difference algorithm to simulate two-dimensional transversely isotropic problem. Recently, Igel et al. [9] presented a finite difference algorithm based on convolution algorithm to simulate general seismic anisotropy.
Fast speed and less storage are the advantages of finite difference method. With the need of large-scale wave field simulation and the development of large-scale parallel computing, finite difference seismic simulation of complex media or high-dimensional (two-dimensional and three-dimensional) models is inevitable. The virtual spectrum method is very popular because its spatial operators accurately reach the Nyquist frequency. But the virtual spectrum method needs Fourier transform, so it is very time-consuming to simulate three-dimensional anisotropy. At the same time, adopting Fourier transform means that each grid point interacts with other points. In a sense, this is inconsistent with the local elastic properties of dynamics. Therefore, it is necessary to consider the locality of difference operator when designing the finite difference scheme for earthquake simulation. On the other hand, it is very important to consider the locality of difference operator to improve the parallelism of the algorithm. Because the information exchange between the nearest grid points is the fastest, it is feasible to simulate the wave field of anisotropic large-scale model.
Based on the above reasons, this paper presents a fast finite difference method for seismic wave propagation in anisotropic media with less memory occupation. In fact, this is a generalization of the algorithm [10, 18].
The numerical simulation of seismic wave propagation usually adopts time marching algorithm. In order to ensure the stability of the algorithm, the time increment is limited by the stability conditions of the algorithm. Choosing a suitable time step can not only ensure the stability of numerical calculation, but also speed up the calculation. Otherwise, it will not only produce non-physical numerical oscillation, but even lead to wrong results. The choice of reasonable time increment depends on the difference scheme and the medium parameters or speed describing the characteristics of the medium, in other words, on the stability conditions of the difference scheme. Therefore, the stability of finite difference scheme is very important in numerical calculation. Although some special cases have been studied in the literature [10, 18], they have not made a detailed and systematic study of the finite difference scheme and its stability for general anisotropic problems. Our purpose is to give a general finite difference scheme for this problem and its stability conditions, and further give some stability formulas under special anisotropic conditions. The results show that our results are consistent with those of Aboudi[ 1] under isotropic conditions.
2 finite difference scheme
2. 1 three-dimensional anisotropy
The equation of motion equilibrium in anisotropic media can be in the following form:
Lithospheric structure and deep action
Where: ρ is the medium density; Fx, fy and fz respectively represent the components of the force source in X, Y and Z directions; Ux, uy and uz represent displacement components in x, y and z directions, respectively.
The stress-strain relationship is, where the elastic parameter matrix and ci, j=cj, I, I, j= 1, 2, …, 6;
Stress matrix σ = (σ xx, σyy, σzz, σyz, σxz, σ xy) t;
Strain matrix ε = (ε xx, εyy, εzz, εyz, εxz, ε xy) t;
The relationship between stress and displacement is:
Lithospheric structure and deep action
manufacture
Lithospheric structure and deep action
Lithospheric structure and deep action
Then equation (1) can be written as:
Lithospheric structure and deep action
Obviously, A, E, Q, B+D, C+G and F+H are real symmetric matrices. Without considering the source term, the following finite difference scheme can be obtained by using the finite difference approximate equation (2):
Lithospheric structure and deep action
Where Δ x, Δ y and Δ z respectively represent the space step in X, Y and Z directions, and Δ t represents the time step.
Lithospheric structure and deep action
2.2 two-dimensional anisotropy
Similarly, the wave equation in two-dimensional anisotropic media can be written as:
Lithospheric structure and deep action
And there are the following different formats:
Lithospheric structure and deep action
3 stability conditions
3. Stability conditions of two-dimensional finite difference schemes in1homogeneous media
Format (5) can be simplified to:
Lithospheric structure and deep action
According to the stability theory analysis of Richtmyer and Morton, we make
Lithospheric structure and deep action
Where U0 is a constant vector, equation (6) becomes:]]
2I—2Pλ+I)U0=0
Where I stands for identity matrix,
Lithospheric structure and deep action
In which:,,,
If the coefficient determinant is 0, then:
Lithospheric structure and deep action
Theorem 1 difference scheme (6) is stable under the following conditions:
Lithospheric structure and deep action
Where the functions F(α, β) and α, β are:
Lithospheric structure and deep action
Lithospheric structure and deep action
Where k = δ z/δ x
According to the stability of difference scheme (6), it is proved that λ in equation (7) satisfies ‖λ‖≤ 1.
According to (7) and Lemma 2 (see Appendix), there are the following inequalities:
Lithospheric structure and deep action
According to the symmetry of A, Q and C+G, the matrix is also symmetric, so from Lemma 3 (see Appendix) and inequality (8), there are
Lithospheric structure and deep action
From the nonnegativity of the elements of matrices A, Q and C+G, we can make 0≤α, and to make (9) hold, we only need to
Lithospheric structure and deep action
manufacture
Lithospheric structure and deep action
Find partial derivative
Lithospheric structure and deep action
According to the characteristics of wave equation, there are
‖A‖0，Q‖0，C+G‖0
So (0,0) and (,) may be the extreme points of the function. Obviously,
f(0，0)=0
Lithospheric structure and deep action
Let's discuss the case of (α, β) ≦ (0,0) or (,):
Through (1 1) and (12), if, then
Lithospheric structure and deep action
Where α and β may be extreme points of f(α, β), so the stability condition is:
Lithospheric structure and deep action
Otherwise, it is the maximum point of the function, and the stability condition is:
Lithospheric structure and deep action
Especially, when δ x ~ δ z, there are the following simplified stability conditions:
Lithospheric structure and deep action
Where α and β are determined by (13) and (14), and k= 1.
3.2 Stability conditions under non-uniform conditions
In the case of non-uniformity, it is difficult to determine the stability condition of the difference scheme (5). However, according to the numerical method theory of partial differential equations [15], we can use the "freezing coefficient" method [14] to analyze its stability. In addition, we give the stability conditions of the difference scheme (5) in the case of inhomogeneous media.
In fact, if the parameter function of the medium is continuous and bounded, we can generally regard a small calculation area as consistent, and then the difference scheme (5) can be degenerated into scheme (6), so that the stability conditions in the small area are also consistent. In addition, according to the continuous boundedness of the medium parameter function, we can obtain the stability conditions of the difference scheme (5) as follows:
If so.
Lithospheric structure and deep action
Otherwise,
Lithospheric structure and deep action
Where H(α, β), α and β are:
Lithospheric structure and deep action
3.3 Stability conditions of difference schemes in some special media
Obviously, by calculating the norms of matrices A, Q and C+G in format (6), the stability conditions in isotropic and transversely isotropic cases can be obtained as follows:
isotropy
Lithospheric structure and deep action
or
Lithospheric structure and deep action
Where λ and μ are Lame constants and P-wave velocities.
Obviously, the stability condition (15) is consistent with the result of Aboudi[ 1].
Transverse homogeneity
If it is ≤p, otherwise.
Lithospheric structure and deep action
Among them,
Lithospheric structure and deep action
Where a, n, l, f and c are elastic constants.
Similarly, we can obtain the stability conditions of inhomogeneous and other special anisotropic media (such as cubic anisotropy, orthogonal anisotropy, etc.).
4 Stability conditions of difference schemes in the case of three-dimensional anisotropy
The following stability conditions can be obtained through the analysis in the previous two-dimensional case:
Theorem 2 If
Max [K 1 ‖ A ‖+K2 ‖ E ‖+K2 ‖ Q ‖, f2(θ2, θ2, θ3)]≤ 1, then the difference scheme (3 where the function F2 (.
Lithospheric structure and deep action
Where:,,,,,;
θ 1, θ2 and θ3 satisfy:
Lithospheric structure and deep action
Where: a = 4k1‖ a ‖ b = 4k2 ‖ e ‖ d = 4k3 ‖ q ‖ g = 4k4 ‖ b+d ‖ e =
For the case of three-dimensional heterogeneous media, the stability condition can be analyzed similarly by the "freezing coefficient" method.
5 Summary and discussion
Finite difference method is an important tool for numerical simulation of seismic wave propagation, and its stability condition is one of the keys to improve the calculation speed. However, for general two-dimensional and three-dimensional anisotropic media, there are few systematic and in-depth studies on their difference schemes and stability conditions. In this paper, a fast difference scheme with less memory occupation is given, and the stability conditions of the difference scheme under general uniform and inhomogeneous anisotropy are systematically analyzed and deduced. We believe that the results of this paper are helpful to the development of anisotropic numerical simulation and provide a theoretical basis for the wide application of finite difference method.
refer to
[1] J. Abdi. Numerical simulation of earthquake source. Geophysics, 197 1, 36,810 ~ 821.
[2]A.Bayliss, K.E.Jordan, B.J.LeMesurier and E.Turkel. Fourth-order accurate finite difference scheme for calculating elastic waves. Bull.Seis.Soc.Am,1986,76,11/.
[3]D.C.Booth and S.Crampin Anisotropic reflectivity technology: theory. Journal of Geophysics, astr SOC, 1983(a), 72,755 ~ 756.
[4]D.C.Booth and S.Crampin Anisotropic reflectivity technology: abnormal arrival waves from anisotropic upper mantle. Journal of Geophysics, American Society of Remote Sensing of the Earth, 1983(b), 72,767 ~ 782.
[5]V.Cerveny. Seismic rays and ray intensity in inhomogeneous anisotropic media. Journal of Geophysics. astr SOC,1972,29.1~13.
[6]V.Cerveny and P.Firbas. Numerical simulation and inversion of seismic body wave travel time in inhomogeneous anisotropic media. Journal of Geophysics. American Geophysical Society,1984,76,41~ 51.
Chen Guanghai. Numerical model of elastic wave propagation in anisotropic inhomogeneous media: finite element method. Extended abstract of the 54th SEG annual meeting, Houston, USA,1984,631~ 632.
[8]E.L.Faria and P.L.Stoffa. Finite difference simulation in transversely isotropic media. Geophysics,1994,59,282 ~ 289.
[9]H.Igel, P.Mora and B.Riollet. Anisotropic wave propagation through finite difference grids. Geophysics,1995,60, 1203~ 12 16.
[10]K.R.Kelly, R.W.Ward, S.Treitel and R.M.Alford. Synthetic seismogram: finite difference method. Geophysics,1976,41,2~27.
[11] b.l.n.kennett. Theoretical reflection seismogram of elastic medium. Journal of Geophysics,1979,27,301~ 321.
[12] b.l.n.kennett. Seismic wave propagation in layered media. Cambridge: Cambridge University Press, 1983.
[13]D.Kosloff and E.Baysal. Forward simulation by Fourier method. Geophysics,1989,47, 1402~ 14 12.
Lu, Gu Lizhen, Chen. Difference methods for partial differential equations. Beijing: Higher Education Press, 1988.
[15]R.D.Richtmyer and K.W.Morton. Difference method for initial value problem. New york: Interscience, 1967.
[16]C.Tsingas, A.Vafidis and E.R.Kanasewich. Calculation of elastic wave propagation in transversely isotropic media by finite difference method. Journal of Geophysics,1990,38,933 ~ 949.
[17]J.Virieus.P-SV wave propagation in heterogeneous media: velocity-stress finite difference method. Geophysics, 1986, 5 1, 889~90 1.
[18] Zhang Zhijun, He Qingde, Teng Junwei, et al .. Simulation of three-component seismic records in two-dimensional transversely isotropic media by finite difference method. Journal of Exploration Geophysics,1993,29,51~ 56.
Appendix: Lemma
Lemma 1 is a necessary and sufficient condition for the real coefficient equation λ 2-2dλ+ 1 = 0, and | d | ≤ 1 is that its root λ satisfies ∣ λ| ≤ 1.
For the proof of lemma 1, see reference [14].
Lemma 2 If A∈Rn×n and A = A', and I is identity matrix, then for the following equation.
∣λ2I—2Aλ+I∣=0，
‖A‖≤ 1 is a necessary and sufficient condition for the root λ of the equation to satisfy ∣ λ | ≤ 1
Evidence: sufficient
Because a is a real symmetric matrix, the existence of T- 1 and T∈Rn×n makes
T- 1AT = diagnosis (d 1, d2, …, dn)
According to the equation in Lemma 2 and the above equation, we can get
(λ2—2d 1λ+ 1)…(λ2—2dnλ+ 1)= 0
‖A‖≤ 1 and▎ (a) ≤ A ‖, ▎ (a) ≤1.
According to, there is
Therefore, for any λ that satisfies the equation, there is ∣∣≤1.
Necessary conditions: Since ∣λ∣≤ 1 and A are real symmetric matrices, we can obtain from lemma 1:
∣di∣≤ 1,i= 1,2,…,n
That is ρ(A)≤ 1.
And because a is a positive matrix, ρ (a) = ‖ A ‖, that is ‖A‖≤ 1.
Lemma 3 if A∈Rn×n and A = A', then
(i) If ‖ i-a ‖≤ 1, ‖ a ‖≤ 2;
(ii) If ‖A‖≥0, then ‖A‖≤2 is a necessary and sufficient condition for ‖ i-a ‖≤ 1
It is proved that: (i) Because A is a real symmetric matrix, there are T- 1 and T∈Rn×n such that T- 1AT=diag(d 1, d2, ..., dn) ≡ d.
According to ‖ I-A ‖≤ 1, there is ‖ T (I-D) T- 1 ‖≤ 1.
Obviously, I-D is a normal matrix, so ‖ t (I-d) t-1‖ = ‖ I-d ‖ ≤1.
that is
So, that is ‖D‖≤2.
Because a is a normal matrix, so
‖A‖= TDT- 1‖= D‖, that is, ‖A‖≤2.
(ii) The necessary condition has been proved in (i), and the sufficient condition is proved below.
According to the proof process of (i):
‖A‖= D‖
Because ‖A‖≤2 and ‖A‖≥0, there is 0 ‖ D ‖≤ 2, that is, the maximum ∣ Di | ≤ 2.
So we can.
That is, I-d ‖ = t-1(I-a) t ‖≤1.
We can get ‖ I-A ‖≤ 1 from the symmetry of a.