Eulerian Gaussian beams for high frequency wave propagation in inhomogeneous media of arbitrary anisotropy
1 Eulerian Gaussian beams for high frequency wave propagation in inhomogeneous media of arbitrary anisotropy
Xingguo Huang University of Bergen, Department of Earth Sciences, Allegaten 41, 5020 Bergen, Norway. Email: [email protected].
Abstract
We present the Eulerian Gaussian beam method in anisotropic media. We derive kinematic and dynamic ray tracing equations based on the level set theory and Eulerian theory using the anisotropic eikonal equation. Compared with the traditional anisotropic Gaussian beam method using ray-centered coordinates, the anisotropic Eulerian Gaussian beam method derived in this work has the following three advantages: (1) it can handle the problem of calculating the distance from the imaging point to the beam point more easily; (2) it allows the travel time and amplitude to be distributed uniformly within the actual computational domain without interpolation; (3) it can handle late arrivals, both theoretically and in calculations, due entirely to ray tracing in the phase space.
Keywords:
Eulerian Gaussian beam; level set method; seismic anisotropy; wave propagation
1. Introduction
The Gaussian beam method is an extension of ray theory. As a result of its use of complex travel times and amplitudes, the Gaussian beam method can deal with the problems that the high-frequency ray theory faces with respect to singularities. The Gaussian beam method has been used for high-frequency wave field calculations (Červený, 2001; Červený et al., 1982; Popov, 1982) as well as for migration imaging (Hill, 1990, 2001; Gray, 2005; Gray and Bleistein, 2009). Much work has been carried out to improve the Gaussian beam method (Alkhalifah, 1995; Červený, 1985; George et al., 1987; Klimes, 1984; Kravtsov and Berczynski, 2007; Nowack, 2003; Nowack and Kainkaryam 2011; Popov et al., 2010; Tanushev, 2008; Zhu et al., 2007). However, the beam is calculated using the ray-centered coordinate system and it is necessary to calculate the distance of the imaging point to the nearest point on the ray or to carry out coordinate transformation (Červený and Pšenčik, 2010; Klimes, 1994). This is a very time-consuming calculation. In contrast, the multi-value travel time is automatically included in the Gaussian beam method based on complex travel times, although a quantitative theoretical proof on late arrivals processed by the Gaussian beam method has not yet been reported. For these reasons, Leung et al. (2007) proposed the Eulerian Gaussian beam method. The Eulerian Gaussian beam method is a high-frequency asymptotic expression of the wave field based on level set theory (Qian and Leung, 2004, 2006; Leung et al., 2004) and is established in the Cartesian coordinate system. To calculate the Eulerian Gaussian beam, ray tracing and dynamic ray tracing are performed and the partial differential equations are solved. However, unlike traditional dynamic ray tracing in the Gaussian beam method using ray-centered coordinates, ray tracing and dynamic ray tracing in the Eulerian Gaussian beam method is based on the level set theory and the eikonal equation in Eulerian theory. Ray tracing of a single ray in the coordinate system is conducted point by point, so there are additional difficulties resulting from the uniform distribution of the phase and amplitude in the entire space and interpolation is required to determine the amplitude. However, due that Eulerian Gaussian beam theory is developed in the Cartesian coordinate system, the calculation is performed on regular grid points. Thus, the calculated traveltime and amplitude are distributed uniformly. In geophysics, to overcome the shortcomings of the ray method in calculating the travel time, eikonal calculations based on the finite difference method (Vidale, 1990; van Trier and Symes, 1991; Sethian, 1996) are widely used. The method is based on the Eulerian theory, but only to calculate the first arrival – that is, the minimum phase – not to calculate later arrivals. Calculation at the level set theory has therefore been used in this work. The aim of the work reported here was to obtain the Eulerian Gaussian beam in anisotropic media. Based on the anisotropic eikonal equation, we determined the derivative formula for ray tracing and then substituted it into the kinetic ray tracing equations to obtain the anisotropic dynamic ray tracing equations. We also modified the phase formula and obtained the anisotropic travel time formula. Finally, we derived weighting functions and superposition formulas for the Eulerian Gaussian beam based on the steepest descent and the stationary phase theory in anisotropic media. As a result of space constraints, this discussion is limited to the theory and there is no discussion and presentation of the numerical methods. We review elastodynamic theory and some basic concepts of high-frequency wave fields. And then we deduce the basic formulas required in ray tracing from the anisotropic eikonal equation to obtain the ray tracing and dynamic ray tracing equations for anisotropic Eulerian Gaussian beam. Finally, we determine the anisotropic Eulerian Gaussian beam and Green's function by using Eulerian Gaussian beam superposition.
2. Basic concepts and equations
The elastic dynamic equation in anisotropic media (Červený, 2001) is as follows (1) where is the wave field, are the Cartesian components of the displacement vector is the elasticity modulus and r is the density. The high-frequency seismic wave field can be expressed as (Červený, 2001) (2) where represents the amplitude. To find the amplitude of the ray, it is necessary to know the initial amplitude at the starting point and the corresponding geometrical spreading factor. T is the travel time of the seismic wave. Substituting this into Equation (1), we obtain the following equation (Červený and Pšenčik, 2010) (3) The Christoffel matrix can be expressed as (Červený and Pšenčik, 2010) (4)
In this formula, , are the slowness parameters and the elastic modulus is (Červený and Pšenčik , 2010) ( ) , , u ijkl k l ij c u r = !! i u , k l u ijkl c ( ) ( ) ( ) x , x exp{ [ x ]} i i i u t U i t T w = - - U ( ) ik ik k U dG - = ( ) , ik m n x p G ( ) x , ik m n ijkl j l p a p p G = j p l p (5)
3. Paraxial anisotropic eikonal equation
Considering the Hamiltonian form of the eikonal equation in two-dimensional anisotropic media (Qian and Leung, 2004): (6) where the slowness parameters , are equal to: (7) is the phase angle and is the wave velocity. To obtain the ray tracing system, we give the relationship between dx,dz and . According to Qian and Leung’s (2004) research, we obtain (8) (9)
On the other hand, to obtain , Qian and Leung (2004) gave the relationship between and as (10) By obtaining the derivation of equation (7), we can also obtain the slowness parameter derivation directly: ( ) ( )( ) ijkl nijkl n n c xa x x r = ( , , , ) 0 F x z p p = p p sin( , , )cos( , , ) p v x zp v x z qqqq == q v t
11 31 3 1 ( ) dx F F Fp pd p p p t - ¶ ¶ ¶ = + ¶ ¶ ¶
11 31 3 3 ( ) dz F F Fp pd p p p t - ¶ ¶ ¶ = + ¶ ¶ ¶ dd qt , dp dp t
11 1 31 3 13 1 31 3 ( )( ) dp F F Fp pd p p xdp F F Fp pd p p z tt -- ¶ ¶ ¶ = + ¶ ¶ ¶¶ ¶ ¶ = + ¶ ¶ ¶ (11) Combining equations (10) with (11), we obtain: (12)
When the rays are nearly level, they are required to meet the following conditions (Symes and Qian, 2003): (13)
At this time (14) (15)
Equations (14) and (15) are used here as the anisotropic ray tracing and dynamic ray tracing equations in next section.
4. Ray tracing and dynamic ray tracing equations
In this work, ray tracing for the Eulerian Gaussian beam has been established based on level set theory. Generally, the level set equation has the form (Qian and Leung, 2004) (16) in the equations (17) (18)
The derivatives (17) and (18) are Equations (14) and (15), respectively where above. According to Equation (16), we can obtain the partial differential equations with ray tracing cos sin sin ( )sin cos cos ( ) vvdp d v dx v dzd v d v x d z dvvdp d v dx v dzd v d v x d z d q q q qqt t t tq q q qqt t t t¶- ¶ ¶¶ = - + ¶ ¶¶- ¶ ¶¶ = - + ¶ ¶
11 31 3 ( ) ( cos sin) d F F F Fp p v vd p p x z qt - ¶ ¶ ¶ ¶ = + -¶ ¶ ¶ ¶
11 31 3 3 ( ) 0 dz F F Fp pd p p p t - ¶ ¶ ¶ = + > ¶ ¶ ¶
11 3 ( ) dx F Fdz p p - ¶ ¶ = ¶ ¶ ( ) ( cos sin ) d F F Fv vdz p x z q q - ¶ ¶ ¶ = - ¶ ¶ ¶ z x d u wdz q f f f f = + + = dxu dz = dw dz q = as follows: (19) The initial conditions for numerical solution (Leung et al., 2007) are: (20)
The dynamic ray tracing equations (Leung et al., 2007) are: (21) where (22)
For isotropic medium , For anisotropic media . Thus (23)
Similar to isotropic media, when the ray tracing equations are solved numerically, the initial value is amended as follows: z xz xz x dx x ux wxdzd u wdzd dxu w p pdz dz q qq J J J Jt t t t = + + == + + == + + = + (0, , )(0, , )(0, , ) 0 x x xxx qJ q qt q === z x xp xxz x pp xp
B uB wB H B H CC uC wC H B H C qq + + = - - + + = - + (v 3 ) pp Txxp T Txx x x x xxx H vuu v NH vu v v v u v vH v == - += u q = cos sin dxu dz q q = +
22 33 (cos sin )(cos sin )(cos sin )( 3 ) (cos sin ) pp Txxp T Txx x x x xxx dxH v dzdx v NdzH v dx dxvv v v v vdz dzH v q qq qq q q q = ++= + - + += (24)
5. Eulerian Gaussian beam
For the Gaussian beam method, the complex travel time is treated with a second-order Taylor expansion (Leung et al, 2007): (25)
The Hessian matrix is , where is the travel time. The value of is obtained by solving the dynamic ray tracing equations, as described later in this paper; the direction of the central beam is : The amplitude value of an arbitrary point on the ray can be expressed as follows (Leung et al., 2007): (26)
Substituting Equations (25) and (26) into Equation (2), the Eulerian Gaussian beam in two-dimensional isotropic media has the form (Leung et al., 2007): (27)
We can obtain the following equation for the travel time in anisotropic media: (28) (0, , ) cos sin(0, , ) cos sin dxB x i dzdxC x dz q q qq q q = ' += + [ ] T Ts x N x X x X M x Xvu t t + - + - - ( z, x) = ( z, M BC - = t , B C N u N = + [ ][ ] (z, x) (0, ) det (0)(z, ) (0, ) (z, ) det (z, x, ) ss s s v t x CA x v x t x C q = [ ]
01 20 ( , ) sinexp{ [ ( , , ) ( )(0, )det ( , , ) cos1 ( , , ) ( , , )]( ) }2 beam s v z xu i z x x xv x C z x vB z x C z x x x qwt qq qq q - = + - + - sin sin ( ) 2 cos sin(z, x) ( , , ) ( )(cos sin )1 ( , , ) ( , , )( )2 dx dxdz dzz x x xdxv dzB z x C z x x x q q q qt t q q qq q - é ù + -ê úë û = + - ++ - In the high-frequency approximation, by using the steepest descent method for stationary phase analysis (Bleistein, 1984), we obtain the first and the second derivative of travel time as follows: (29) (30)
By using the stationary phase method, we obtain: (31) where are the eigenvectors corresponding the Christoffel matrix [Equation (4)]. The Green's function based on the general theory of high-frequency radiation in anisotropic media is (Červený, 2001): (32)
We compared Equations (32) with (31) to obtain the weighting function in Eulerian Gaussian beam superposition: (33)
Therefore the Green’s function based on Eulerian Gaussian beam superposition can be written as: (34) where the integration region is:
0, ( )/ s N v x z
T x = - =
0, 4( )/ ( ) s NN v x z dxdzx rT v a q q = - + - = [ ] i NN ki k NN g g g TT g g Tdx Tdz p wwp ww q q = - é ùë ûY = - é ùë û + i k g g [ ] i G g g Tv C wp = -
0, 20 0 det 1 (cos sin )2 4 NN T dxv dz w q qp pé ùë ûY = + [ ]
20 ( , ) ( , ) 210 sin [sin ( ) 2cos sin ]( , ) (cos sin ) exp{ [ ( , , )(0, )det ( , , ) (cos sin )1( ) ( , , ) ( , , ) ]}02 ( ) s x z s dx dxv z x dx dz dzG d i z x dxx v C z x dz v dzx x B z x C z x x x x q q q q qq q J wt qq q qq q ÎG - + - = Y + + + - + ò - (35)
6. Conclusions
The Gaussian beam method has an important role in high-frequency wave field calculations and seismic migration. However, the difficulties during the calculation process faced by the Gaussian beam method in high-frequency wave fields and seismic imaging have not yet been solved – that is, the time-consuming coordinate transformation and also the phase distribution. The Eulerian Gaussian beam method has been proposed to solve these two problems. We derived the equations for the anisotropic Eulerian Gaussian beam. The directional derivative of the ray parameters was obtained from the Hamiltonian form of the eikonal equation and the x , z coordinates of the derivative of the travel time. We derived anisotropy travel time formulas based on the physical meaning of the direction factor ,as well as the relationship between the angle and .In addition, we obtained the ray tracing and dynamic ray tracing equations. The Eulerian Gaussian beam method is within the scope of ray theory and therefore we used the stationary phase method to obtain the weighting function and the changes in the weighting function with the source point of the travel time. In future studies, we plan to develop numerical algorithms for the high-frequency wave field and seismic migration using the anisotropic Eulerian Gaussian beam. References
Alkhalifah T (1995) Gaussian beam depth migration for anisotropic media. Geophysics 60: 1474 - { } (x, , ) : X( , , ) s z z x x q qG = u q u Y media: Gaussian beam approach. Geophysics Journal of the Royal Astronomical Society 88: 43–79. George T, Virieux J, Madariaga R (1987) Seismic wave synthesis by Gaussian beam summation: a comparison with finite differences. Geophysics 52: 1065–1073. Gray S H (2005) Gaussian beam migration of common shot records. Geophysics 70: S71–S77. Gray S H, Bleistein N (2009) True-amplitude Gaussian-beam migration. Geophysics 74: S11–S23. Hill N (1990) Gaussian beam migration. Geophysics 55: 1416–1428. Hill N R (2001) Prestack Gaussian beam depth migration. Geophysics 66: 1240–1250. Klimes L (1984) Expansion of high-frequency time harmonic wavefield given on an initial surface into Gaussian beams. Geophysics Journal of the Royal Astronomical Society 79: 105–118. Klimes L (1994) Transformations for dynamic ray tracing in anisotropic media. Wave Motion 20: 261–272. Kravtsov Yu A, Berczynski P (2007) Gaussian beams in inhomogeneous media: a review. Studies in Geophysics and Geodesy 51: 1–36. Leung S, Qian J, Osher S (2004) A level set method for three dimensional paraxial geometrical optics with multiple sources. Communications in Mathematical Sciences 2: 657–686. Leung S J, Qian J, Burridge R (2007) Eulerian Gaussian beams for high-frequency wave propagation. Geophysics 72: SM61–SM76. Nowack R L (2003) Calculation of synthetic seismograms with Gaussian beams. Pure and Applied Geophysics 160: 487 - - -
97. Popov M M, Semtchenok N M, Popov P M et al. (2010), Depth migration by the Gaussian beam summation method. Geophysics 75: S81–S93.11