Multiscale reverse-time-migration-type imaging using the dyadic parabolic decomposition of phase space

# Multiscale reverse-time-migration-type imaging using the dyadic parabolic decomposition of phase space

Fredrik Andersson Mathematics LTH, Centre for Mathematical Sciences, Lund Institute of Technology, Lund University, Sweden (fa@maths.lth.se)    Maarten V. de Hoop Department of Mathematics, Purdue University, West Lafayette, IN 47907, USA (mdehoop@purdue.edu).    Herwig Wendt CNRS, Institut de Recherche en Informatique de Toulouse (IRIT), Toulouse, France (herwig.wendt@irit.fr)
###### Abstract

We develop a representation of reverse-time migration in terms of Fourier integral operators the canonical relations of which are graphs. Through the dyadic parabolic decomposition of phase space, we obtain the solution of the wave equation with a boundary source and homogeneous initial conditions using wave packets. On this basis, we develop a numerical procedure for the reverse-time continuation from the boundary of scattering data and for RTM migration. The algorithms are derived from those we recently developed for the discrete approximate evaluation of the action of Fourier integral operators and inherit from their conceptual and numerical properties.

## 1 Introduction

Reflection seismology is a commonly used method to study the properties of Earth’s subsurface in geophysical exploration. Point sources are placed on Earth’s surface which generate acoustic waves in the subsurface that are reflected where the medium properties vary discontinuously. These reflections are recorded at Earth’s surface by (arrays of) point receivers. The goal of seismic imaging is to reconstruct the singular variations in medium properties from the reflected waves recorded at the surface [10, 4, 29]. The most common formulation for seismic inverse scattering takes the form of a linearized inverse problem for the medium coefficient in the acoustic wave equation, where the linearization is performed about a smoothly varying background. Here, the background model is assumed to be known. However, via a formulation as a separable inverse problem one can also proceed with determining this background model. The linearization defines a single scattering operator that maps the coefficient contrast to the data, i.e., the restriction of the scattered wave field to the acquisition set. The adjoint defines the process of imaging.

We consider reverse-time migration (RTM) [23, 32, 19, 3, 28], and the RTM-based inverse scattering transform developed and analyzed by Op’t Root et al. [20]. Through an appropriate formulation based on wave field continuation [16], we arrive at a representation of RTM in terms of a Fourier Integral Operator (FIO) associated with a canonical graph. Indeed, we use such a representation. The key contribution of this paper lies in the development of an algorithm for the solution of the wave equation with a boundary source and homogeneous initial conditions using the dyadic parabolic decomposition of phase space in the limit of fine scales. This algorithm is then composed with an imaging condition to yield the RTM-type imaging or inverse scattering. We explicitly admit the formation of caustics.

Viewing wave packets as localized plane waves, our approach has connections to methods in which one designs sources that favor (directional) illumination of particular subdomains of the subsurface. We mention plane-wave migration and beam-wave migration. In plane-wave migration one synthesizes plane-wave source experiments [33]. Given a plane-wave source one can then introduce tilted coordinates to carry out the wave field extrapolation with a limited accuracy propagator [24]. In beam-wave migration, Brandsberg-Dahl and Etgen [5] use a rotating coordinate system and essentially couple wave field methods with band-limited properties to ray-geometric methods. Furthermore, we mention the use of coherent states in this context by Albertin et al. [1]. Instead of tilted coordinates, one can use curvilinear coordinates in combination with a paraxial propagator [22]; the curvilinear coordinates may be generated as geodesics initiated from a point source or a plane wave. We note, however, that these methods are downward-continuation based whereas our approach is based on reverse-time continuation; also, we decompose the data although we could incorporate a synthesis of wave packet sources as well.

Our numerical solution is derived from the algorithm that we developed for FIOs [2]. Computational efficiency arises from organizing the decomposition and propagation by directions associated with frequency boxes instead of individual wave packets. The superposition of wave packets is complete and their propagation, as well as the corresponding imaging, converges in the limit of fine scales in smooth velocity models. Our formulation is insensitive to specific choices of (absorbing) boundary conditions, which is in contrast to PDE based solutions, including beam-wave migration. Moreover, it naturally conveys angular information which can be used in the imaging process, for instance, for computing restricted angle transforms. Candès, Demanet & Ying developed a fast butterfly algorithm for the application of Fourier integral operators associated with canonical graphs [9], which presents an interesting alternative to the propagation component of our algorithm.

Our algorithm is particularly well suited for application to (limited aperture) seismic array data providing a way of partial (in phase space) imaging possibly with a small set of sources. Moreover, if one need not generate a ‘global’ image, we do not need to evaluate the relevant wave field solutions at ‘all’ times, unlike algorithms based on numerically solving the wave equation, which enables computationally efficient target-oriented imaging. Target-oriented imaging can be used effectively, for example, with available arrays and earthquakes in studying heterogeneities and discontinuities in Earth’s mantle [18, 31]. A key element of our algorithm is finding a low-rank separated representation of the amplitude of the relevant FIO, which we do using Prolate Spheroidal Wave Functions (PSWFs) [2]. Demanet & Ying [12] proposed a method of finding such a representation based on the randomized sampling algorithm for constructing factorizations for low-rank matrices.

Our algorithm involves the propagation of high-frequency waves. To compare its complexity with the computational complexities of RTM algorithms based on numerically solving the wave equation, one can essentially compare the complexities of the backpropagation of a boundary source (from single-source data). Considering (backward) solving the wave equation in dimension on a grid of side length and a propagation time of order , the CFL condition implies that the time step is of order . The applications of a short-time propagator, with a presumed complexity , then yields a complexity . The time step in our algorithm is in principle, while the application of the Fourier integral operator representing the propagator is per frequency box (see [2]); the number of time steps to be computed is affected by the size of the target. Demanet & Ying [12] already pointed out the time upscaling of an approach to wave propagation using a FIO.

The outline of the paper is as follows. In Section LABEL:sec:2, we summarize the parametrix construction of the wave equation, and introduce the relevant Hamilton system and linearized Hamilton-Jacobi equations describing the geometry of the imaging process. In Section LABEL:sec:3 we formulate reverse-time continuation from the boundary and obtain a particular oscillatory integral representation for the kernel of this process, to which the algorithm for FIOs that we developed in an earlier paper applies. Subsection LABEL:sec:33 contains this key new result, and Subsection LABEL:sec:comprtc its computational counterpart. These are also the main components of the asymptotic form of the RTM-based inverse scattering transform and imaging algorithm, which we develop in Section LABEL:sec:4. In Section LABEL:sec:5, we give numerical examples of reverse-time continuation and inverse scattering also in the presence of caustics. We end with a discussion in Section LABEL:sec:6.

### Dyadic parabolic decomposition of phase space

We briefly discuss the (co)frame of curvelets and wave packets [8, 14, 25]. We will implicitly suppose that the data are decomposed into wave packets below, and we will develop wave packet based algorithms with accuracy [2].
Let represent a (seismic) velocity field, and let be the Fourier transform. One begins with covering the positive -axis () by overlapping boxes of the form

 Bk=[ξ′k−L′k2,ξ′k+L′k2]×[−L′′k2,L′′k2]n−1. \hb@xt@.01(1.1)

Here, both the centers and the side lengths , follow parabolic scaling

 ξ′k∼2k,L′k∼2k,L′′k∼2k/2,as k→∞.

Next, for each , let vary over a set of approximately uniformly distributed unit vectors111By convention, we let be aligned with the -axis.. Let denote a choice of rotation matrix which maps to , and

 Bν,k=Θ−1ν,kBk. \hb@xt@.01(1.2)

In the (co)frame construction, one encounters two sequences of smooth functions, and , on , each supported in , that form a copartition of unity

 ^χ0(ξ)^β0(ξ)+∑k≥1∑ν^χν,k(ξ)^βν,k(ξ)=1, \hb@xt@.01(1.3)

and satisfy the estimates

 |⟨ν,∂ξ⟩j∂αξ^χν,k(ξ)|+|⟨ν,∂ξ⟩j∂αξ^βν,k(ξ)|≤Cj,α2−k(j+|α|/2).

One now forms

 ^ψν,k(ξ)=ρ−1/2k^βν,k(ξ),^φν,k(ξ)=ρ−1/2k^χν,k(ξ), \hb@xt@.01(1.4)

where is the volume of . These functions satisfy the estimates

 ∀N:|φν,k(x)||ψν,k(x)|}≤CN2k(n+1)/4(2k|⟨ν,x⟩|+2k/2∥x∥)−N. \hb@xt@.01(1.5)

To obtain a (co)frame, one introduces the integer lattice: , the dilation matrix

 Dk=12π(L′k01×n−10n−1×1L′′kIn−1),detDk=(2π)−nρk,

and the points . The frame elements are now defined in the Fourier domain as

 ^φγ(ξ)=^φν,k(ξ)exp[−i⟨xν,kj,ξ⟩],γ=(j,ν,k),k≥1, \hb@xt@.01(1.6)

and similarly for . Thus, one obtains the transform pair

 uγ=∫u(x)¯¯¯¯¯¯¯¯¯¯¯¯¯ψγ(x)dx,u(x)=∑γuγφγ(x) \hb@xt@.01(1.7)

with the property that for each .

## 2 Parametrix

Here, we summarize the parametrix construction for the wave equation. We consider the Cauchy initial value problem

 [∂2∂t2+A(x,Dx)]u = 0,A(x,Dx)=c(x)D2xc(x), \hb@xt@.01(2.1) u(x,0) = 0,∂u∂t(x,0)=h(x); \hb@xt@.01(2.2)

we have normalized the pressure: .

To evaluate the parametrix, we use the first-order system for that is equivalent to this wave equation,

 ∂∂t(u∂u∂t)=(01−A(x,Dx)0)(u∂u∂t). \hb@xt@.01(2.3)

This system can be decoupled, namely, by the matrix-valued pseudodifferential operators

 V(x,Dx)=(11−iB(x,Dx)iB(x,Dx)),Λ(x,Dx)=12(1iB(x,Dx)−11−iB(x,Dx)−1),

where is a pseudodifferential operator of order 1.

The use of a general symbol in our presentation facilitates the extension of our algorithm to the imaging with elastic waves [6].

The principal symbol of is given by . Then

 u±=12u±12iB(x,Dx)−1∂u∂t, \hb@xt@.01(2.4)

satisfy the two first-order (“half wave”) equations

 P±(x,Dx,Dt)u±=0, \hb@xt@.01(2.5)

where

 P±(x,Dx,Dt)=∂∂t±iB(x,Dx), \hb@xt@.01(2.6)

supplemented with the initial conditions

 u±|t=0=h±,h±=±12iB(x,Dx)−1h. \hb@xt@.01(2.7)

We construct operators that solve the initial value problem (LABEL:eq:decoupled_firstorder), (LABEL:eq:hpolIVs): . Then . The operators are Fourier integral operators. Their construction is well known, see for example Duistermaat [17, Chapter 5]. Microlocally, the solution operator associated with (LABEL:eq:first_order_time) can be written in the matrix form

 S(t)=V(S+(t)00S−(t))Λ;

in this notation, .

For the later analysis, we introduce the operators and : solves the problem

 [∂∂t−(01−A(x,Dx)0)]S(t,s) = 0, S(⋅,s)|t=s = 0,∂S∂t(⋅,s)|t=s=Id,

so that the solution of

 [∂2∂t2+A(x,Dx)]u=f,u(t<0)=0,

is given by

 u(y,t)=∫t0P1S(t,s)(0f(⋅,s))(y)ds=∬G(y,x,t−s)f(x,s)dxds,

where we identified the causal Green’s function . Here, is the projection, . Likewise, solves (for ) the problem

 P+(x,Dx,Dt)S+(⋅,s) = 0, S+(⋅,s)|t=s = Id,

so that the causal solution of

 P+(x,Dx,Dt)u+=f+,f+=12iB(x,Dx)−1f,

is given by

 u+(y,t)=∫t−∞(S+(t,s)f+(⋅,s))(y)ds=∬G+(y,x,t−s)f+(x,s)dxds,

while the anticausal solution is given by

 u+(y,t)=−∫∞t(S+(t,s)f+(⋅,s))(y)ds=∬G+(y,x,s−t)f+(x,s)dxds.

A similar construction holds with replaced by .

### 2.1 Oscillatory integral representation

For sufficiently small (in the absence of conjugate points), one obtains the oscillatory integral representation,

 (S±(t)h±)(y)=(2π)−n∬a±(y,t,ξ)exp[iϕ±(y,t,x,ξ)]h±(x)dxdξ, \hb@xt@.01(2.8)

where

 ϕ±(y,t,x,ξ)=α±(y,t,ξ)−⟨ξ,x⟩. \hb@xt@.01(2.9)

We note that . Singularities are propagated along the bicharacteristics, which are determined by Hamilton’s equations generated by the principal symbol

 dytdt=±∂Bprin(yt,ηt)∂η,dηtdt=∓∂Bprin(yt,ηt)∂y. \hb@xt@.01(2.10)

We denote the solution of (LABEL:eq:Hamilton) with the sign and initial values at by . The solution with the sign is found upon reversing the time direction and is given by . Away from conjugate points, and determine and ; we write and . (We also use the parametrization in which the roles of and are interchanged.) Then

 α+(y,t,ξ)=⟨ξ,˜xt(y,ξ)⟩.

To highest order,

 a+(y,t,ξ)=∣∣\rule{0.0pt}{14.226378pt}∂(yt)∂(x)∣∣∣x=˜xt(y,ξ),ξ\rule{0.0pt}{14.226378pt}∣∣−1/2. \hb@xt@.01(2.11)

We consider the perturbations of with respect to the initial conditions ,

 Wt(x,ξ)=(Wt1(x,ξ)Wt2(x,ξ)Wt3(x,ξ)Wt4(x,ξ))=(∂xyt(x,ξ)∂ξyt(x,ξ)∂xηt(x,ξ)∂ξηt(x,ξ)). \hb@xt@.01(2.12)

This matrix solves the (linearized) Hamilton-Jacobi equations,

 dWtdt(x,ξ)=(∂ηyBprin(yt,ηt)∂ηηBprin(yt,ηt)−∂yyBprin(yt,ηt)−∂yηBprin(yt,ηt))Wt(x,ξ), \hb@xt@.01(2.13)

subject to initial conditions . We note that away from conjugate points, the submatrix is invertible. Because

 ˜xt=∂α+∂ξ,˜ηt=∂α+∂y,

integration of (LABEL:equ:perturbS) along yields:

 ∂2α+∂y∂ξ(yt(x,ξ),t,ξ) = (Wt1(x,ξ))−1, \hb@xt@.01(2.14) ∂2α+∂ξ2(yt(x,ξ),t,ξ) = (Wt1(x,ξ))−1Wt2(x,ξ), \hb@xt@.01(2.15) ∂2α+∂y2(yt(x,ξ),t,ξ) = Wt3(x,ξ)(Wt1(x,ξ))−1, \hb@xt@.01(2.16)

which we evaluate at . It follows that

 a+(y,t,ξ)=∣∣\rule{0.0pt}{8.535827pt}detWt1|x=˜xt(y,ξ),ξ\rule{0.0pt}{8.535827pt}∣∣−1/2.

The amplitude of , then becomes

 a+(y,t,ξ)12iBprin(˜xt(y,ξ),ξ)−1

to leading order; we denote this amplitude by . The amplitude follows from time reversal: .

In the case of conjugate points, we use the semigroup property of and decompose the time step into smaller time steps such that in each step the formation of caustics is avoided. Numerically, the size of the smaller time steps can be determined by monitoring the rank-deficiency of , see [11] for a more general point of view and Subsection LABEL:sec:comprtc for an application.

### 2.2 The source field

In the absence of caustics, we can change phase variables in the oscillatory integral representation of according to

 G+(y,x,t)=(2π)−1∫∫(2π)−n∫a+(y,t′,ξ)exp[iϕ+(y,t′,x,ξ)]dξexp[iτ(t−t′)]dt′dτ=(2π)−1∫a′+(y,x,τ)exp[iτ(t−T(y,x))]dτ. \hb@xt@.01(2.17)

By applying the method of stationary phase in the variables , one can show that the source field can be written in the form [6]

 G(x,~x,t)=(2π)−1∫a′(x,~x,τ)exp[iτ(t−T(x,~x))]dτ. \hb@xt@.01(2.18)

Here, is the source location and is the travel time satisfying the eikonal equation

 Bprin(x,−∂xT(x,~x))=−1 \hb@xt@.01(2.19)

and to highest order with

 |A(x,~x,τ)|=(2π)−(n−1)/2∣∣\rule{0.0pt}{12.8% 0374pt}det∂(x,ξ,t)∂(y,x,τ)% \rule{0.0pt}{12.80374pt}∣∣1/2 \hb@xt@.01(2.20)

see [6] for details. We introduce

 n~x(x)=∂xT(x,~x)|∂xT(x,~x)|; \hb@xt@.01(2.21)

in view of (LABEL:eq:eik),

 |∂xT(x,~x)|=1Bprin(x,n~x(x)).

We note that through we obtain the incidence angle of the source field at . In Section LABEL:sec:angle, we will arrange and study the images with respect to incidence angle. We also note that can be estimated from the Poynting vector at [35, 36] or from (possibly normalized by the autocorrelation, ; note that this normalization is primarily applied to suppress the dependency on ), for instance in the PDE solution formulation of RTM.

## 3 Reverse-time continuation from the boundary

The key results we obtain in this section are the formulation of an oscillatory integral representation and its computation using dyadic parabolic decomposition and wave packets for reverse-time continuation with a boundary source. These are also central in the formulation and computation of the inverse scattering and imaging operators presented in Section LABEL:sec:4. We introduce Euclidean boundary normal coordinates, ; that is, , and defines the boundary. We let denote a bounded open subset of . We denote the restriction to the boundary by .

We let be an anticausal solution to

 [∂2∂t2+A(x,Dx)]wr(x,t)=δ(xn)g(x′,t); \hb@xt@.01(3.1)

we have with

 wr,+(y,t)=−∫∞t(S+(t,s)12iB−1R∗xn˜ΨΣg(⋅,s))(y)ds

noting that

 R∗xng(x,t)=δ(xn)g(x′,t)

for any functions defined on . Here, is a pseudodifferential cutoff designed to remove grazing rays. The relation between contributions from negative frequencies and positive frequencies is

 wr,−(y,t)=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯wr,+(y,t). \hb@xt@.01(3.2)

We now introduce principal parts of symbols, , as the solutions for of

 Aprin(x′,xn,ξ′,ζ)=τ2.

We write . In the further analysis we will need the operator,

 C(x′,xn,D−1tDx′,1)at the surface, xn=0

with principal symbol .

### 3.1 Conjugate points

In the case of conjugate points, we introduce a partition of unity into (with overlap in time). Incorporating this partition of unity in , we obtain a set of cutoffs, . The first index signifies a subdivision in while the second index identifies intervals in time.

To describe the use of the semigroup property, we fix . Let and assume, without loss of generality, that we need two smaller time intervals, and , say, to avoid conjugate points in the smaller time intervals. Then

 ∫T1t(S+(t,s)12iB−1R∗xn(˜ΨΣ,11+˜ΨΣ,12)g(⋅,s))(y)ds=∫t+t1t(S+(t,s)12iB−1R∗xn˜ΨΣ,11g(⋅,s))(y)ds+S+(t,t+t1−δ)∫T1t+t1−δ(S+(t+t1−δ,s)12iB−1R∗xn˜ΨΣ,12g(⋅,s))(y)ds. \hb@xt@.01(3.3)

We now focus on representations for

 ∫t+t1t(S+(t,s)12iB−1R∗xn˜ΨΣ,11g(⋅,s))(y)dsand∫T1t+t1−δ(S+(t+t1−δ,s)12iB−1R∗xn˜ΨΣ,12g(⋅,s))(y)ds,

in the absence of conjugate points.

### 3.2 Oscillatory integral representations

We have

 χnwr,+(y,t)=(2π)−n∫∫∞t∫χna(bkd)(x′,s−t,y,η)exp[i(−α+(x′,0,s−t,η)+⟨η,y⟩)]g(x′,s)dx′dsdη, \hb@xt@.01(3.4)

where

 a(bkd)(x′,s−t,y,η)=∣∣\rule{0.0pt}{14.226378pt% }∂(ys−t)∂(x)∣∣∣η,x=xs−t(x′,0,η)\rule{0.0pt}{14.226378pt}∣∣−1/212iτ−1˜ΨΣ(x′,s,ξ′,τ) \hb@xt@.01(3.5)

up to terms of lower order, that is, the error (expressed in ) is of order , and is a cutoff function which removes contributions for (the expressions for and in terms of are given in (LABEL:eq:dalp2) and (LABEL:eq:dalp3) below). The operator is a FIO, the canonical relation of which is a subset of

 {(y,η;(ys−t)′(y,η),s−t,(ηs−t)′(y,η),−Bprin(y,η)) | (ys−t)n(y,η)=0}.

The dyadic parabolic decomposition of phase space enters in the reverse-time continuation as

 χn(yn)wr,+(y,t)=χn(yn)∑ν,k∬{\rule{0.0pt}{14.226378pt}(2π)−n∫ˆβν,k(η)ˆχν,k(η)a(bkd)(x′,s−t,y,η)exp[−iα+(x′,0,s−t,η)]dη\rule{0.0pt}{14.% 226378pt}}g(x′,s)dx′ds. \hb@xt@.01(3.6)

Fixing corresponds with (directional) “controlled illumination.”

### 3.3 Boundary source decomposition; wave packets in space-time

We change phase variables in the representation for . We could do this in two steps, changing parametrizations from to and then to . Here, we carry out this change in a single step:

 χnwr,+(y,t)=(2π)−2n∬∫∞t∬a(bkd)(x′,s−t,y,η)exp[i(−α+(x′,0,s−t,η)+⟨η,y⟩)]exp[i(τs+⟨ξ′,x′⟩)]dηdx′ds ˆg(ξ′,τ)dξ′dτ; \hb@xt@.01(3.7)

applying the method of stationary phase in yields solving

 ∂ηα+(x′,0,s−t,η) = y, \hb@xt@.01(3.8) ∂x′α+(x′,0,s−t,η) = ξ′, \hb@xt@.01(3.9) ∂sα+(x′,0,s−t,η) = τ \hb@xt@.01(3.10)

for given and fixed (which is viewed as a parameter here). The solutions, , are the stationary points of . We have . These equations imply that

 y=˜ys0−t(x′0,0,η0)ξ′=˜ξs0−t′(x′0,0,η0)} that is,  (x′0,0,˜ξs0−t′,C(x′0,0,˜ξs0−t′,τ))Φs0−t→(˜ys0−t,η0).

For given , is determined since , and determine a unique ray, in view of the absence of conjugate points. Thus we need to solve

 η0 = ˜ηs0−t(y,ξ′,C(x′0,0,ξ′,τ)), \hb@xt@.01(3.11) x′0 = ˜xs0−t′(y,ξ′,C(x′0,0,ξ′,τ)), \hb@xt@.01(3.12) 0 = ˜xs0−tn(y,ξ′,C(x′0,0,ξ′,τ))   or  s0=T(x′0,0,y)+t, \hb@xt@.01(3.13)

for . To obtain a unique solution, in general, we need to localize , which we do by substituting a wave packet contribution, that is, for . Then

 −α+(x′0(y,