SEISMIC INVERSE SCATTERING IN THE ‘WAVE-EQUATION’ APPROACH

CHRISTIAAN C. STOLK AND MAARTEN V. DE HOOP

Abstract. Seismic data are commonly modeled by a high-frequency single scattering approximation. This amounts to a linearization in the medium coefficient about a smooth background. The discontinuities are contained in the medium perturbation. The wave solutions in the background medium admit a geometrical optics representation. Here we describe the wave propagation in the background medium by a one-way wave equation. Based on this we derive the double-square-root equation, which is a first order pseudodifferential equation, that describes the continuation of seismic data in depth. We consider the modeling operator, its adjoint and reconstruction based on this equation. If the rays in the background that are associated with the reflections due to the perturbation are nowhere horizontal, the singular part of the data is described by the solution to an inhomogeneous double-square-root equation. We derive a microlocal reconstruction equation. The main result is a characterization of the angle transform that generates the common image point gathers, and a proof that this transform contains no artifacts. Finally, pseudodifferential annihilators based on the double-square-root equation are constructed. The double-square-root equation approach is used in seismic data processing.

1. Introduction

In reflection seismology one places point sources and point receivers on the earth’s surface. The source generates acoustic waves in the subsurface, that are reflected where the mediumproperties vary discontinuously. The recorded reflections that can be observed in the data are used to reconstruct these discontinuities.

The data are commonly modeled by a high-frequency single scattering approximation. This amounts to a linearization in the medium coefficient about a smooth background. The discontinuities are contained in the medium perturbation [1]. Thus a linear operator, the modeling operator, depending on the background, that maps the perturbation to the data is obtained. Both the smooth background and the perturbation are in general unknown and have to be reconstructed jointly. In this paper we analyze this reconstruction in the wave-equation approach.

The reconstruction of the perturbation given the background is essentially done by applying the adjoint of mentioned linear map (seismic imaging). The solutions in the background medium admit a geometrical optics representation. Thus the modeling operator is Fourier integral operator [15] (for a general reference of Fourier integral operators see [6]). If the composition of adjoint and modeling operator, the normal operator, is pseudodifferential, then the position of the singularities of the perturbation are recovered by applying

The authors thank The Mathematical Sciences Research Institute for partial support, through NSF grant DMS-9810361. They also thank the members of the institute, and in particular Gunther Uhlmann, for providing a very stimulating environment during the Inverse Problems program in Fall 2001.

the adjoint to the data, and a microlocal reconstruction can be carried out. Under various assumptions on the background, concerning the presence of caustics and the geometry of the rays, results concerning the normal operator have been obtained [1, 8, 12, 10, 16] (section 2.)

In the Kirchhoff approach an approximation to the adjoint operator is constructed that depends on the background only through travel times and complex amplitudes that appear in the geometrical optics approximation. In the so called wave-equation approach [4, 3, 13] the data are downward continued, leading to data from fictitious experiments below the surface at various depths. To form the adjoint a restriction is applied to the downward continued data (imaging condition). Using a one-way wave equation (section 4) for the propagation of waves in the background, the downward continuation is described by the so called double-square-root equation.

The data is formally redundant. There exist invertible Fourier integral operators that generate a set of reconstructions [18] from which the common image point gathers are obtained. Under Beylkin’s conditions [1], in particular in the absence of caustics, this can be done by using different subsets of the data. In the presence of caustics the mentioned set can be parameterized by angle between in- and out-going rays at the image point. If the background medium is correct the reconstructions in the set should be the same. This is a criterion that is used in the reconstruction of the background (migration velocity analysis). The redundancy leads to the existence of pseudodifferential operators that annihilate the singular part of the data [18].

In the Kirchhoff approach common image point gathers parameterized by angle can be generated by a generalized Radon transform. In the presence of caustics artifacts were observed in numerical examples in [2]. By microlocal analysis of the Kirchhoff approach the presence of artifacts was shown in [17]. In section 3 an approach to suppress these artifacts is discussed.

In this paper, we consider the modeling, adjoint and reconstruction based on the doublesquare-root equation approach. The double-square-root equation is a pseudodifferential equation. If the rays in the background that are associated with the reflections due to the perturbation are nowhere horizontal, the singular part of the data is described by the solution to an inhomogeneous double-square-root equation (section 5). We derive a microlocal reconstruction equation in Proposition 6.1. The main result, Theorem 7.1 and Proposition 7.2, is a characterization of the angle transform that generates the common image point gathers, and a proof that this transform contains no artifacts. Annihilators based on the double-square-root equation are constructed in Corollary 8.1.

2. High-frequency Born modeling and imaging

We consider the scalar wave equation for acoustic waves in a constant density medium in ℝ^{n}. In preparation of the later analysis, we distinguish the vertical coordinate z ∈ ℝ from the horizontal coordinates x ∈ ℝ^{n−1} and write (x, z) ∈ ℝ^{n}. In these coordinates the scalar acoustic wave equation is given by

P u = f , P = c(x, z)^{−2} \frac{∂}{∂t}^{2} + ∑_{j=1}^{n−1} D_{x_{j}}^{2} + D_{z}^{2} , (1)

where D_{x} = −i\frac{∂}{∂x}, D_{z} = −i\frac{∂}{∂z}. The equation is considered in a time interval ]0, T [.

If c ∈ C^{∞} the solution operator of (1) propagates singularities along bicharacteristics. These are the solutions of a Hamilton system with Hamiltonian given by the principal symbol of P

P (x, z, ξ, ζ, τ) = −c(x, z)^{−2}τ^{2} + ∥ξ∥^{2} + ζ^{2}.

The Hamilton system is given by

\frac{∂(x, z, t)}{∂λ} = \frac{∂P}{∂(ξ, ζ, τ)}, \frac{∂(ξ, ζ, τ)}{∂λ} = −\frac{∂P}{∂(x, z, t)}. (2)

Its solutions will be parameterized by initial position (x_{0}, z_{0}), take-off direction α ∈ S^{n−1} and frequency τ,

x = x(x_{0}, z_{0}, α, τ, t)

and similarly for z, t, ξ, ζ; τ is invariant along the Hamilton flow. Here the evolution parameter is the time t. In section 4 we will change the evolution parameter to z, and use a similar notation to denote the bicharacteristics parameterized by z.

By Duhamel’s principle, a causal solution operator for the inhomogeneous equation (1) is given by

(3)  u(x, z, t) = ∫_{0}^{t} ∫ G(x, z, t − t_{0}, x_{0}, z_{0})f(x_{0}, z_{0}, t_{0}) dx_{0}dz_{0}dt_{0},

where G is a Fourier integral operator with canonical relation that is essentially a union of bicharacteristics. Its kernel can be written as a sum of contributions

(4) G(x, z, t, x_{0}, z_{0}) = ∑_{i} ∫_{ℝ^{N^{(i)}}} a^{(i)}(x, z, t, x_{0}, z_{0}, θ) exp[iϕ^{(i)}(x, z, x_{0}, z_{0}, t, θ)] dθ,

where the ϕ^{(i)} are non-degenerate phase functions and the a^{(i)} suitable symbols, see [6, chapter 5].

We adopt the linearized scattering approximation. The linearization is in the coefficient c around a smooth background c_{0}, c = c_{0}+δc. The perturbation δc may contain singularities. We assume that its support is contained in z > 0. The perturbation in G at the acquisition surface z = 0 is given by

(5) δG(r, 0, t, s, 0) = ∫_{ℝ^{n−1}×ℝ_{+}} ∫_{0}^{t} G(r, 0, t − t_{0}, x_{0}, z_{0}) 2c_{0}^{−3}(x_{0}, z_{0})δc(x_{0}, z_{0})

∂_{t_{0}}^{2}G(x_{0}, z_{0}, t_{0}, s, 0) dt_{0}dx_{0}dz_{0},

where both r, s ∈ ℝ^{n−1}. The singular part of δG is obtained by substituting (4) into (5). This defines the data modeling map

F = F [c_{0}] : δc ↦ R_{0}δG,

where R_{0} is the restriction defined by

R_{0} : D′_{Γ}(ℝ^{2n+1}) → D′(Y) , u(x, z, t, x_{0}, z_{0}) ↦ (R_{0}u)(x_{0}, x, t) = u(x, 0, t, x_{0}, 0)

with acquisition manifold Y a bounded open subset of ℝ^{2n−2} × ℝ_{+} that contains the range of values of (s, r, t). The restriction is defined on distributions D′_{Γ} with wavefront sets contained in an open conic set Γ such that Γ̄ does not contain points of the form (x, 0, t, x_{0}, 0, 0, ζ, 0, ζ_{0}, 0). Since τ ≠ 0 on WF(δG) the restriction is well defined. Since Y is bounded and the waves propagate with finite speed we may assume that δc is supported in a bounded open subset X of ℝ^{n−1} × ℝ_{+}. We assume that \overline{X} ∩ \{z = 0\} = ∅.

Assumption 1. There are no rays from (s, 0) to (r, 0) with travel time t such that (s, r, t) ∈ Y. For all ray pairs connecting (r, 0) via some (x, z) ∈ X to (s, 0) with total time t such that (s, r, t) ∈ Y, the rays intersect the plane z = 0 transversally at r and s.

Theorem 2.1. [15, 10] With Assumption 1 the map F is a Fourier integral operator D′(X) → D′(Y) of order (n − 1)/4 with canonical relation

(6){(x(x, z, β, τ, t_{s}), x(x, z, α, τ, t_{r}), t_{s} + t_{r}, ξ(x, z, β, τ, t_{s}), ξ(x, z, α, τ, t_{r}), τ; x, z, ξ, ζ) |

t_{s}, t_{r} > 0, z(x, z, β, τ, t_{s}) = z(x, z, α, τ, t_{r}) = 0, (ξ, ζ) = −τc_{0}(x, z)^{−1}(α + β),

(x, z, α, β, τ) ∈ subset of X × (S^{n−1})^{2} × ℝ\0}.

Assumption 1 is microlocal. One can identify the set of points (s, r, t, σ, ρ, τ) ∈ T^{∗}Y \0 where this assumption is violated. If the symbol ψ = ψ(s, r, t, σ, ρ, τ) vanishes on a neighborhood of this set, then the composition ψF of the pseudodifferential cutoff ψ = ψ(s, r, t, D_{s}, D_{r}, D_{t}) with F is a Fourier integral operator as in the theorem.

We assume that ψ is as before and in addition vanishes outside Y. To image the singularities of δc from the data we consider the adjoint F^{∗}ψ, which is a Fourier integral operator also.

Assumption 2. [7] The projection of the canonical relation (6) on T^{∗}Y \0 is an embedding.

Since (6) is a canonical relation that projects submersively on the subsurface variables (x, z, ξ, ζ), the projection of (6) on T^{∗}Y \0 is immersive [9, 25.3.6]. Therefore only the injectivity in the assumption needs to be verified [10].

The following theorem describes the reconstruction of δc modulo a pseudodifferential operator with principal symbol that is nonzero at (x, z, ξ, ζ) whenever there is a point (s, r, t, σ, ρ, τ; x, z, ξ, ζ) in the canonical relation (6) with (s, r, t, σ, ρ, τ) in the support of ψ (i.e. whenever there is illumination or insonification).

Theorem 2.2. With Assumption 2 the operator F^{∗}ψF is pseudodifferential of order n − 1.

3. Generalized Radon transform in scattering angle

Consider the projection of the canonical relation (6) on the (x, z, s, r, τ) variables. Where this projection is locally diffeomorphic, the canonical relation (6) can be described by a phase function of the form

τ(T^{(m)}(x, z, s, r) − t)

where T^{(m)} is the value of the time variable in (6). There is a set {T^{(m)}}_{m∈M} that describes the canonical relation except for a neighborhood of the subset of the canonical relation where mentioned projection is degenerate. Each T^{(m)} is defined on a subset D^{(m)} of ℝ_{(x,z,s,r)}^{3n−2}. We define F^{(m)} to be a contribution to F with phase function given by T^{(m)}, such that on a subset of the canonical relation where the projection is nondegenerate F is given microlocally by ∑_{m∈M} F^{(m)}.

We can use (x, z, ξ, ζ) ∈ T^{∗}(ℝ^{n−1}×ℝ_{+})\0 as local coordinates on the canonical relation (6). In addition, we need to parameterize the subsets of the canonical relation given by (x, z, ξ, ζ) = constant; we denote such parameters by e. The canonical relation (6) was parameterized by (x, z, α, β, τ). We relate (x, z, ξ, ζ, e) by a coordinate transformation to (x, z, α, β, τ): A suitable choice when α ≠ β is the scattering angles given by

(7)  e(x, z, α, β) = (arccos(α · β), \frac{−α + β}{2 sin(arccos(α · β)/2)}) ∈]0, π[×S^{n−2}.

The migration dip ν is defined as

(8)  ν(α, β) = \frac{α + β}{∥α + β∥} ∈ S^{n−1}.

On D^{(m)} there is a map (x, z, α, β) ↦ (x, z, s, r). We define e^{(m)} = e^{(m)}(x, z, s, r) as the composition of e with the inverse of this map. Likewise, we define ν^{(m)} = ν^{(m)}(x, z, s, r).

We define the generalized Radon transform in scattering angle or the Kirchhoff angle transform via a restriction in F^{∗} of the mapping e^{(m)} to a prescribed value e, i.e. the distribution kernel of each contribution F^{(m)}^{∗} is multiplied by δ(e − e^{(m)}(x, z, s, r)) . Its kernel is given by

(9) L(x, z, e, r, s, t) = ∑_{m∈M}(2π)^{−(n−1)} ∫ \overline{A^{(m)}(x, z, s, r, τ)}e^{iΦ^{(m)}(x,z,e,s,r,t,ε,τ)} dτdε,

where A^{(m)} is a symbol for the m-th contribution to F, supported on D^{(m)}, and

Φ^{(m)}(x, z, e, s, r, t, ε, τ) = τ(T^{(m)}(x, z, s, r) − t) + ⟨ε, e − e^{(m)}(x, z, s, r)⟩.

Here, ε is the cotangent vector corresponding to e. Let ψ_{L} = ψ_{L}(D_{s}, D_{r}, D_{t}) be a pseudodifferential cutoff such that ψ(σ, ρ, τ) = 0 on a conic neighborhood of τ = 0. Then [17] ψ_{L}L is a Fourier integral operator with canonical relation

(10) ∪_{m∈M} {(x, z, e^{(m)}(x, z, s, r), ξ^{(m)}(x, z, s, r, τ, ε), ζ^{(m)}(x, z, s, r, τ, ε), ε;

s, r, T^{(m)}(x, z, s, r), ρ^{(m)}(x, z, s, r, τ, ε), σ^{(m)}(x, z, s, r, τ, ε), τ) |
(x, z, s, r) ∈ D^{(m)}, ε ∈ ℝ^{n−1}, τ ∈ ℝ\0}

where

(11) ξ^{(m)}(x, z, s, r, τ, ε) = ∂_{x}Φ^{(m)} = τ∂_{x}T^{(m)}(x, z, s, r) − ⟨ε, ∂_{x}e^{(m)}(x, z, s, r)⟩,

and likewise expressions for ζ^{(m)}, σ^{(m)} and ρ^{(m)}.

Let d be the Born modeled data. To reveal any artifacts generated by L, i.e. singularities in Ld at positions not corresponding to an element of WF(δc), we consider the composition L F. This composition is equal to the sum of a smooth e-family of pseudodifferential operators and, in general, a non-microlocal operator the wavefront set of which contains no elements with ε = 0 [17, theorem 6.1].

We discuss the practical considerations of the suppression of artifacts generated by L by modification of transform (9). The correct way to remove artifacts from Ld is by a Fourier cut-off in ε/∥(ξ, ζ)∥ about 0. For bandlimited data, this is done approximately by local averaging in e at every (x, z). Such procedure has been applied to both synthetic and real data examples in [2]. The artifacts were not fully suppressed. The following modification of the kernel of L (9) was applied,

L_{χ}(x, z, e, r, s, t) = ∑_{m∈M}(2π)^{−(n−1)} ∫ \overline{A^{(m)}(x, z, s, r, τ)}χ(x, z, ν^{(m)}(x, z, s, r))

× e^{iΦ^{(m)}(x,z,e,s,r,t,ε,τ)} dτdε,

where χ(x, z, ν) is a smooth cutoff function on ℝ^{n} × S^{n−1}. Observe that ν^{(m)} is the direction of \frac{∂T^{(m)}}{∂(x,z)} .

Remark 3.1. The transform L_{χ} restricts the wavefront set of operator L. If WF(δc) is contained in {(x, z, λ(α_{x}(x, z), α_{z}(x, z))) | (x, z, λ) ∈ ℝ^{n−1} × ℝ_{+} × ℝ}, where (α_{x}(x, z), α_{z}(x, z)) is a smooth covector field on ℝ^{n−1} × ℝ_{+}, there is a χ with a small support in ν such that L_{χ} generates the true image. Artifacts with singular direction (ξ^{(m)}, ζ^{(m)}) such that the direction of \frac{∂T^{(m)}}{∂(x,z)} is outside the support of χ, are suppressed.

The main difficulty arises when the background medium is not (accurately) known. Without knowledge of the background medium there is no criterion to distinguish artifacts from the true image. The approach following the double-square-root equation to be introduced below does not generate artifacts. This is of particular importance in the development of tomographic methods for reconstructing the background medium.

4. The one-way wave or single-square-root equation

In this section we discuss the problem of solving the wave equation by evolution in the vertical, z direction. This problem is in general not well posed, but microlocal solutions can be obtained.

Consider the wave equation rewritten as a first-order system

(12)  \frac{∂}{∂z} ( u_{\frac{∂u}{∂z}}) = (\begin{matrix}0 & 1 \\ −A(x, z, D_{x}, D_{t}) & 0 \end{matrix}) ( u_{\frac{∂u}{∂z}}) + (\begin{matrix}0 \\ f \end{matrix}) ,

where A(x, z, ξ, τ) = c_{0}(x, z)^{−2}τ^{2} − ∥ξ∥^{2}. Microlocally, away from the zeroes of A(x, z, ξ, τ), this system can be transformed into diagonal form modulo a smoothing operator [20]. There is a family of pseudodifferential operator matrices Q(z) = Q(x, z, D_{x}, D_{t}) such that, microlocally,

(\begin{matrix}u_{+} \\ u_{−} \end{matrix}) = Q(z) ( u_{\frac{∂u}{∂z}}) ,  (\begin{matrix}f_{+} \\ f_{−} \end{matrix}) = Q(z) (\begin{matrix}0 \\ f \end{matrix}) ,

satisfy the one-way wave or single-square-root (SSR) equations

(13)  (\frac{∂}{∂z} ± iB_{±}(x, z, D_{x}, D_{t})) u_{±} = f_{±},

see [20]. The principal symbol b of the B_{±} are given by b(x, z, ξ, τ) = \sqrt{A(x, z, ξ, τ)} = τ\sqrt{\frac{1}{c_{0}(x,z)^{2}} − τ^{−2}∥ξ∥^{2}}. For (x, t, ξ, τ) such that the symbol B_{±} is real, the equation is of hyperbolic type, corresponding (microlocally) to propagating waves. To describe the associated bicharacteristics the Hamiltonian P can be used, as well as ζ ∓b(x, z, ξ, τ). The solution of a one-way wave equation describes the propagation of singularities along rays in intervals where ±\frac{∂z}{∂t} > 0. The case where the symbol B_{±} is imaginary corresponds to either evanescent waves, or waves that blow up like in a backward heat equation.

We choose the normalization of Q(z), such that (13) is selfadjoint microlocally where the symbol is real,

u = Q_{+}^{∗}u_{+} + Q_{−}^{∗}u_{−},
f_{±} = ± \frac{1}{2}iQ_{±}f,

where Q_{±} = Q_{±}(z) = Q_{±}(x, z, D_{x}, D_{t}) are z-families of pseudodifferential operators with principal symbols τ^{−1/2}(\frac{1}{c_{0}(x,z)^{2}} − τ^{−2}∥ξ∥^{2})^{−1/4}.

At the zeroes of A(x, z, ξ, τ) the operators B_{±}(z), Q_{±}(z) are not yet defined. To this end, we regularize the problem by replacing A(x, z, ξ, τ) in (12) by

(14)  c_{0}(x, z)^{−2}τ^{2} − ∥ξ∥^{2} − iτ^{2}ϕ(x, z, ξ, τ),

where ϕ(x, z, ξ, τ) is positive, small, homogeneous of order 0 in (ξ, τ) and is supported on a small neighborhood of the set of zeroes of A(x, z, ξ, τ). With this modification the operators B_{±}, Q_{±} are z-families of pseudodifferential operators, defined on the entire cotangent space T^{∗}ℝ_{(x,t)}^{n}\0. With this choice of sign for the regularizing imaginary term there is a well defined solution operator G_{−}(z, z_{0}), z_{0} > z, of the initial value problem for u_{−} given by (13) with f_{−} = 0, see [22, theorem XI.2.1]. The adjoint G_{−}(z, z_{0})^{∗} describes the propagation from z to z_{0} of (13), regularized in accordance with (14), but with opposite sign of the imaginary part. The operator G_{−} is a Fourier integral operator with complex phase [11], [9, chapter XXV], [22, chapters X and XI]. By Duhamel’s principle a microlocal solution operator for the inhomogeneous equation is given by

(15)  u_{−}(·, z) = ∫_{−∞}^{z} G_{−}(z, z_{0})f_{−}(·, z_{0}) dz_{0}.

Microlocally −\frac{1}{2}iQ_{−}^{∗}(z)G_{−}(z, z_{0})Q_{−}(z_{0}) is equal to the Green’s function for singularities propagating along bicharacteristics (cf. (2)) with −\frac{∂z}{∂t} > ϵ, for some ϵ > 0 that depends on c_{0} and the support of ϕ in (14).

The operator G_{−} propagates singularities at (x_{0}, ξ_{0}, τ, z_{0}) along the bicharacteristics for z in an interval ]Z(x_{0}, ξ_{0}, τ, z_{0}), z_{0}], which is the maximal interval such that the regularized symbol b is real valued. As a consequence the bicharacteristic in this interval is nowhere horizontal. From now we use z as the evolution parameter for bicharacteristics, and denote them as

(16) (x(x_{0}, z_{0}, ξ_{0}, τ, z), z, t(x_{0}, z_{0}, ξ_{0}, τ, z), ξ(x_{0}, z_{0}, ξ_{0}, τ, z), b(x(x_{0}, z_{0}, ξ_{0}, τ, z), z, ξ(x_{0}, z_{0}, ξ_{0}, τ, z), τ), τ).

Remark 4.1. The symbols of the operators B_{±}(z), Q_{±}(z) can be written as an infinite sum of elementary symbols, ∑_{i=1}^{∞} v_{i}(x, z, τ)w_{i}(z, ξ, τ) say, that is rapidly converging where they are smooth, i.e. away from the set of zeroes of A(x, z, ξ, τ), in accordance with [21, equation (2.1.11)]. The elementary symbols correspond to multiplications either in horizontal position space (v_{i}) or in horizontal wavenumber (Fourier) space (w_{i}). This is exploited in fast numerical solvers of (13), see [5].

5. Modeling, the double-square-root equation

We show that the Born modeling operator can be written, modulo smoothing terms, in terms of the solution operator to the double-square-root (DSR) equation (21) below. We assume that the rays that connect source and receiver to a scattering point in X have nowhere horizontal tangent directions.

Assumption 3. (DSR assumption) If (x, z) ∈ X and α, β ∈ S^{n−1}, t_{s}, t_{r} > 0 depending on (x, z, α, β) are such that z(x, z, β, τ, t_{s}) = z(x, z, α, τ, t_{r}) = 0 and (x(x, z, β, τ, t_{s}), x(x, z, α, τ, t_{r}), t_{s} + t_{r}) ∈ Y (cf. (6)), then

\frac{∂z}{∂t} (x, z, β, τ, t) < −ϵ, t ∈ [0, t_{s}],
\frac{∂z}{∂t} (x, z, α, τ, t) < −ϵ, t ∈ [0, t_{r}],

where ϵ > 0 was introduced below (15).

This assumption is stronger than Assumption 1. This assumption is microlocal, and, given the background medium, a pseudodifferential cutoff can be applied to the data to remove microlocally the part of the data where Assumption 3 is violated.

Under Assumption 3 and the assumption that δc = 0 on a neighborhood of z = 0, the singular part of the Born data is unchanged when G in (5) is replaced by −\frac{1}{2}iQ_{−}^{∗}(x, z, D_{x}, D_{t})G_{−}Q_{−}(x, z, D_{x}, D_{t}). Define the operator H(z, z_{0}), z < z_{0} by

(17) (H(z, z_{0}))(s, r, t, s_{0}, r_{0}, t_{0}) =

∫_{ℝ}(G_{−}(z, z_{0}))(s, t − t_{0} − t′, s_{0}, 0)(G_{−}(z, z_{0}))(r, t′, r_{0}, 0) dt′.

Here (G_{−}(z, z_{0}))(r, t′, r_{0}, 0) denotes the distribution kernel of G_{−}(z, z_{0}), and similarly for H(z, z_{0}). Define the maps I_{1}, I_{2} by

(18)  I_{1} : D′(ℝ^{n}) → D′(ℝ^{2n−1}) : u(x, z) ↦ δ(r − s)u(\frac{r+s}{2}, z),
(19)  I_{2} : D′(ℝ^{2n−1}) → D′(ℝ^{2n}) : u(r, s, z) ↦ δ(t)u(\frac{r+s}{2}, z).

The operators G, G_{−}, Q_{−} all are of convolution type in the time variable. It follows that the singular part of the Born approximated data (5) is given by

(20) F δc = ∫_{0}^{z_{max}} Q_{−,s}^{∗}(0)Q_{−,r}^{∗}(0)H(0, z)Q_{−,s}(z)Q_{−,r}(z) \frac{1}{2}D_{t}^{2}(I_{2}I_{1} c_{0}^{−3}δc)(·, z) dz.

where Q_{−,s}(z) is short for Q_{−}(s, z, D_{s}, D_{t}) and similarly for Q_{−,r}(z). Here z_{max} is the maximum depth illuminated from acquisition manifold Y given the background medium c_{0}.

Define the inhomogeneous double-square-root (DSR) equation by

(21)  (\frac{∂}{∂z} − iB_{−}(s, z, D_{s}, D_{t}) − iB_{−}(r, z, D_{r}, D_{t}))u = g(s, r, t, z).

It follows from the definition of G_{−}(z, z_{0}) that the operator H(z, z_{0}) is a solution operator for the Cauchy initial value problem for the (regularized) DSR equation. The Born approximated data is modulo smooth functions equal to the solution of the DSR equation at z = 0 with inhomogeneous term g(s, r, t, z) = \frac{1}{2}D_{t}^{2}(I_{2}I_{1} c_{0}^{−3}δc)(s, r, t, z) .

The operator H(z, z_{0}) propagates singularities along the bicharacteristics for z in an interval ]Z_{min}(s_{0}, r_{0}, σ_{0}, ρ_{0}, τ), z_{0}], which is the intersection of the one-way intervals associated with the source and receiver bicharacteristics, in the notation of (16),

(22) (x(s_{0}, z_{0}, σ_{0}, τ, z), x(r_{0}, z_{0}, ρ_{0}, τ, z), t_{0} + t(s_{0}, z_{0}, σ_{0}, τ, z) + t(r_{0}, z_{0}, ρ_{0}, τ, z), z, ξ(s_{0}, z_{0}, σ_{0}, τ, z), ξ(r_{0}, z_{0}, ρ_{0}, τ, z), τ, Γ(x(s_{0}, z_{0}, σ_{0}, τ, z), x(r_{0}, z_{0}, ρ_{0}, τ, z), ξ(s_{0}, z_{0}, σ_{0}, τ, z), ξ(r_{0}, z_{0}, ρ_{0}, τ, z), τ, z)),

where

(23)  Γ(s, r, σ, ρ, τ, z) = b(s, z, σ, τ) + b(r, z, ρ, τ).

6. Double-square-root reconstruction

Let N := F^{∗}ψF, where ψ = ψ(s, r, t, D_{s}, D_{r}, D_{t}) is a suitable pseudodifferential cutoff as in section 2. By Theorem 2.1, with Assumption 2, N = N(x, z, D_{x}, D_{z}) is a pseudodifferential operator. Therefore

(24)  N(x, z, D_{x}, D_{z})δc = F^{∗}ψ(s, r, t, D_{s}, D_{r}, D_{t})d,

and δc can be reconstructed microlocally where the principal symbol of N is non-zero. Assumption 3 implies the injectivity in Assumption 2 and is hence stronger than Assumption 2 by the remark below Assumption 2. We present a reconstruction formula similar to (24) based on the DSR Born modeling (20). This leads, again, to reconstruction modulo a pseudodifferential operator for which an explicit expression is given.

The adjoint operator H(0, z)^{∗} propagates the data downward and backward in time. We consider H(0, z)^{∗}d as a function of (s, r, t, z). The operator H(0, z)^{∗} is microlocally a Fourier integral operator with real phase and canonical relation which follows from (22) with z_{0} = 0, for z in an interval [0, Z_{max}(s_{0}, r_{0}, σ_{0}, ρ_{0}, τ)[

(25) {(x(s_{0}, 0, σ_{0}, τ, z), x(r_{0}, 0, ρ_{0}, τ, z), t_{0} + t(s_{0}, 0, σ_{0}, τ, z) + t(r_{0}, 0, ρ_{0}, τ, z), z, ξ(s_{0}, 0, σ_{0}, τ, z), ξ(r_{0}, 0, ρ_{0}, τ, z), τ, Γ(s_{0}, r_{0}, σ_{0}, ρ_{0}, τ, z); s_{0}, r_{0}, t_{0}, σ_{0}, ρ_{0}, τ) |

(s_{0}, r_{0}, t_{0}, σ_{0}, ρ_{0}, τ) ∈ T^{∗}ℝ^{2n−1}\0, z ∈ [0, Z_{max}(s_{0}, r_{0}, σ_{0}, ρ_{0}, τ)[}.

The adjoint of the operator I_{2} (cf. (19)) is given by the restriction R_{2} defined by

g(r, s, t, z) ↦ (R_{2}g)(r, s, z) = g(r, s, 0, z).

If u = u(s, r, z) then Ku = (Ku)(s, r, t) will be defined by

Ku = ∫_{−∞}^{0} H(0, z)I_{2}u(·, z) dz.

Let ψ be a pseudodifferential cutoff in the (s, r, t) variables supported in

{(s_{0}, r_{0}, t_{0}, σ_{0}, ρ_{0}, τ) | t_{0} ∈

[0, −t(s_{0}, 0, σ_{0}, τ, Z_{max}(s_{0}, r_{0}, σ_{0}, ρ_{0}, τ)) − t(r_{0}, 0, ρ_{0}, τ, Z_{max}(s_{0}, r_{0}, σ_{0}, ρ_{0}, τ))[}.

Then ψK is a Fourier integral operator with real phase. The adjoint K^{∗} of K follows from the equality

⟨ψd, ∫_{−∞}^{0} H(0, z)I_{2}u(·, z) dz⟩_{(s,r,t)} = ∫_{−∞}^{0}⟨R_{2}H(0, z)^{∗}ψd, u(·, z)⟩_{(s,r)} dz

and is given by

K^{∗}ψ = R_{2}H(0, z)^{∗}ψ.

Since \frac{∂t}{∂z} < 0, t depends strictly monotone on z for each (s_{0}, r_{0}, t_{0}, σ_{0}, ρ_{0}, τ). The DSR bicharacteristic with initial values (s_{0}, r_{0}, t_{0}, σ_{0}, ρ_{0}, τ) hence intersects the t = 0 hyperplane at most once, and the intersection is transversal. From this, it follows that the composition K^{∗}ψ is a Fourier integral operators with a locally invertible canonical relation.

The canonical relation maps points (s_{0}, r_{0}, t_{0}, σ_{0}, ρ_{0}, τ) in a subset of the cotangent acquisition space T^{∗}ℝ^{2n−1}\0 to (r, s, z, σ, ρ, ζ) in a subset of T^{∗}ℝ^{2n−1}\0. This map converts time to depth. We define the symbol Ψ(s, r, z, σ, ρ, ζ) as the pull back of ψ(s_{0}, r_{0}, t_{0}, σ_{0}, ρ_{0}, τ) under the inverse of this map. Starting from the Born modeling (20), we find the following reconstruction proposition.

Proposition 6.1. There are pseudodifferential operators Φ = Φ(x, z, D_{x}, D_{z}) of order n − 1 with principal symbol

(26)  Φ(x, z, ξ, ζ) = ∫_{ℝ^{n−1}} Ψ(x, x, z, \frac{1}{2}ξ − θ, \frac{1}{2}ξ + θ, ζ) dθ,

and Ξ(z) = Ξ(s, r, t, D_{s}, D_{r}, D_{t}, z) of order 0 with principal symbol

Ξ(s, r, t, σ, ρ, τ, z) = |\frac{∂Γ}{∂τ}(s, r, σ, ρ, τ, z)|

= c_{0}(s, z)^{−2}(c_{0}(s, z)^{−2} − τ^{−2}∥σ∥^{2})^{−1/2} + c_{0}(r, z)^{−2}(c_{0}(r, z)^{−2} − τ^{−2}∥ρ∥^{2})^{−1/2},

such that

(27) Φ(x, z, D_{x}, D_{z})δc

= 2c_{0}^{3}R_{1}R_{2}Ξ(z)Q_{−,s}^{∗}(z)^{−1}Q_{−,r}^{∗}(z)^{−1}H(0, z)^{∗}Q_{−,s}(0)^{−1}Q_{−,r}(0)^{−1}D_{t}^{−2}ψd,

where d = F δc is the Born modeled data.

Proof. We calculate microlocally the principal symbol of K^{∗}K. The kernel of the operator H(0, z) has microlocally an oscillatory integral representation with a phase function associated with generating function, S = S(z, s, r, t, y_{0I}, η_{0J}) say, where y_{0} = (s_{0}, r_{0}, t_{0}) and η_{0} is the corresponding cotangent vector, and {I, J} is a partition of {1, …, 2n − 1},

H(0, z)(y_{0}, s, r, t) = (2π)^{−(2n−1+|I|)/2}

× ∫ A(z, s, r, t, y_{0}, η_{0J})e^{i(S(z,s,r,t,y_{0I},η_{0J})−⟨η_{0J} ,y_{0J} ⟩)} dη_{0J}, The adjoint H(0, z)^{∗} has amplitude \overline{A(z, s, r, t, y_{0}, η_{0J})} and phase −S(z, s, r, t, y_{0I}, η_{0J}) + ⟨η_{0J}, y_{0J}⟩. Hence, the kernel of the composition H(0, z)^{∗}H(0, z) has the oscillatory integral representation

(2π)^{−(2n−1)} ∫ \overline{A(z, s′, r′, t′, y_{0}, η_{0J})}A(z, s, r, t, y_{0}, η_{0J})

× e^{i[−S(z,s′,r′,t′,y_{0I},η_{0J})+S(z,s,r,t,y_{0I},η_{0J})]} dy_{0I}dη_{0J}.

We expand the phase in a Taylor series around (s′, r′, t′) = (s, r, t) and identify the gradient

− \frac{∂S}{∂(s, r, t)}(z, s, r, t, y_{0I}, η_{0J}) =

(σ(z, s, r, t, y_{0I}, η_{0J}), ρ(z, s, r, t, y_{0I}, η_{0J}), τ(z, s, r, t, y_{0I}, η_{0J})).

Applying a change of variables, (y_{0I}, η_{0J}) ↦ (σ, ρ, τ), the phase takes the form

⟨(σ, ρ, τ), (s′ − s, r′ − r, t′ − t)⟩.

Since H(0, z)^{∗}H(0, z) is a pseudodifferential operator with symbol 1, microlocally in the support of the cutoff ψ, we conclude that the principal part a of the amplitude A satisfies

|a(z, s, r, t, y_{0}, η_{0J})| = |\frac{∂(σ, ρ, τ)}{∂(y_{0I}, η_{0J})}|^{1/2} .

The kernel of operator K has an oscillatory integral representation similar to the one of H(0, z),

K(y_{0}, s, r, z) = (2π)^{−(2n−1+|I|)/2} ∫ A(z, s, r, 0, y_{0}, η_{0J})e^{i(S(z,s,r,0,y_{0I},η_{0J})−⟨η_{0J} ,y_{0J} ⟩)} dη_{0J}

(we have applied I_{2} at z). It follows that the composition K^{∗}K, carrying out an analysis similar to the one for H(0, z)^{∗}H(0, z), is a pseudodifferential operator, microlocally. Its amplitude has principal part

|\frac{∂(ζ, σ, ρ)}{∂(y_{0I}, η_{0J})}|^{−1} |\frac{∂(σ, ρ, τ)}{∂(y_{0I}, η_{0J})}| .

For fixed (s, r, σ, ρ, z) the map τ ↦ Γ(s, r, σ, ρ, τ, z) is invertible on a set given by |τ| sufficiently large. This map will be denoted by Γ^{−1} = Γ^{−1}(s, r, z, σ, ρ, ζ). It follows that the principal part of K^{∗}K is microlocally given by

|\frac{∂Γ}{∂τ}(s, r, σ, ρ, Γ^{−1}(s, r, z, σ, ρ, ζ), z)|^{−1} .

We finally consider the composition of operators R_{1}ΨI_{1}. Its kernel has an oscillatory integral representation,

(2π)^{−(2n−1)} ∫_{ℝ^{2n−1}} Ψ(x, x, z, σ, ρ, ζ) e^{i⟨(x,x,z)−(x′,x′,z′),(σ,ρ,ζ)⟩} dρdσdζ =

(2π)^{−(2n−1)} ∫_{ℝ^{n}} ∫_{ℝ^{n−1}} Ψ(x, x, z, \frac{1}{2}ξ − θ, \frac{1}{2}ξ + θ, ζ) dθ e^{i(⟨(x−x′),ξ⟩+⟨(z−z′),ζ⟩)} dξdζ, upon changing variables of integration, σ = \frac{1}{2}ξ − θ, ρ = \frac{1}{2}ξ + θ. The domain of the θ integral is bounded depending on (ξ, ζ), since Ψ is a cutoff in (s, r, t, σ, ρ, τ). Hence R_{1}ΨI_{1} is a pseudodifferential operator of order n − 1 with principal symbol (26).

Applying the above results to the expression (20) for the Born modeled data leads to the statement of the proposition. □

Remark 6.2. Note that in (27) all the operators on the right hand side, except the downward continuation H, act at depth z only. On the contrary, the operator Φ(x, z, D_{x}, D_{z}) on the left hand side depends on the Hamiltonian flow associated with the background medium in the depth interval [0, z]. In the usual wave-equation imaging algorithms it is hence straightforward to include the pseudodifferential factors on the right hand side of (27). To account for the operator Φ(x, z, D_{x}, D_{z}) on the left hand side one requires an additional ray computation.

Remark 6.3. Depending on the background medium, the reconstruction can also be done using data on a submanifold Y′ of Y. Let R′ be the restriction of a function on Y to Y′, so that the forward map for this case is given by R′F. In suitable local coordinates (y′, y′′) on Y such that y′′ = 0 defines Y′, the adjoint I′ of R′ is given by the map (I′f)(y′, y′′) = f(y′)δ(y′′). Conditions such that F^{∗}I′ψ′R′F is pseudodifferential are given in [12], where ψ′ is a suitable pseudodifferential cutoff. Reconstruction modulo a pseudodifferential operator is done in this case by first applying the map I′ to the data, and then applying the previous procedure. Applying I′ to the data simply means adding zeroes where there is no data in Y.

7. The wave-equation angle transform

We define the wave-equation angle transform A_{WE} by the following integral of the downward continued data, H(0, z)^{∗}d,

(28)  (A_{WE}d)(x, z, p) = ∫_{ℝ^{n−1}}(H(0, z)^{∗}ψd)(x − \frac{h}{2}, x + \frac{h}{2}, ph)χ(x, z, h) dh,

(cf. [14]), where h ↦ χ(x, z, h) is a compactly supported cutoff function the support of which contains h = 0.

Theorem 7.1. Let C_{0} be an upper bound for c_{0}. Assume that

(29)  |p| < p_{max} < \frac{1}{2}C_{0}^{−1}.

Then A_{WE} is a Fourier integral operator such that A_{WE}F is a smooth p-family of pseudodifferential operators in (x, z). Let C_{1} be an upper bound for \frac{∂c_{0}^{−2}}{∂x}. If in addition the function h ↦ χ(x, z, h) is supported in B(0, R), where R depends on C_{1} and ψ, then the canonical relation of A_{WE} corresponds to an invertible map from a subset of T^{∗}ℝ_{(s,r,t)}^{2n−1} to a subset of T^{∗}ℝ_{(x,z,p)}^{2n−1} that has nonempty intersection with the set θ = 0 (where θ is the p-covector).

Proof. The map d ↦ H(0, z)^{∗}ψd is a Fourier integral operator with canonical relation given in (25).

The Schwarz kernel of the map H(0, z)^{∗}ψd ↦ A_{WE}d equals

δ(x − \frac{s+r}{2}) δ(p(r − s) − t) δ(z − z′) χ(x, z, r − s)

= (2π)^{−n−1} ∫ e^{i(⟨ξ,x− \frac{s+r}{2} ⟩+τ(p(r−s)−t)+ζ(z−z′))} dξ dτ dζ.

It is a Fourier integral operator with canonical relation that is contained in T^{∗}ℝ_{(x,z,p)}^{2n−1}\0 × T^{∗}ℝ_{(s,r,t,z)}^{2n}\0 and given by

(30) \{(\frac{s+r}{2}, z, p, ξ, ζ, (r − s)τ; s, r, p(r − s), z, \frac{ξ}{2} + pτ, \frac{ξ}{2} − pτ, τ, ζ) |

(s, r, z, p, ξ, ζ, τ) ∈ (subset of )ℝ^{4n−1}}.

This canonical relation can be parameterized by the coordinates of T^{∗}ℝ_{(s,r,t,z)}^{2n}\0 except t, that is (s, r, z, σ, ρ, τ, ζ). The projection of (30) on T^{∗}ℝ_{(s,r,t,z)}^{2n}\0 is a hypersurface defined by

(31)  t = ⟨\frac{σ − ρ}{2τ}, (r − s)⟩ .

The canonical relation of the map d ↦ H(0, z)^{∗}ψd, considered as a function of (s, r, t, z), is parameterized by (s_{0}, r_{0}, t_{0}, σ_{0}, ρ_{0}, τ, z). The canonical relation (25) is time translation invariant and the line in T^{∗}ℝ_{(s,r,t,z)}^{2n}\0 parameterized by t_{0} for fixed (s_{0}, r_{0}, σ_{0}, ρ_{0}, τ) intersects the hypersurface (31) transversally. It follows that the composition of the canonical relations (25) and (30) is transversal. The composition is parameterized by (s_{0}, r_{0}, σ_{0}, ρ_{0}, τ, z). It follows that A_{WE} is a Fourier integral operator.

The composition H(0, z)^{∗}ψF, that maps δc = δc(x, z′) to the downward continued data as a function of (s, r, t, z), is a Fourier integral operator with canonical relation

(32) {(x(x, z′, σ′, τ, z), x(x, z′, ρ′, τ, z), t(x, z′, σ′, τ, z) + t(x, z′, ρ′, τ, z), z, ξ(x, z′, σ′, τ, z), ξ(x, z′, ρ′, τ, z), τ, ζ; x, z′, σ′ + ρ′, ζ′) |

(x, z′, σ′, ρ′, ζ′, z) ∈ a subset of ℝ^{3n}, τ such that ζ′ = Γ(x, x, σ′, ρ′, τ, z′), ζ = Γ(x(x, z′, σ′, τ, z), x(x, z′, ρ′, τ, z), ξ(x, z′, σ′, τ, z), ξ(x, z′, ρ′, τ, z), τ, z)}.

The propagation of singularities upward by F and downward by H(0, z)^{∗} is along the same DSR bicharacteristics.

We show that the composition A_{WE}F is a Fourier integral with canonical relation contained in

(33) {(x, z, p, ξ, ζ, 0; x, z, ξ, ζ) | (x, z, ξ, ζ) ∈ T^{∗}ℝ_{(x,z)}^{n}\0, p ∈] − p_{max}, p_{max}[}.

From that it follows that A_{WE}F is a p-family of pseudodifferential operators. The projection of (32) on T^{∗}ℝ_{(s,r,t,z)}^{2n} intersects the hypersurface (31) at z′ = z, since then r − s = 0 and t = 0, leading to elements in (33). Since singularities propagate with speed less than C_{0}, if (v_{s}, v_{r}, v_{t}, v_{z}, v_{σ}, v_{ρ}, 0, v_{ζ}) is a tangent vector to the DSR bicharacteristic, then \frac{v_{s}−v_{r}}{v_{t}} ≤ 2C_{0}. Therefore, by (29), the composition of (32) with (30) is transversal and contains no elements outside (33).

The projection of the canonical relation of A_{WE} on the second component T^{∗}ℝ_{(s,r,t)}^{2n−1}\0 is invertible if each DSR bicharacteristic with initial values (s_{0}, r_{0}, t_{0}, σ_{0}, ρ_{0}, τ), parameterized by z, intersects the hypersurface (31) at most once and transversally. Let p(z) denote \frac{σ−ρ}{2τ} along a certain DSR bicharacteristic and let t(z) denote the time and h(z) denote the value of r − s. The elements of the canonical relation of A_{WE} correspond to solutions of t(z) − ⟨p(z), h(z)⟩ = 0. To estimate the derivative of the left hand side we observe that

\frac{∂t}{∂z} − ⟨p(z), \frac{∂h}{∂z}(z)⟩ < −ϵ_{0}

for some ϵ_{0} > 0 depending on the cutoff ψ (or ϕ in (14)) and on the value p_{max}. Since

\frac{∂ξ}{∂z} = \frac{∂b}{∂x} = −\frac{τ}{\sqrt{c_{0}^{−2} − τ^{−2}∥ξ∥^{2}}} \frac{∂c_{0}^{−2}}{∂x} ,

there is a bound on \frac{∂p}{∂z} in terms of C_{1} and on ψ (to bound the square root from below). It follows that for some ϵ_{1} < ϵ_{0}

|⟨\frac{∂p}{∂z}(z), h(z)⟩| < ϵ_{1} < ϵ_{0}

if ∥h∥ < C_{2}C_{1}^{−1} for some constant C_{2} depending on ψ and p_{max}. This implies that the function z ↦ t(z) − ⟨p(z), h(z)⟩ is monotone. Hence the projection of the canonical relation of A_{WE} on T^{∗}ℝ_{(s,r,t)}^{2n−1}\0 is invertible. It follows from the above reasoning that the projection on T^{∗}ℝ_{(x,z,p)}^{2n−1}\0 is invertible as well. This establishes the last statement of the theorem. □

To conclude this section we determine, at the principal symbol level, the modification of (28) that leads to microlocal reconstruction.

Proposition 7.2. Define Ã_{WE} by

(Ã_{WE}d)(x, z, p) = ∫_{ℝ^{n−1}} χ(x, z, h)2c_{0}(x, z)^{3}

×(Ξ(z)Q_{−,s}^{∗}(z)^{−1}Q_{−,r}^{∗}(z)^{−1}H(0, z)^{∗}Q_{−,s}(0)^{−1}Q_{−,r}(0)^{−1}D_{t}^{−2}ψd)(x− \frac{h}{2}, x+ \frac{h}{2}, ph) dh.

Suppose that χ(x, z, 0) = 1 and h ↦ χ(x, z, h) is supported in B(0, R) (cf. Theorem 7.1), then Ã_{WE} is an invertible Fourier integral operator. Let the symbol Ψ_{WE} = Ψ_{WE}(x, z, p, ξ, ζ, θ) (where θ is the p-covector) be given by the pull back of ψ under the map from T^{∗}ℝ_{(x,z,p)}^{2n−1}\0 to T^{∗}ℝ_{(r,s,t)}^{2n−1}\0 induced by the canonical relation of A_{WE}. The composition Ã_{WE}F is a p-family of pseudodifferential operators with principal symbol Ψ_{WE}(x, z, p, ξ, ζ, 0).

Proof. It is sufficient to show that the map

(34)  δc ↦ ∫ (Ξ(z)H(0, z)^{∗} ∫ H(0, z′)I_{2}I_{1}δc dz′) (x − \frac{h}{2}, x + \frac{h}{2}, ph) dh,

microlocally has principal symbol equal to 1. In this proof we will omit the cutoff functions that are part of the symbols; the calculations will be valid microlocally on the support of a cutoff.

Using an oscillatory integral representation of H similar to the one in the proof of Proposition 6.1, we find that the principal contribution to the kernel of this map, as a function of (x, z, p; x′, z′), can be written as

(2π)^{−(2n−1)} ∫ Ξ(x − \frac{1}{2}h, x + \frac{1}{2}h, ph, −\frac{∂S}{∂s} , −\frac{∂S}{∂r} , −\frac{∂S}{∂t} , z)

× \overline{A(z, x − \frac{1}{2}h, x + \frac{1}{2}h, ph, y_{0}, η_{0J})}A(z′, x′, x′, 0, y_{0}, η_{0J})

× e^{i[−S(z,x− \frac{1}{2} h,x+ \frac{1}{2} h,ph,y_{0I},η_{0J})+S(z′,x′,x′,0,y_{0I},η_{0J})]} dy_{0I}dη_{0J} dh.

We expand the phase in a Taylor series around (x′, z′, h) = (x, z, 0) and identify the gradient at (x, z, 0),

− \frac{∂S}{∂x} (z, x, x, 0, y_{0I}, η_{0J}) = σ(z, x, x, 0, y_{0I}, η_{0J}) + ρ(z, x, x, 0, y_{0I}, η_{0J}),

− \frac{∂S}{∂z} (z, x, x, 0, y_{0I}, η_{0J}) = ζ(z, x, x, 0, y_{0I}, η_{0J}),

− \frac{∂S}{∂h} (z, x, x, 0, y_{0I}, η_{0J})

= −\frac{1}{2}σ(z, x, x, 0, y_{0I}, η_{0J}) + \frac{1}{2}ρ(z, x, x, 0, y_{0I}, η_{0J}) + pτ(z, x, x, 0, y_{0I}, η_{0J}).

Applying a change of variables, (y_{0I}, η_{0J}) ↦ (σ, ρ, ζ), the phase takes the form

⟨σ + ρ, x − x′⟩ + ζ(z − z′) + ⟨\frac{1}{2}(ρ − σ), h⟩.

The amplitude factor Ξ\overline{A}A becomes equal to one by the calculations in the proof of Proposition 6.1. Upon changing integration variables, σ = \frac{1}{2}ξ − θ, ρ = \frac{1}{2}ξ + θ, the oscillatory integral takes the leading-order form (microlocally)

(2π)^{−(2n−1)} ∫ e^{i(⟨ξ,x−x′⟩+ζ(z−z′)+⟨θ,h⟩)} dξdζdθdh.

It follows by integrating the θ, h variables that (34) indeed has principal symbol equal to 1 microlocally. □

8. Annihilators

It was observed in [18] that there are pseudodifferential operators that annihilate the data, due to the fact that the inverse problem is formally overdetermined. On transformed data Ã_{WE}d, these are given by \frac{∂}{∂p_{i}}, i = 1, …, n − 1 hence the annihilators are given by

⟨Ã_{WE}⟩^{−1} \frac{∂}{∂p_{i}} Ã_{WE},

Where ⟨Ã_{WE}⟩^{−1} is a regularized inverse for Ã_{WE}, that is a microlocal inverse for a subset of T^{∗}ℝ_{(x,z,p)}^{2n−1} where Ã_{WE} is invertible.

We define the operator K̃^{∗} that maps d to a function (s, r, z) by

(K̃^{∗}d)(s, r, z) = 2c_{0}((s + r)/2, z)^{3}

× (R_{2}Ξ(z)Q_{−,s}^{∗}(z)^{−1}Q_{−,r}^{∗}(z)^{−1}H(0, z)^{∗}Q_{−,s}(0)^{−1}Q_{−,r}(0)^{−1}D_{t}^{−2}d)(s, r, z).

Using (27), we observe that the operator K̃^{∗} acting on d = F δc yields

K̃^{∗}ψd = Ψ(s, r, z, D_{s}, D_{r}, D_{z})I_{1}δc

(for the definition of Ψ see remark above Proposition 6.1). Applying the operator M given by the multiplication by r − s to this equation yields

MK̃^{∗}ψd = [M, Ψ(s, r, z, D_{s}, D_{r}, D_{z})] I_{1}δc,

i.e. there is only a lower order contribution, since (r − s)δ(r − s) = 0 and hence MI_{1} = 0 (cf. (18)).

For the subset in T^{∗}ℝ_{(s,r,z)}^{2n−1} where the operator Ψ(s, r, z, D_{s}, D_{r}, D_{z}) is microlocally elliptic, the operator

M̃:= M − [M, Ψ(s, r, z, D_{s}, D_{r}, D_{z})]Ψ(s, r, z, D_{s}, D_{r}, D_{z})^{−1}

is an annihilator of K̃^{∗}ψd to all orders.

Corollary 8.1. A pseudodifferential annihilator of the data is given by

W = KM̃K̃^{∗}.

Note that W = W [c_{0}] depends on the background medium. The semi-norm ∥W [c_{0}]d∥ can be viewed as the wave equation analog of the differential semblance functional [19].

References

[1] G. Beylkin. Imaging of discontinuities in the inverse scattering problem by inversion of a causal generalized Radon transform. J. Math. Phys., 26(1):99–108, 1985.

[2] S. Brandsberg-Dahl, M. V. De Hoop, and B. Ursin. Focusing in dip and ava compensation on scatteringangle/azimuth common image gathers. Geophysics, 2002. In print.

[3] Jon F. Claerbout. Imaging the Earth’s Interior. Blackwell Scientific Publications, 1985.

[4] R. W. Clayton. Common midpoint migration. Technical Report SEP-14, Stanford University, 1978.

[5] Maarten V. de Hoop, Jérôme H. Le Rousseau, and Ru-Shan Wu. Generalization of the phase-screen approximation for the scattering of acoustic waves. Wave Motion, 31(1):43–70, 2000.

[6] J. J. Duistermaat. Fourier Integral Operators. Birkhäuser, Boston, 1996.

[7] Victor Guillemin. On some results of Gel′fand in integral geometry. In Pseudodifferential operators and applications (Notre Dame, Ind., 1984), pages 149–155. Amer. Math. Soc., Providence, RI, 1985.

[8] S. Hansen. Solution of a hyperbolic inverse problem by linearization. Comm. Partial Differential Equations, 16:291–309, 1991.

[9] Lars Hörmander. The Analysis of Linear Partial Differential Operators, volume 4. Springer-Verlag, Berlin, 1985.

[10] A. P. E. ten Kroode, D. J. Smit, and A. R. Verdel. A microlocal analysis of migration. Wave Motion, 28:149–172, 1998.

[11] Anders Melin and Johannes Sjöstrand. Fourier integral operators with complex-valued phase functions. In Fourier integral operators and partial differential equations (Colloq. Internat., Univ. Nice, Nice, 1974), pages 120–223. Lecture Notes in Math., Vol. 459. Springer, Berlin, 1975.

[12] C. J. Nolan and W. W. Symes. Global solution of a linearized inverse problem for the wave equation. Comm. Partial Differential Equations, 22(5-6):919–952, 1997.

[13] A. M. Popovici. Prestack migration by split-step dsr. Geophysics, 61(5):1412–1416, 1996.

[14] Marie Prucha, Biondo Biondi, and William Symes. Angle-domain common image gathers by waveequation migration. In 69th Ann. Internat. Mtg, pages 824–827. Soc. Of Expl. Geophys., 1999.

[15] Rakesh. A linearized inverse problem for the wave equation. Comm. Partial Differential Equations, 13(5):573–601, 1988.

[16] Christiaan C. Stolk. Microlocal analysis of a seismic linearized inverse problem. Wave Motion, 32:267– 290, 2000.

[17] Christiaan C. Stolk. Microlocal analysis of the scattering angle transform. Preprint, Dept. of Computational and Applied Mathematics, Rice University, http://www.caam.rice.edu/ cstolk/angle.ps, March 2000.

[18] Christiaan C. Stolk and Maarten V. de Hoop. Modeling and inversion of seismic data in anisotropic elastic media. Preprint 1143, Department of Mathematics, Utrecht University, March 2000.

[19] William W. Symes. A differential semblance algorithm for the inverse problem of reflection seismology. Comput. Math. Appl., 22(4-5):147–178, 1991.

[20] Michael E. Taylor. Reflection of singularities of solutions to systems of differential equations. Comm. Pure Appl. Math., 28:457–478, 1975.

[21] Michael E. Taylor. Pseudodifferential Operators and Nonlinear PDE, volume 100 of Progress in Mathematics. Birkhäuser, Boston, 1991.

[22] François Treves. Introduction to Pseudodifferential and Fourier Integral Operators, volume 2. Plenum Press, New York, 1980.

Christiaan C. Stolk, The Mathematical Sciences Research Institute, 1000 Centennial Drive #5070, Berkeley CA 94720-5070 USA

  E-mail address: cstolkcaam.rice.edu

Maarten V. de Hoop, The Mathematical Sciences Research Institute, 1000 Centennial Drive #5070, Berkeley CA 94720-5070 USA

  E-mail address: mdehoopMines.EDU