A complex point-source solution of the acoustic eikonal equation

for Gaussian beams in transversely isotropic media

Running head: Complex point-source solutions

Xingguo Huang

1

, Hui Sun

2,3

, Zhangqing Sun

4

and Nuno Vieira da Silva

5

1

University of Bergen, Department of Earth Science, Allegaten 41, 5020 Bergen, Norway. Emails:

Xingguo.Huang@uib.no,xingguo.huang19@gmail.com.

2

Southwest Jiaotong University, Faculty of Geosciences and Environmental Engineering, Chengdu

610031, China; (Corresponding author) Email: sunhui@swjtu.edu.cn.

3

University of California, Earth and Planetary Sciences, Modeling and Imaging Laboratory, 1156 High

Street, Santa Cruz, CA 95064, USA.

4

Jilin University, College for Geo-exploration Science and Technology, 938 Ximinzhu Street,

Changchun 130026, China. Email: sun_zhangq@jlu.edu.cn.

5

Formerly Imperial College London, Department of Earth Science and Engineering, South Kensington

Campus, London SW7 2AZ, UK. Email: nuno.vdasilva@gmail.com.

ABSTRACT: The complex traveltime solutions of the complex eikonal equation are the basis of

inhomogeneous plane-wave seismic imaging methods, such as Gaussian beam migration and tomography.

We present the analytic approximations for complex traveltime in transversely isotropic media with a

titled symmetry axis (TTI), which is defined by a Taylor series expansion over the anisotropy parameters.

The formulation for complex traveltime is developed using perturbation theory and the complex point

source method. The real part of the complex traveltime describes the wavefront and the imaginary part

of the complex traveltime describes the decay of amplitude of waves away from the central ray. We

derive the linearized ordinary differential equations for the coefficients of the Taylor-series expansion

2

using perturbation theory. The analytical solutions for the complex traveltimes are determined by

applying the complex point source method to the background traveltime formula and subsequently

obtaining the coefficients from the linearized ordinary differential equations. We investigate the influence

of the anisotropy parameters and of the initial width of the ray tube on the accuracy of the computed

traveltimes. The analytical formulas, as outlined, are efficient methods for the computation of complex

traveltimes from the complex eikonal equation. In addition, those formulas are also effective methods

for benchmarking approximated solutions.

3

INTRODUCTION

The complex eikonal equation plays an important role in wave propagation problems of

inhomogeneous plane waves (Červený, 2001). The inhomogeneous plane wave is a solution of the

elastodynamics equation if the traveltime is complex-valued (Krebes and Le, 1994). Studies on the

complex eikonal equation fall into two groups. In the first group, the complex eikonal equation is used

to describe the traveltimes in attenuating media. The complex eikonal equation is used for ray-based

imaging and inversion (Huang et al, 2019) for an attenuating medium (i.e. the stiffness coefficients in

the frequency-domain are complex-valued). Early studies on the complex traveltimes are based on the

complex ray tracing method (Hearn and Krebes, 1990a, 1990b; Zhu and Chun, 1994; Thomson, 1997;

Chapman et al., 1999; Kravtsov et al., 1990; Egorchenkov and Kravtsov, 2001). The computation of the

complex ray tracing is very expensive because it is implemented in a high dimensional space. The

introduction of the real ray tracing method based on perturbation theory (Vavryčuk, 2008a, 2008b, 2010,

2012; Klimeš and Klimeš, 2011) into the complex traveltimes computation enables computing complex

traveltimes numerically.

The research work carried out in the 1980s is principally focused on developing a linear system of

equations for asymptotic solutions concentrated in the vicinity of a ray utilizing a complex eikonal

(Babich and Nomofilov, 1981; Nomofilov, 1982). That early reported work shows Hamiltonian

formulations for the complex solution in some special systems of local coordinates. Although the theory

is well established, numerical implementation for Hamiltonian formulations is not known. Effectively,

solutions for Hamiltonian formulations require solving a system of differential equations numerically.

In the second group, the complex eikonal equation is used to describe the complex traveltimes of

Gaussian beams. Since Deschamps’s (1971), Wang and Deschamps’ (1974) and Felsen’s (1976) original

formulations of the complex phase and the complex eikonal equation method, the complex traveltimes

method has long attracted interest. Magnanini and Talenti (1999, 2003, 2006) use the Bäcklund transform

to conduct some extensive research on the complex eikonal equation. The Bäcklund transform constitutes

the basis of soliton theory (Terng and Karen, 2000). More recently, Li et al. (2011) develop a fast

marching method to solve the complex eikonal equation without the limitation of the paraxial ray

approximation. Accordingly, Huang et al. (2016b) design a local algorithm to compute the local complex

traveltimes by combining the fast marching method with the nonuniform grid-based finite difference

method. The complex traveltimes of the evanescent wave, in which their real part describes the wavefront

while their imaginary part describes the decay away from the central ray, play an important role in local

plane wave. Because of its advantages for handling caustics and multipath resulting from ray tracing in

complex media, the complex traveltimes have been widely used for Gaussian beam modelling (Červený,

4

1982; Popov, 1982; Huang et al., 2016a; Huang, 2018) and seismic imaging (Hill, 1990, 2001; Gray,

2005, 2009).

The work outlined herein sits in the second group. In this paper, we focus on solving the Gaussian

beam problem in a non-attenuating medium (i.e., the medium behaves elastically, and the stiffness

coefficients are real-valued). We use a “real-valued” eikonal equation, but we force the traveltimes to be

complex-valued, while Hearn and Krebs (1990a, 1990b), Zhu and Chun (1994), Thomson (1997),

Chapman et al. (1999), Kravtsov et al. (1999), Egorchenkov and Kravtsov (2001) and Hao and

Alkhalifah (2017a, b) use complex-valued eikonal equation without any enforcement of the traveltimes.

Previously, the condition that the gradient of the real part of the traveltimes is orthogonal to that of the

imaginary part is taken to solve the complex eikonal equation. That condition does not hold true for

traveltimes determined from the complex eikonal equation for an attenuating isotropic medium. A key

advantage of basing our discussion on Gaussian beam theory is that our expressions are readily applicable

to complex heterogeneous media which is paramount for accurate seismic imaging (Hill, 2001).

When media are anisotropic, the isotropic approximation is no longer adequate (Huang and

Greenhalgh, 2019). Areas with thin layering and fractures are particularly common in the geophysical

exploration context, therefore the effect of anisotropy needs to be considered in realistic applications of

the complex eikonal equation. This paper outlines the computation of the traveltimes for TTI media.

Červený (2001) introduces the complex ray tracing method for general anisotropic elastic media. It is

important to note that the application of such formulations is not usually feasible as the vast majority of

the recorded datasets do not contain complete information on the elements of the full anisotropic stiffness

tensor. In addition, it is unfeasible to consider the full anisotropic characterization of elastic media . It is

often verified that disregarding shear waves from the constitutive law for an anisotropic medium yields

an accurate representation of the kinematics of waves at the expense of sacrificing the accuracy of the

computed amplitudes. However, it is important to note that modeling seismic waves with the true

amplitude is generally not feasible as the exact source wavelet is unknown. In addition, the dynamic

response of an acoustic medium is different from that of an elastic medium. The acoustic approximation

has been widely used in seismic inversion and imaging (Alkhalifah, 2000, Grechka, 2004, Zhang et al.,

2011, Duveneck and Bakker, 2011, da Silva et al, 2016). The work follows the same rationale regarding

the constitutive law. The real eikonal approximation method is originally proposed as a means of

computing the real-valued traveltimes for anisotropic media (Alkhalifah, 2011a, 2011b, 2013; Stovas

and Alkhalifah, 2012; Waheed et al., 2013; Hao and Alkhalifah, 2017a, 2017b; Hao et al., 2018). In these

methods, a Taylor expansion of the traveltimes with respect to the anisotropy parameters and to the

background traveltimes for tilted elliptically isotropic media is carried out. The major advantage of these

methods is that if we know the background traveltimes and anisotropy parameters, we can obtain

traveltimes for general anisotropic media directly at a relatively low computational cost. More recently,

5

Huang et al. (2018) extend the mentioned above perturbation theory to the real and imaginary parts of

the complex traveltimes, and they solve successfully the complex eikonal equation in VTI and

orthorhombic media (Huang and Greenhalgh, 2018). Unlike the real eikonal approximation method, we

take Taylor expansions of both the real and imaginary parts of the complex traveltimes. Then, we

construct the linearized complex eikonal equations by substituting the Taylor’s expansions into the

complex eikonal equations.

According to the recent results of Huang et al. (2016b), solving the complex eikonal equation leads

to an inverse process for computing the imaginary slowness. This process is especially complicated for

anisotropic media. This is because the complex eikonal equation for anisotropic media includes several

parameters: the normal moveout (NMO) velocity, velocity along the symmetry axis, and the anellipticity

parameter, .

This paper focus on two principal aspects. The first aspect is deriving the linearized complex eikonal

equations based on perturbation theory (Alkhalifah, 2011a, 2011b, 2013; Stovas and Alkhalifah, 2012;

Waheed et al., 2013; Hao et al., 2016; Hao and Alkhalifah, 2017a, 2017b; Hao et al., 2018). A key

advantage of the linearized complex eikonal equations is providing possible method of solving the TTI

complex eikonal equation using a finite-difference method, which can be used to address the computation

of traveltimes in general heterogeneous media. Once the TTI complex eikonal equation can be solved

directly, the paraxial ray approximation can be avoided while solving the complex eikonal equation.

Following the recent methods (Huang and Greenhalgh, 2018; Huang et al., 2018), we develop a linear

partial differential equation for the TTI complex eikonal equation. However, we use a tilted elliptically

isotropic background medium for the perturbation expansion associated with the tilt angle, which allows

us to perform the perturbation expansion in η only. This reduces uncertainty in the computed complex

traveltimes.

The second aspect is obtaining the analytical formulas for homogeneous TTI media. The analytical

solutions can be used to benchmark the accuracy of alternative methods, as well as, providing insight on

the kinematics of wave propagation in TTI media. As pointed out earlier the latter has found many

relevant applications in applied geophysics, especially in seismic imaging and inversion. Even though

the theory developed in this paper applies to the general complex traveltimes for TTI media in the

complex space, our closed form (analytical) solution is for a specific Gaussian beam. This is based on

the fact that, for a complex point source (Choudhary and Felsen, 1974; Wang and Deschamps, 1974;

Felsen, 1976, 1984; Wu, 1985), the beam has an initial constant phase front and a Gaussian field profile,

and the wavefield has a Gaussian decay along the direction perpendicular to the central ray. In this case,

the contours of the real part of the traveltimes define the equiphase contours, while the contours of the

imaginary part define the phase paths. In addition, the theoretical development of the analytical formulas

is based on another assumption that the distance is limited to the region of the paraxial ray approximation.

η