ON THE MEASUREMENT OPERATOR FOR SCATTERING IN LAYERED MEDIA

. We describe new mathematical structures associated with the scattering of plane waves in piecewise constant layered media, a basic model for acoustic imaging of laminated structures and in geophysics. Using explicit formulas for the reﬂection Green’s function it is shown that the measurement operator satisﬁes a system of quasilinear PDE with smooth coeﬃcients, and that the sum of the amplitude data has a simple expression in terms of inverse hyperbolic tangent of the reﬂection coeﬃcients. In addition we derive a simple geometric description of the measured data, which, in the generic case, yields a natural factorization of the inverse problem.


1.
Introduction. There is an extensive literature concerning acoustic scattering in layered media, motivated by applications to imaging technologies such as seismic and ultrasound. Nevertheless, one of the simplest models, that of piecewise constant structure where the impedance profile is a step function, has a richer mathematical structure than appears to have been previously recognized. The purpose of the present paper is to describe this structure briefly, and to mention some possible implications, both for solving the inverse problem and for wider theory. The key insight is provided by a relatively recent combinatorial analysis of scattering sequences.
Traditional approaches to piecewise constant layered media have centred predominantly on: discretization or assumptions of constant layer thickness [3] [4]; approximation, involving linearization or other forms of multiple suppression [1, Chapter 2] [10]; probabilistic analysis, aimed at extracting bulk properties rather than exact structure [11] [5] [6]; or abstract scattering theory [16]. We depart from the traditional approach and use the formulas in [7] to analyze the exact structure of the full nonlinear measurement operator from the deterministic perspective, uncovering several unexpected features. Firstly, the measurement operator itself satisfies a system of quasilinear PDE in which physical parameters occur as variables. Unlike the underlying wave equation which has discontinuous (piecewise constant) coefficient, the equations governing the measurement operator have smooth coefficients. This is noteworthy from the general perspective of PDE with discontinuous coefficients. A second unexpected feature of the measurement operator is a beautiful trace formula whereby the sum of the scattering amplitudes converges to a function of the underlying reflectivities that can be expressed very simply in terms of hyperbolic tangent. Thirdly, we present a simple geometric description of the scattering experiment that separates out the effects of the two types of physical parameters at play, reflectivities and travel time thicknesses. The scattering experiment is shown to consist of translation of a fixed function, given explicitly a priori, followed by a pushforward and truncation. The translation vector is simply the list of reflectivities, and the pushforward stems directly from the vector of travel time thicknesses, with truncation corresponding to time limited measurement. These three results give new qualitative and analytic insights into the inverse problem for layered media, and suggest possible features of scattering in other, more complicated settings.
The paper is organized as follows. We detail the technical setting in §1.1 and distinguish two classes of layered media according to the associated impedance profile: piecewise constant media, whose impedance profile is a step function; and continuous layered media whose impedance profile has a square integrable logarithmic derivative. The point of this distinction is that the latter class has a well-developed deterministic theory [15] [14] [12], whereas the former class, which is the subject of this paper, does not. We then state the main results in §1.2.
The main results are proved in §2, and the paper concludes with a brief discussion in §3.
1.1. Scattering in layered media. Consider layered media whose structure varies in the depth direction z only, orienting coordinates such that z increases downward into the medium and decreases upward. Given initial conditions that depend only on depth, the evolution of particle velocity u(z, t) in a layered medium is governed by the one dimensional wave equation where ζ(x) > 0 denotes acoustic impedance, and x is travel time depth Suppose further that ζ(x) = C 0 is constant for x < 0, and that lim x→∞ ζ(x) = 1. A basic problem in acoustic imaging is to determine ζ(x) for x > 0 from scattering data. That is, suppose a source waveform W ∈ C ∞ c (R) is supported on the negative half line, and initially travels downward, corresponding to initial conditions The resulting echoes d(t) = u(0, t) recorded at x = 0 during a certain time interval 0 < t < T comprise scattering data that carries information about ζ(x) for 0 < x < T 2 . The (forward) scattering problem is to determine d given W and ζ; the inverse scattering problem is to determine ζ given d, or given d and W . Two particular classes of impedance functions have been the focus of considerable research: (i) piecewise constant ζ, which, apart from the case of evenly spaced jump points, have been studied mainly from the probabilistic standpoint; and (ii) impedance functions whose logarithmic derivative is square integrable. Note that the logarithmic derivative of ζ arises naturally as the only non-constant coefficient when (1) is rewritten in the form The condition ζ ζ ∈ L 2 (R) excludes the possibility of ζ having a jump discontinuity (as does the closely related condition log ζ ∈ H 1 (R + )), so the aforementioned two classes are disjoint. At first blush one might guess that the latter class is more general and therefore more difficult to analyze. But, curiously, historical evidence suggests the opposite. The scattering problem for impedance functions whose logarithmic derivative is square integrable has a rich theory, including a remarkable nonlinear Plancherel theorem, whereas there appears to be no comparable analysis established for the case where ζ is a step function.
In the present paper we focus on the piecewise constant class of layered media, analyzing the nonlinear measurement operator ζ → d restricted to impedance functions of the form Here H(x) denotes the Heaviside function, and jump points are indexed according to their natural order The number of layers n ≥ 2 is assumed to be fixed; but the particular value of n can be arbitrarily large. Choosing x * j ∈ (x j−1 , x j ) arbitrarily for 1 ≤ j ≤ n + 1, set one can easily reconstruct ζ in (4) from the pair (τ, r), giving a correspondence We may thus regard the reflectivity vector r and the travel time per layer vector τ as physical parameters characterizing a given piecewise constant layered medium. It turns out that the measurement operator ζ → d for the scattering problem can be expressed very simply in terms of (τ, r); for convenience we therefore work with (τ, r) instead of ζ.
For piecewise constant layered media the measured scattering data has the form where G (τ,r) denotes the impulse response, or Green's function, at the boundary. The latter is a delta train, expressible in the form where, for each i, 0 < t i < t i+1 and a i = 0. We refer to the points t i as arrival times (always indexed according to their natural order), and the corresponding a i as amplitudes.
The nature of the source waveform W occurring in (9), and whether or not it is known, are contingent on instrumentation and the experimental context at hand, be it ultrasonic imaging of laminated structures or seismic imaging of geological strata. The problem of removing W from the measured data when W is given (deconvolution) or unknown (blind deconvolution) is of course essential in any practical context, but it is contingent on the context, and we shall not treat it here. The present paper concentrates on the impulsive limit W → δ, in which the measured data is a time limited truncation of the Green's function G (τ,r) itself. The results obtained nevertheless apply to the more realistic case where W is nontrivial. It is evident from (11) that d is represented by the sequence of pairs (a i , t i ) for 1 ≤ i ≤ i max ; in the sequel we shall make use of this representation whenever it is convenient, without further comment. Henceforth we refer to the mapping (τ, r) → d as the measurement operator (so, in effect, the measurement operator is just the measured data viewed as a function of the physical parameters).

1.2.
Main results. We have three main results. The first is that the measurement operator (τ, r) → d is itself governed by a system of quasilinear second order equations in which the physical parameters (τ, r) occur as independent variables. This system of equations essentially comprises a model of the acoustic reflection experiment in the context of piecewise constant layered media. Call a travel time vector τ , or a pair (τ, r), generic if the τ j (1 ≤ j ≤ n) are linearly independent over the integers. Note that the set of nongeneric τ ∈ R n + is meagre.
Theorem 1.1. Referring to (11), on a generic open subset of parameter space (τ, r) ∈ R n + × (−1, 1] n , for every 1 ≤ i ≤ i max and 1 ≤ j ≤ n, the pair a i , t i satisfies the system of quasilinear equations = 0, a i being bounded in a neighbourhood of r = 0.
Note that the expression ∂ ∂τn+1 denotes the zero operator in the above statement. The second result is a trace formula that expresses the sum of the amplitude data very simply in terms of hyperbolic tangent.
The third main result was first presented in [8]-but we include it here for completeness. It represents the measurement operator in simple geometric terms as a composition of operations which depend on the physical parameters (τ, r) in a transparent way. The operations are: translation by the vector 1 2 r (of an explicitly given function ψ : R n → R); sampling on the integer lattice Z n ; pushing forward onto R via τ * (the natural pushforward associated with τ ); and multiplication by χ [0,T ] . More formally, given w ∈ R n let T w denote the translation operator Sampling a function f ∈ C(R n ) on the integer lattice Z n may be represented by multiplication with the tempered distribution the sampling operator being the distribution valued mapping The pushforward associated with a given v ∈ R n maps a subclass of distributions on R n to distributions on R as follows. Let S(R n ) denote Schwartz space and F ∈ S (R n ) a tempered distribution. Where possible, define v * F ∈ S (R) by the formula This requires that F (ρ) be well defined for arbitrary ϕ ∈ S(R), which limits the domain of v * to a proper subset of S (R n ); the subset evidently includes distributions of the form f S when f ∈ C(R n ) is bounded. In more concrete terms, temporarily add a superscript to the delta function to indicate the dimension of the underlying space, writing δ n ∈ S (R n ); the pushforward v * is linear and maps translates of δ n to translates of δ 1 according to the formula

Theorem 1.3 ([8]
). There exists a locally polynomial function ψ : R n → R, defined a priori, independently of (τ, r) ∈ R n + × (−1, 1] n , such that (20) d = χ [0,T ] τ * ((T 1 2 r ψ)S). The explicit form of ψ is given in the next section. Theorem 1.3 factors the measurement operator into geometrically simple operations that separate the respective effects of r and τ . In particular, we define the covering data, which depends only on r, to be the distribution  Let Z + denote the non-negative integers, including 0; when referring to entries of an n-tuple k, we adopt the convention that k n+1 = 0.  (10) is given by the formula Note that each y ∈ R n can be uniquely expressed in the form y = k + r/2, where k ∈ Z n and r ∈ (−1, 1] n . Based on this, define ψ : R n → R in terms of the Green's function coefficients in (23) by the formula (24) Then, by construction, Observe that the order relation (26) is constant in a neighbourhood of τ . Moreover, the set of k ∈ Z n + such that τ, k < T is constant in a neighbourhood of any generic τ ∈ R n + , as is the ordering of the corresponding (finite) set of numbers τ, k . Letting k (i) ∈ Z n + denote the lattice point such that t i = τ, k (i) , it follows that the sequence is constant in a sufficiently small neighbourhood of any generic τ . Theorem 2.1 implies that in the generic case each amplitude and arrival time pair (a i , t i ) has the form for some k ∈ {1} × Z n−1 + , and that the value of k is constant in a sufficiently small neighbourhood of τ by the foregoing argument. Thus in order to prove Theorem 1.1 it suffices to verify that the pair (28) satisfies equation (12) and boundary condition (13). Given that t i = τ, k , where k ∈ {1} × Z n−1 + is constant in a neighbourhood of τ , it follows that (29) ∂t i ∂τ j = k j and ∂t i ∂τ j+1 = k j+1 , whereby equation (12) reduces to a linear equation in a i , For each fixed j in the range 1 ≤ j ≤ n, the representation (28) implies that, viewed as a function of r j , a i is proportional to f (kj ,kj+1) (r j ). Thus a i satisfies the equation (30) provided f (kj ,kj+1) (r j ) does. Relabelling the pair (k j , k j+1 ) as (p, q) and r j as x, we have shown that a i satisfies (12) provided f (p,q) (x) satisfies the ODE Straightforward (but tedious) evaluation of the formula (22), plugged into the lefthand side (31), verifies that f (p,q) satisfies (31), including in the degenerate cases where min{p, q} < 1.
Given the representation (28), the boundary condition (13) reduces to the equation pq f (p,q) (1) = 0, which can be read directly from the formula (22). This completes the proof of Theorem 1.1.
Somewhat more is true. Standard solution methods applied to equation (31) show that (31) together with boundedness of g near 0 and the boundary condition pq g(1) = 0 determine the functions f (p,q) up to a scalar multiple. It follows that the quasilinear system (12) and boundary condition (13) determine the functions a i up to a scalar multiple-given local linearity of the t i in τ .
Given a pair of n-tuples (w, z) ∈ D n × T n , let Ψ wj zj : D → D denote the linear fractional transformations which if w ∈ D n are automorphisms of the Poincaré disk. (Note that if w j ∈ T is on the disk boundary then Ψ wj zj (v) = z j w j collapses to a constant function.) The following well-known result (cf. [2, §2.5]) is sometimes referred to as a backward recurrence formula.  G (τ,r) (σ) = Ψ r1 e iτ 1 σ • Ψ r2 e iτ 2 σ • · · · • Ψ rn e iτn σ (0). Taking the Fourier transform of (23) with respect to t yields the key identity behind Theorem 1.2.
Proof. Comparing Theorem 2.2 to the Fourier transform with respect to t of Theorem 2.1 yields that for every (τ, r) ∈ R n Since the mapping (τ, σ) → (e iτ1σ , . . . , e iτnσ ) is a surjection from R n + × R onto T n , (35) follows immediately from (36).
Let Ψ(r, z) denote the composition of disk automorphisms on the left-hand side of (35), The next step is to extend the domain of Ψ(r, z), as defined by the formula (37), to z ∈ C n , and to recognize (35) as the Taylor series expansion of a holomorphic function.
There exists a sequence of positive numbers δ = (δ 1 , . . . , δ n ) such that for every r ∈ [−ε, ε] n the function Ψ(r, z) is holomorphic on the polydisk Proof. By Hartogs' Theorem it suffices to prove that for each z j taken separately, there exists a positive δ such that Ψ(r, z) is holomorphic in z j for all r ∈ [−ε, ε] n and all z ∈ P δ . Using the fact that each disk automorphism Ψ rν zν is holomorphic on a disk centred at the origin of radius larger than 1, independently of r ∈ [−ε, ε] n , one can prove the the desired univariate statement by a straightforward inductive argument.
It follows from Lemma 2.3 that the domain of convergence of the Taylor series of Ψ(r, z) (expanded about z = 0) includes P δ . By Corollary 1 and the Fourier transform of Theorem 2.1, the Taylor coefficients are precisely the numbers (39) a i = n j=1 f (kj ,kj+1) (r j ).
Since 1 = (1, . . . , 1) belongs to the interior of P δ we may evaluate the Taylor series at z = 1 to yield and are guaranteed that the series converges absolutely, uniformly in r ∈ [−ε, ε] n . Theorem 1.2 then follows from the observation that 3. Discussion. Theorem 1.1 shows that in the special case of piecewise constant layered media, the measurement operator is itself the solution to a PDE that has a different character from the wave equation governing sound propagation through physical media. The PDE (12) models the reflection experiment as a function of medium parameters, as opposed to (1), which models a physical process evolving in time and space. Referring to (12) as a measurement model, there are several points to be made. It is a remarkable theoretical fact that a measurement model exists in the first place. But more than this, the measurement model has smooth coefficients, making it more amenable to direct analysis than the original wave equations, whose discontinuous coefficients place it outside the scope of much standard theory. From a broader theoretical perspective Theorem 1.1 suggests two things. Firstly, is the existence of a measurement model a general phenomenon? Secondly, the measurement model gives qualitative information about the measurement operator beyond the observation of nonlinearity (nonlinearity being a typical feature of inverse problems for linear PDE, where one seeks to determine unknown coefficients from measurement of a restricted solution). But how should the precise nature of the nonlinearity be characterized? If the existence of a measurement model turns out to be a more general phenomenon, then the classification of the measurement model as a PDE provides a natural way to categorize the measurement operator.
As mentioned at the end of the proof of Theorem 1.1, the equation (12) and boundary condition (13) determine the the functions a i up to a scalar multiple, given the linear structure t i = τ, k of the arrival times. It is natural to ask whether the trace formula (14) then determines the a i exactly, given their value up to a scalar. We have not addressed this delicate question because, apart from considerable technical difficulties, it can be answered much more cleanly by way of complexification. That is, one can show that the Fourier transform of a complexified version of the covering data (21) is determined uniquely by a system of linear equations having smooth coefficients and a trace condition that complexifies (14). Details of the complexification are rather involved and comprise a separate paper [9], which is complementary to the present work.
The inverse problem for scattering in layered media factors naturally according to the statement of Theorem 1.3 into two steps. The first step is to determine the pushforward τ * . That is, given the sequence of arrival times t 1 < · · · < t imax , determine τ such that t i = τ, k is the projection of the appropriate part of the integer lattice {1} × Z n−1 . Provided the cutoff time T > 0 is sufficiently large and τ is generic, there is a straightforward algorithm by which this first step can be solved uniquely (see [8] for details). The solution to the first step determines the lattice point k (i) associated to each amplitude a i = T 1 2 r ψ(k (i) ), for 1 ≤ i ≤ i max , and hence the restriction of the covering data D to the set of lattice points k such that 0 < τ, k < T . The second step of the inverse problem is to determine r given this restriction of the covering data. Here the simple geometric dependence of the covering data on r comes into play. The fact that the covering data is obtained by sampling the translate of ψ by 1 2 r on the integer lattice means that the second step of the inverse problem can be viewed as a global optimization problem involving all the amplitudes a i on essentially equal footing. There is no need to discard information by suppressing multiple reflections, as is traditionally done in seismic signal processing. On the contrary each a i is the value of a translate of ψ at a known lattice point k (i) and so constrains the possible translate 1 2 r. One can carry this further by analyzing the actual formula for ψ, but we shall not do so in the present paper. The main point is that there is a simple geometric picture implicit in Theorem 1.3 whereby the inverse scattering problem generically factors into two steps, an inverse projection problem to determine τ , followed by a global optimization problem to determine r, which is the translate of ψ whose evaluation at the appropriate lattice points best approximates the measured amplitude data.
In non-generic cases such as where τ is proportional to 1, the sequence of arrival times t 1 < · · · < t imax need not uniquely determine τ , and the inverse problem, though still tractable, fails to decouple.
In practice the energy of the source waveform (3) is sometimes not known, and so, in place of d in the form (11), one has an unknown scalar multiple of thisuncalibrated data, as described in [13]. We remark that the optimization problem of determining r given the restriction of the covering data D to the set of lattice points k such that 0 < τ, k < T can still be solved in this uncalibrated setting. And if the series in Theorem 1.2 converges sufficiently rapidly (equivalently, if the cutoff time T > 0 is sufficiently large), the partial sum To summarize, we have used the explicit formula of Theorem 2.1 to show that the measurement operator for scattering in piecewise constant layered media satisfies a system of quasilinear equations and an elegant trace formula, and that the inverse problem has an simple geometric description that separates out the effects of layer thickness τ and reflectivity r.