{'fraction_non_alphanumeric': 0.05562800165389142, 'fraction_numerical': 0.010400432556216406, 'mean_word_length': 4.2056291390728475, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multiscale reverse-time-migration-type imaging using the dyadic parabolic decomposition of phase space

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.


Introduction
In reflection seismology one places point sources and point receivers on the earth's surface. A source generates acoustic waves in the subsurface, which are reflected where the medium properties vary discontinuously. In seismic imaging, one aims to reconstruct the properties of the subsurface from the reflected waves that are observed at the surface [8,4,23]. There are various approaches to seismic imaging, each based on a different mathematical model for seismic reflection data with underlying assumptions. In general, seismic scattering and inverse scattering have been formulated in the form of a linearized inverse problem for the medium coefficient in the acoustic wave equation. The linearization is around a smoothly varying background, called the velocity model, which is a priori also unknown. However, in the inverse scattering setting considered here, we assume the background model to be known. The linearization defines a single scattering operator mapping the model contrast (with respect to the background) to the data, that consists of the restriction to the acquistion set of the scattered field. The adjoint of this map defines the process of imaging in general.
We consider reverse-time migration (RTM) [17,24,14,3,22], and the RTM-based inverse scattering transform developed and analyzed by Op't Root et al. [15]. Through an appropriate formulation of the wave field continuation approach, we arrive at a representation of RTM in terms of a Fourier Integral Operator (FIO) associated with a canonical graph. We revisit the original reverse-time imaging procedure. An RTM migration algorithm constists of three main parts: The modeling of source wave propagation in forward time, the modeling of receiver or reflected wave propagation in reverse time, and the application of the so called imaging condition [8,4].
The key contribution of the present work lies in the solution of the wave equation with a boundary source and homogeneous initial conditions using wave packets in the limit of fine scales. We employ the dyadic parabolic decomposition of phase space and develop an efficient procedure for the numerical evaluation of the solution operator. As an application, we consider reversetime imaging by incorporation of the imaging condition in the receiver wave field boundary source reverse time solution.
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 planewave migration one synthesizes plane-wave source experiments [25]. Given a plane-wave source one can then introduce tilted coordinates to carry out the wavefield extrapolation relatively accurately with a limited accuracy propagator [18]. In beam-wave migration, Brandsberg-Dahl and Etgen [5] use a rotating coordinate system and essentially couple wavefield 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 [16]; 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 approach is derived from the algorithm that we developed for Fourier integral operators [2]. 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.
The outline of the paper is as follows. In Section 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 3 we discuss reverse time continuation from the boundary and obtain a particular oscillatory integral representation for the kernel of this process, to which the algorithm for Fourier integral operators that we developed in an earlier paper applies. In Section 4 we give the asymptotic form of the RTMbased inverse scattering transform. In Section 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 6.

Dyadic parabolic decomposition of phase space
We briefly discuss the (co)frame of curvelets and wave packets [7,11,19]. We will implicitly suppose that the data are decomposed into wave packets below, and we will develop wave packet based algorithms with accuracy O(2 −k/2 ) [2].
Let u ∈ L 2 (R n ) represent a (seismic) velocity field, andû(ξ) = u(x) exp[−i x, ξ ] dx be the Fourier transform. One begins with covering the positive ξ 1 axis (ξ = ξ 1 ) by overlapping boxes of the form Here, both the centers ξ k and the side lengths L k , L k follow parabolic scaling Next, for each k ≥ 1, let ν vary over a set of approximately 2 k(n−1)/2 uniformly distributed unit vectors 1 . Let Θ ν,k denote a choice of rotation matrix which maps ν to e 1 , and In the (co-)frame construction, one encounters two sequences of smooth functions,χ ν,k andβ ν,k , on R n , each supported in B ν,k , that form a co-partition of unitŷ and satisfy the estimates One now formsψ where ρ k is the volume of B k . These functions satisfy the estimates ∀N : To obtain a (co)frame, one introduces the integer lattice: X j := (j 1 , . . . , j n ) ∈ Z n , the dilation matrix The frame elements are now defined in the Fourier domain asφ and similarly forψ γ (ξ). Thus, one obtains the transform pair with the property that γ : k =k, ν =ν u γ φ γ (ξ) =û(ξ)β ν,k (ξ)χ ν,k (ξ) for each ν, k.

Parametrix
We consider the Cauchy initial value problem we have normalized the pressure: u = c −1 p.
To evaluate the parametrix, we use the first-order system for u that is equivalent to this wave equation, ∂ ∂t This system can be decoupled, namely, by the matrix-valued pseudodifferential operators is a pseudodifferential operator of order 1. The principal symbol of B(x, D x ) is given by B prin (x, ξ) = A prin (x, ξ). Then satisfy the two first-order ("half wave") equations where supplemented with the initial conditions We construct operators S ± (t) that solve the initial value problem (12), (14): The operators S ± (t) are Fourier integral operators. Their construction is well known, see for example Duistermaat [13], Chapter 5. Microlocally, the solution operator associated with (10) can be written in the matrix form For the later analysis, we introduce the operators S(t, s) and S ± (t, s): S(t, s) solves the problem so that the solution of is given by where we identified the causal Green's function G(y, x, t − s). Here, P 1 is the projection, P 1 u ∂u ∂t = u. Likewise, S + (t, s) solves (for t ∈ R) the problem S + (·, s)| t=s = Id, so that the causal solution of while the anticausal solution is given by A similar construction holds with + replaced by −.

Oscillatory integral representation
For sufficiently small t (in the absence of conjugate points), one obtains the oscillatory integral representation, where We note that α − (y, t, ξ) = −α + (y, t, −ξ). Singularities are propagated along the bicharacteristics, that are determined by Hamilton's equations generated by the principal symbol ±B prin (x, ξ) We denote the solution of (17) with the + sign and initial values (x, ξ) at t = 0 by (y t (x, ξ), η t (x, ξ)) = Φ t (x, ξ). The solution with the − sign is found upon reversing the time direction and is given by (y −t (x, ξ), η −t (x, ξ)). Away from conjugate points, y t and ξ determine η t and x; we write x = x t (y, ξ) and η t = η t (y, ξ). (We also use the parametrization in which the roles of (y, ξ) and (x, η) are interchanged.) Then α + (y, t, ξ) = ξ, x t (y, ξ) .
To highest order, We consider the perturbations of (y t , η t ) with respect to the initial conditions (x, ξ), This matrix solves the (linearized) Hamilton-Jacobi equations, subject to initial conditions W t=0 = I. We note that away from conjugate points, the submatrix W t 1 is invertible. Because integration of (20) along (y t , η t ) yields: which we evaluate at x = x t (y, ξ). It follows that The amplitude of S + (t) 1 2 iB −1 , then becomes a + (y, t, ξ) 1 2 iB prin ( x t (y, ξ), ξ) −1 to leading order; we denote this amplitude by a + (y, t, ξ). The amplitude a − follows from time reversal: a − (y, t, ξ) = a + (y, t, −ξ).
In the case of conjugate points, we use the semi-group property of S(t, s) and decompose the time step into smaller time steps such that in each step the formation of caustics is avoided.

The source field
In the absence of caustics, we can change phase variables in the oscillatory integral representation of G according to (24) We find the leading-order contribution, A + say, to a + by applying the method of stationary phase in the variables (ξ, t ): , we find that Furthermore, φ + (y, T (y, x), x, ξ(y, x, τ )) = 0. Thus the source field can be written in the form Here,x is the source location and T is the travel time satisfying the eikonal equation, to highest order, a = A with We also write Ax(x, τ ) = A(x,x, τ ). We introduce in view of (29), .
We note that through nx(x) we obtain the incidence angle of the source field at x. We also note that nx(x) can be estimated from the Poynting vector (at t = T (x,x)), [27,28] or from (possibly normalized by the autocorrelation, , for instance in the PDE solution formulation of RTM. In Section 5.3, we will arrange and study the images with respect to incidence angle.

Reverse time continuation from the boundary
We introduce Euclidean boundary normal coordinates, x = (x , x n ); that is, x = (x 1 , . . . , x n−1 ) and x n = 0 defines the boundary. We let Σ denote a bounded open subset of {(x, t) ∈ R n x × R t | x n = 0}. We denote the restriction to the boundary by R xn .
We let w r be an anti-causal solution to we have w r = w r,+ + w r,− with for any functions g defined on R n−1 is a pseudodifferential cutoff designed to remove grazing rays. The relation between contributions from negative frequencies and postive frequencies is For the later analysis, we introduce the square-root operator,

Conjugate points
In the case of conjugate points, we introduce a partition of unity in Σ ⊂ R n−1 x ×R t (with overlap δ in time). Incorporating this partition of unity in Ψ Σ , we obtain a set of cutoffs, Ψ Σ,ij . The first index signifies a subdivision in R n−1 x while the second index identifies intervals in time.
To describe the use of the semi-group property, we fix i. Let i = 1 and assume that we need two time intervals, [t, t + t 1 ] and [t + t 1 , T 1 ], say. Then We now focus on representations for in the absence of conjugate points.

Oscillatory integral representations
Microlocally, mod S −2 , and χ n is a cutoff function which removes contributions for x n < 0. The operator The dyadic parabolic decomposition of phase space enters in the reverse-time continuation as Fixing (ν, k) corresponds with (directional) "controlled illumination".

Boundary source decomposition, wave packets in space-time
We change phase variables in the representation for w r,+ . We could do this in two steps, changing parametrizations from ((x , x n ), η) to (y, (ξ , ξ n )) and then (s, ξ n ) to (x n , τ ). Here, we carry out this change in a single step: applying the method of stationary phase in (η, x , s) yields solving for given (y, ξ , τ ) and t fixed (which is viewed as a parameter here). The solutions, (η 0 , x 0 , s 0 ), are the stationary points of −α + (x , 0, s − t, η) + η, y + τ s + ξ , x ]. We have s 0 > t. These equations imply that for given t, s 0 is determined ((x 0 , 0), η 0 and s 0 − t determine a unique ray, in view of the absence of conjugate points). Thus we need to solve for (η 0 , x 0 , s 0 ). To obtain a unique solution, in general, we need to localize g, which we do by substituting a wave packet contribution, that is, g γ ϕ γ (ξ , τ ) for g(ξ , τ ). Then We evaluate and subject to the substitutions according to (39)-(41), and then and From these expressions we form the Hessian which is used to transform a (bkd) (x , s − t, y, η) to the amplitude a (y, t, ξ , τ ), so that Essentially, this representation corresponds with local coordinates (x n , y, ξ , τ ) for the canonical relation of the solution operator with t fixed.

Computation
We adapt the "box algorithm" for the multiscale discrete approximation of Fourier integral operators developed in [2] to (50), with accuracy O(2 −k/2 ) at frequency scale k. In the general case, the medium can be laterally varying at the boundary. Then we need to employ compactly supported cutoff functions in x , realized by the partition of unity Ψ Σ,ij . Within each cutoff i the lateral variation of the (smooth) velocity model is negligible at the boundary, and the algorithm outlined below can then be applied for each cutoff i separately. We now assume that the wave speed does not vary laterally at the surface and fix i = 1. Without loss of generality, we suppose that we need N s time intervals [t + (n s − 1)t 1 , t + n s t 1 ], n s = 1, . . . , N s , of duration t 1 = (T 1 − t)/N s in order to avoid the formation of conjugate points.
The solution of the corresponding Hamilton-Jacobi system yields the propagator matrix W s from which we obtain the quantities We can now apply the box-algorithm to each time interval n s and obtain w (ns) where α To obtain the final reverse time continued wave fields w (ns) r,+ (y, t), we construct a parametrix for the Cauchy initial value problems for the half-wave equation with initial data w (ns) r,+ (y, t + (n s − 1)t 1 ), n s = 2, . . . , N s , initial time t + (n s − 1)t 1 and final time t. We compute these parametrices using the box algorithm (this has been studied in detail in [2]). We make use of the semi-group property and obtain the parametrix for the reverse time interval [t + (n s − 1)t 1 , t] as the composition of the parametrices for the time intervals [t + (n s − n p + 1)t 1 , t + (n s − n p )t 1 ], n p = 2, . . . , n s . Finally, we have The different steps involved in modeling receiver wave propagation from the boundary in reverse time are summarized in Algo. 1 and illustrated in Fig. 2 for a numerical example detailed in Section 5.
We assume that We let w r be an anti-causal solution of (32) with let operators L and K be defined by Kw(y) = w(y, Tx(y)).
Operator K is a restriction to a hypersurface in R n+1 . Operator ∂ − n+1 2 t is to be read as the pseudodifferential operator with symbol τ →σ(τ )(iτ ) − n+1 2 in whichσ is a smooth function, valued 1 except for the origin where it is 0. The imaging operator, H, is then defined as Hd Σ (y) = (KL(w r,+ + w r,− ))(y).
To leading order symbols, we get i [∂ t θ + (y, t, ξ , τ ) + c(y)nx(y) · ∂ y θ + (y, t, ξ , τ )] exp[i θ + (y, t, ξ , τ )] g(ξ , τ ) dξ dτ. (55) We can use (33) in the computations. Through a simple modification, we can incorporate the imaging condition in the box-algorithm for reverse time continuation from the boundary detailed in Section 3.4, yielding a RTM imaging algorithm. Without loss of generality, we assume here that the source signature is a delta function; general discrete source signatures can be accommodated for in a straight-forward way by viewing them as a weighted sum of delta functions shifted by the time step of the Hamilton-Jacobi equations. Suppose that the source travel time Tx(y) and amplitude Ax(y) have been evaluated for the image region (here, by evaluating the Hamiltonian and Hamilton-Jacobi system of the half-wave equation, i.e. raytracing, cf. Section 2.2). We begin with the evaluation of H for partial reverse time continuation from the boundary (cp. Algo. 1, Part I). We obtain a contribution of time interval n s to the image at y * = y * (x , 0, s * ) if (n s − 1)t 1 + t ≤ Tx(y * ) ≤ n s t 1 + t − s * . Subject to this condition, the coordinate transform for image point y * is given by T * (ns) ν,k (y * ) = (x 0 , Tx(y * ) + s * ) and the propagator matrices are given by W s * . Application of the box algorithm with T * (ns) ν,k and W s * with Ax incorporated in the amplitudeã * yields the partial image ∆ * (ns,t+(ns−1)t1) (y) = ν,kã * (y, ν) In the case of conjugate points (n s > 1), we proceed with the evaluation of H for the subsequent half-wave equation reverse time continuation of the wave fields w (ns) r,+ (y, t + (n s − 1)t 1 ) (cp. Algo. 1, Part II). Consider continuation of w (ns) r,+ (y, t + (n s − 1)t 1 ) to w (ns) r,+ (y, t + (n s − 2)t 1 ) (n p = 2). In this process, we compute the bicharacteristics (ỹ s (y 0 , η 0 ),η s (y 0 , η 0 )) for s ∈ (0, t 1 ]. We can easily monitor the condition s * = t + (n s − 1)t 1 − Tx(ỹ s * ) under which we obtain an image contribution atỹ s * . The associated coordinate transform is given byT * (ns) ν,k (ỹ s * ) = y 0 , and the propagator matrices byW * = W (y 0 , η 0 , s * ). With these ingredients, application of the box algorithm yields the partial image∆ * (ns,np) (y); we obtain the final image ∆ dΣ (y) = Ns ns=1 ∆ * (ns,t+(ns−1)t1) (y) + Ns np=2 Ns ns=np∆ * (ns,np) (y).
The structure of the inverse scattering procedure is summarized in Algo. 2. Note that in the evaluation of the partial images ∆ * and∆ * , we can gather the incident angles η(y * ; ν, k; n s , n p ) of the reverse time continued wave field, which we can, for instance, use for monitoring scattering angles, see Section 5.3.

Numerical examples
We illustrate the theoretical and algorithmic developments in this work with two numerical examples. Although applicable in general dimension, we restrict ourselves here to n = 2 dimensions for clarity and convenience in illustration of the underlying computational structure. In a first example, we illustrate reverse time continuation from the boundary of a reflected wave field in the presence of conjugate points. Then, we use reverse time continuation for imaging of conormal singularities.
The data in these examples, g(x , t) and d Σ (x , t), respectively, are generated using time domain finite difference. The computational domain is fixed to N × N = 512 × 512 samples. r,+ (y, t + (n s − 1)t 1 ), n s = 1, . . . , 4 (top row); reverse time continued wave fields w (ns) r,+ (y, t + (n s − n p )t 1 ), n s = n p , . . . , 4 for n p = 2, 3 and 4, respectively (rows 2 to 4). Full reverse time continued wave field w r,+ (y, t = 0) (bottom left corner). 1 (center)). Despite several re-decomposition steps using the discrete wave packet transform involved in computing the reverse time continuation (semi-group), there are no visible amplitude fluctuation artifacts. In particular, we note that the edges of the cusp in the data are well focussed.
Finally, note that time intervals 1 and 2 do not contain any significant energy. While with the proposed procedure, it is possible to only compute the wave field for time intervals 3 and 4 (hence, no computation time and memory would be used for time slice 1 and 2), we chose to compute time intervals 1 and 2 and include them in Fig. 2 for completeness and clarity.

Imaging of conormal singularities
We proceed with a numerical illustration for imaging of conormal singularities by reverse time continuation from the boundary using the wave-packet based computational procedure developed in Section 4. The velocity model is plotted in Fig. 3 (top left). It consists of a decentered Gaussian low velocity (30% peak contrast with respect to the background velocity) and contains several horizontal line reflectors and one deep tilted line reflector. The (normal incident) reflectivity of the line reflectors varies with location and is documented in Fig. 5 (left). The data are generated using time domain finite difference and a Ricker wavelet with a peak frequency of 7Hz. The single source is located at the center of the boundary,x = 0. In Fig. 4, we plot the wave fields generated in the sub-surface for several time instances (for better visibility, we subtracted the wave field that is obtained when the line reflectors are not present). Despite the simple model, we observe relatively complex wave fields and, for late time instances, the formation of caustics. Furthermore, we observe artifacts from non-ideal absorbing boundaries, from multiple reflections, and in particular ringing artifacts for the wave fronts after the formation of caustics at large time instances. These are also present in the data, which we plot in Fig. 3 (center).
We need N s = 4 time intervals to avoid conjugate points within each semi-group step in the computational procedure described in Section 4 and outlined in Algo. 2. This partitioning in time intervals is indicated in Fig. 3 (right) on the data d Σ (x , t) after de-recomposition using the discrete wave packet transform. Again, as mentioned in the previous subsection, we can regularize and pre-process the data during this de-recomposition step.
We approximate the source signature with a single delta function at its temporal maximum and compute the source wave field by evaluating the Hamiltonian and Hamilton-Jacobi equations (dynamic) ray tracing). In Fig. 6, we plot the the partial images and reverse-time continued wave fields produced by the procedure described in Section 4, organized according to its hierarchical semi-group structure (cp. Fig. 2). Each column corresponds with one time interval of the data (from left to right, data slice n s = 1, . . . , 4, respectively), and transition from (group of) row(s) Figure 4: The wave field generated by a source atx = 0 for different time instances after subtraction of the reference wave field which is obtained when the line reflectors (indicated in white) are not present.
i to (group of) row(s) i + 1 corresponds with a semi-group re-decomposition and subsequent half-wave equation reverse time continuation and partial imaging step: The top row shows the snapshots w (ns) r,+ (y, t + (n s − 1)t 1 ), n s = 1, . . . , 4 produced by partial reverse time continuation of the 4 data slices (Algo. 2, Part I). The corresponding partial image ns ∆ * (ns,t+(ns−1)t1) (y) evaluated during this step is plotted in the bottom left corner of Fig. 6. At this stage, data slice 1 is fully reverse time continued (t = 0) while data slices 2 to 4 will be further reverse time continued after a semi-group re-decomposition (and enter Part II of Algo. 2). The second and third group of rows plot the output of Part II (cp. Algo. 2) of the procedure for n p = 2 and n p = 3, respectively: w (ns) r,+ (y, t+(n s −n p )t 1 ) (top rows) and∆ * (ns,np) (y) (bottom rows). We stop the semigroup iteration at n p = 3 because the energy of the data wave fields w 4 ns=np∆ * (ns,np) (y) we obtain, which we plot in Fig. 5 (right). We observe that all the reflectors are imaged correctly and well focussed, regardless of their depth, dip angle and background velocity. Note that we could further focus the image by using the full source signature (here, we use a delta source approximation). Certain reflectors are partially outside of the zone of illumination (e.g. the two rightmost reflectors at depth y 2 = 3.2 and y 2 = 6.5) and hence produce smiling "tails" due to the truncation of the wave field in the data (cp. Kirchhoff migration). Similarly, the corners of the line reflectors act as point diffractors and produce tails according to partial illumination and restricted geometry. We note that the ringing artifacts in the data components stemming from the two deepest reflectors are also present in the imagethe algorithm images what is in the data, and the data are imperfect. This is also the case for the artifact at depth y 2 = 1.3 in the image, which results from an imperfectly removed direct arrival (cp. Fig. 3 (right), (y 1 , t) = (3, 2.8)).
We note that the proposed wave packet based reverse time continuation procedure provides a natural and versatile framework for partial reconstruction by enabling phase-space localized control for the data (scale, orientation, position of the data wave packets) as well as directly for the image (scale, orientation and position of reverse time continued data wave packets; full angular information such as scattering angle, reflector dip angle).

Restricted angle transform
Having based our reverse time continuation and imaging procedure on dyadic parabolic decomposition and wave packets, we can extract angular information for the source and scattered wave fields. Indeed, for a given frequency boxχ ν,k , we have in each image point knowledge of the incident angles of the wave fronts. This information can be directly used in the imaging process, for instance for the restriction of the image to small scattering angles, or to certain target dip angles.
First, for a single source, we display the images obtained in the previous subsection as a function of incidence angle at the image points. Whereas geometrically the image of a singularity at one surface location is significant at one incidence angle only, the finite-frequency content of the wave packets result in a slight spread around the specular reflection angle. We indicate the incidence angle at the specular reflection with a red dot; see Fig. 7.
Secondly, we evaluate images of the singularities for multiple sources and re-arrange them in terms of local incidence angle [10,26,27,29]. In case the correct background velocity function is used, up to illumination effects, the images generated at different angles are the same; this reflects a redundancy in such data. If we perturb the background by moving the smooth lens, we still obtain a coherent images; however, the singularities move with changing incidence angle, see Fig.8. This behavior can be exploited to develop a procedure for reflection tomography [20,6]. r,+ (y, t + (n s − n p + 1)t 1 ) for n p = 2, n s = 2, . . . , 4 (center top rows) and for n p = 2, n s = 2, . . . , 4 (center top rows): snapshots w (ns) r,+ (y, t + (n s − n p )t 1 ) and partial images∆ * (ns,np) (y); partial images produced by Part I (bottom row, left) and by Part II for n s = 2, 3, 4 (bottom right).

Discussion
We have obtained a representation of RTM in terms of a FIO associated with a canonical graph and developed a solution to the wave equation with a boundary source and homogeneous initial conditions using dyadic parabolic decomposition of phase space and wave packets. We explicitly admit the formation of caustics. We devised a numerical procedure for its discrete evaluation which is derived from the algorithm for the evaluation of the action of Fourier integral operators through approximations [2], yielding accuracy O(2 −k/2 ) at frequency scale k. Our algorithm is organized by frequency boxesχ ν,k of the wave packet transform. Computations can be performed independently per individual box and are hence embarrassingly parallel. We obtain an effective one-step multi-scale procedure for reverse time continuation from the boundary for a given time interval, from T 1 to t, say. In this process, we can apply the imaging condition and obtain a reverse-time-migration imaging algorithm. Note that reverse-time continuation from the boundary yields also the foundations for an algorithm for the single-layer potential operator [9]. While numerical illustrations have been devised here for n = 2 dimensions, the concepts and computational procedures are valid for arbitrary dimension.
In the presence of conjugate points, we split the time interval for reverse time continuation into a sequence of smaller time intervals and reverse time continue partial wave fields subsequently for these time intervals using the semi-group property of the RTM operator. Numerically, this implies one discrete wave packet transform re-decomposition of the wave fields for each transition point from one time interval to another. After the first semi-group re-decomposition, reverse-time continuation essentially reduces to the evaluation of the wave equation for the propagation of an initial wave field, and any of the algorithms developed in [2] could be used as a computational basis. Here, we proposed a "box-algorithm" due to its favorable computational complexity and practical accuracy.
The computational complexity of our algorithm is of the order O(N n log(N )) per frequency box for each semi-group step (with a total number of O(N (n−1)/2 ) boxes; note however that depending on the data and the imaging target not all boxes will need to be computed). The computational complexity arises essentially from the complexity of the unequally spaced FFTs involved in the box algorithm (cf. [2] for details). Reverse-time-migration imaging increases the computational cost roughly by a factor two with respect to reverse time continuation of the wave field from the boundary because of the additional unequally spaced FFTs that need to be evaluated in order to create the image. Note that with the exception of the source wave field travel times and amplitudes and one snapshot during the semi-group re-decomposition, our procedure does not require the computation and storage of snapshots. Its computational cost and memory requirements remain of the order of the one-step evaluation of Cauchy initial value problems for evolution equations [2].
Evaluation of the RTM operator for all wave packets associated with a given frequency boxeŝ χ ν,k at once requires the existence of a homogeneous boundary layer near the acquisition surface. When the wave speed is not constant near the boundary, we need to localize computations and either introduce a partitioning of the acquisition surface or use wave packets as individual local data quanta, yielding a wave packets based procedure at the price of increased computational complexity with respect to a frequency box driven algorithm.
We note that by viewing wave packets as localized plane waves, our method can be connected to plane-wave and beam-wave migration [5]. Here, we can construct "beams" as reverse time continued data wave packets based on phase-space localized paraxial approximation in geodesic coordinates.
Reverse time continuation from the boundary can in principle be generalized to extended imaging using multi-source data based on downward continuation [21]. The corresponding evolution equation replacing (2.6) can be found in [12,Eq. (17)]. In this case, the evolution equation is defined in 2n − 1 dimensional extended space.