A time-domain probe method for three- dimensional rough surface reconstructions

The task of this paper is to develop a Time-Domain Probe Method for the 
reconstruction of impenetrable scatterers. The basic idea of the method is 
to use pulses in the time domain and the time-dependent response of the scatterer 
to reconstruct its location and shape. The method is based on the 
basic causality principle of time-dependent scattering. The method is independent 
of the boundary condition and is applicable for limited aperture scattering data. 
 
   In particular, we discuss the reconstruction of the shape of a rough surface in 
three dimensions from time-domain measurements of the scattered field. 
In practise, measurement data is collected where the incident field is given by 
a pulse. We formulate the time-domain field reconstruction problem 
equivalently via frequency-domain 
integral equations or via a retarded boundary 
integral equation based on results of Bamberger, Ha-Duong, Lubich. In contrast 
to pure frequency domain methods here we use a time-domain characterization 
of the unknown shape for its reconstruction. 
 
   Our paper will describe the Time-Domain Probe Method and relate it to previous 
frequency-domain approaches on sampling and probe methods 
by Colton, Kirsch, Ikehata, Potthast, Luke, Sylvester et al. The approach 
significantly extends recent work of 
Chandler-Wilde and Lines (2005) and Luke and Potthast (2006) on the time-domain 
point source method. 
We provide a complete convergence analysis for the method for the rough surface 
scattering case and provide numerical simulations and examples.


Introduction
We consider the scattering of a time-dependent pulse from a rough surface in three dimensions. It is usually easy to send waves into some area of the earth or the human body and measure their scattered field. The task of inverse scattering theory is to explore or visualize the structure of the unknown region and to evaluate its properties. This includes localization and shape reconstruction of the objects and interfaces and the determination of material parameters like density, conductivity or permittivity.
A rough surface denotes a surface which is usually a non-local perturbation of an infinite flat surface such that the surface lies within a finite distance of the original plane. In particular, we expect the rough scattering surface Γ to be the graph of some bounded continuous function f : R 2 → R, i.e. (1) Γ := x = (x 1 , x 2 , x 3 ) ∈ R 3 : x 3 = f (x 1 , x 2 ) . All these points will then be on the unknown surface Γ. Reconstructions of U s (x, t) can be carried out by the time-domain point source method [3] or the potential method [1].
Our goal is the reconstruction of the shape of the surface by given measurements in time.
To this end, we describe the problem settings in the time-domain and its relation to the frequency-domain settings.
Frequency-domain problems are widely studied, for the rough surface case we refer to [4], [5], [6], [8], [9], [10], [11], [12]. There exist many methods to solve the inverse problem, most of which have been worked out for the bounded obstacle case. For example, iterative methods update some reconstruction using gradients or the Fréchet derivative with respect to the unknown boundaries [7], [26], [13]. The Point Source Method [3] and the Kirsch-Kress Method [13] reconstruct the full field and then use the boundary condition to find the unknown shape. Probe Methods as introduced by Ikehata, [20], Potthast, Nakamura and others, compare [30], usually define some indicator function via particular incident fields which can be used to construct or visualize the unknown shapes or objects. As examples for further approaches we name Sampling Methods, [14], [31], Range Tests, [21], [29], and Factorization Methods, [22], a survey is given in [30] and in [27].
The inverse rough surface problem in two dimension in the frequency domain is for example discussed by DeSanto and Wombell, [16], [17], and for the periodic case by Elschner and Yamamoto, [18]. The three-dimensional rough surface case for a single frequency is worked out in [1]. For the numerical realization of the rough surface case we employ the multi-section approach presented by Heinemeyer, Lindner, Potthast, [19].
In acoustic applications time-domain measurements are usually relatively easy to obtain whereas pure frequency domain methods do not take the full advantage of the data which are available. Here we suggest a method taking full time-measurements into account.
It is a significant step forward compared to earlier work by Chandler-Wilde and Lines [3] and by Luke and Potthast [25]. The Time-Domain Probe Method is a full time-domain scheme: it is based on causality which cannot be directly used in the frequency domain. Thus, it naturally incorporates the knowledge of scattered fields for many frequencies.
The method can be seen as an extension of the probing methods of Ikehata, Potthast, Nakamura, Sini and others. Our approach here will not rely on a particular probe in the time domain, but will in principle work with a large variety of incident time-dependent fields. In particular, the probing fields here do not need to have a singularity of any type. This avoids numerical instabilities which is one of the key problems for the realization of the time-domain probe schemes. In contrast to engineering schemes from travel-time tomography here we use a full reconstruction of the time-dependent scattered field U s which is exploited for reconstructing the unknown surfaces.
As a part of the probing procedure we employ the time-domain field reconstruction problem by frequency-domain inverse methods and Fast Fourier Transform (FFT) via a single layer potential approach as first used by Kirsch and Kress in 1986, [13]. An alternative has been developed with the Point Source method by Chandler-Wilde and Lines, [3]. Here, we need to employ the Dirichlet Green's function for the half-space instead of the standard free-space fundamental solution for the Helmholtz equation as it has been carried out for the forward problem, [4], [5].
For the time-domain probe method the actual boundary condition is not explicitly used in the reconstruction algorithm. All arguments will be analogous for other type of boundary conditions. Thus, we expect the method to be independent of the particular physical nature of the object under consideration. Here, we will only investigate the case of the Dirichlet boundary condition in detail and leave other boundary conditions to future research.
Our presentation plan is as follows. In Section 2 we study the time-domain scattering problem and describe the retarded potential boundary integral equation which is often used to model and simulate the scattering process. Section 3 serves to summarize results about the frequency-domain scattering problem from [4], [5], [6], [8], [9], [10]. In Section 4, we present the Time-Domain Probe Method and in Section 5 we analyse its convergence. Section 6 shows numerical results for the rough surface problem.
For 0 < β ≤ 1, let BC 1,β (R 2 ) denote the set of those bounded continuous functions f : R 2 → R whose gradients are bounded and uniformly Hölder continuous with index β > 0. We require that the scattering surface Γ is the graph of a function f ∈ BC 1,β (R 2 ) which is bounded by two constants f − , f + > 0 such that and that our surfaces Γ can be represented as The domain of propagation is given by To indicate the dependence on f we also use the notation Ω f and note that Γ = ∂Ω.

The Time-Domain Problem
Let σ > 0 a constant controlling the temporal distribution of a time-dependent pulse. We study the scattering of a pulse given by on a rough surface with the free-space fundamental solution where we assume that g is chosen such that U i (x, t) is compactly supported and C 3 -smooth with respect to time. Here, we consider the Fourier transform F with respect to frequency or time, respectively, using capital letters for the time-dependent fields and small letters in the frequency domain, i.e. we define Throughout this paper we set the wave speed equal to one. We formulate the time-domain problem as follows.
Problem 1 (The Direct Problem in Time-Domain). Given an incident pulse U i (x, z, t) which is C 3 -smooth and compactly supported with respect to time, find a solution U s ∈ We assume that the Fourier transform u s (x, κ) of U s (x, t) uniformly satisfies the boundedness condition (compare [4]) with some constant c for any fixed κ, (κ) ≥ 0, and we demand u s to satisfy a limiting absorption principle (see (19)).
Taking the Fourier transform with respect to time on both sides of the wave equation (7) we obtain the Helmholtz equation (10) ∆u s + κ 2 u s = 0. and we observe Under the above conditions we obtain uniqueness and solvability of the time-domain problem by combining uniqueness and solvability of the frequency problems and the bounds on the inverse as worked out in [5] by an application of the Fourier transform.
For a real-valued mapping U : Ω×R → R it holds that (FU (x, ·))(κ) = (FU (x, ·))(−κ). Further, in three-dimensions (in contrast to the two-dimensional case) there are no difficulties arising at κ = 0 as the fundamental solution Φ(x, y, κ) has no singularity in κ = 0 and in the following we can assume that κ ∈ R.
Often, time-domain solutions to scattering problems are represented by time-domain integral equations [2,24,15]. Let Γ * be some surface with Γ * ⊂ Ω. By an application of the Fourier Transform to a single-layer potential representation in the frequency domain, every solution of Problem 1 which allows such a representation for allmost all wave numbers κ with a density which is in L 2 (R) with respect to the wave number can be represented as a potential For the limit x approaching the boundary we derive that any density ϕ of the above representation must satisfy the boundary integral equation where we define is known as retarded potential. For a investigation of numerical methods to solve the retarded boundary integral equation for bounded surfaces, i.e. Γ * = ∂Ω where Ω is a bounded subset of R 3 , we refer for example to [15]. A complete theoretical background using Laplace transforms can be found in [2] and in [24].
Here, for simulation and inversion we employ calculations using FFT and frequency domain methods. An application of the FFT with respect to time to equation (13) leads to the standard single-layer boundary integral equation in the frequency domain for all frequencies κ ∈ R. We solve this problem for κ in a uniform grid of points and employ the inverse FFT to obtain an approximate solution to the retarded integral equation. For more details about equivalence and estimates we refer to the arguments worked out in [25]. We also refer to [28] where this equivalence has been used to construct a time-domain filter for field reconstruction from a family of frequency-domain filters. In the next section we summarize some relevant results on the frequency problem for the rough surface setting.

The Direct Scattering Problem in the Frequency Domain
In the frequency domain problem we consider the scattering of an acoustic field from the rough surface Γ. The incident field is due to a point source in z ∈ R 3 defined by u i (x, κ) = Φ(x, z, κ), where Φ is the fundamental solution of the Helmholtz equation (5). As we have seen we can limit our attention to the case κ ≥ 0. The direct problem is to find the scattered field u s (·, κ) ∈ C 2 (Ω)∩C(Ω) such that the total field u(·, κ) = u i (·, κ)+u s (·, κ) is a solution of the Helmholtz equation (16) ∆u(x, κ) + κ 2 u(x, κ) = 0 for x ∈ Ω.
Furthermore, the total field is required to satisfy the Dirichlet boundary condition (17) u(x, κ) = 0 for all x ∈ Γ, and the scattered field is supposed to be bounded in space, i.e.
for some constant c > 0 only depending on k. In the case that the wave number is a positive real number, we follow [4] and require the limiting absorbing principle, i.e. that for sufficiently small > 0 the solution with wave number κ = k 0 + i exists and the pointwise limit is satisfied for every x ∈ Ω.
Problem 2 (Direct Rough Surface Scattering Problem). Let u i (·, κ) be an incident field due to a point source at the point z ∈ Ω, i.e.
We can convert this scattering problem into a boundary value problem seeking the scattered field in the form Then, the Dirichlet boundary condition of the direct scattering problem with the incident field (20) yields and thus, the remainder v satisfies the boundary condition We remark that the total field u(·, κ) satisfies the direct scattering problem if and only if v(·, k) solves the following boundary value problem, see [4]. Furthermore, the direct scattering problem can be reformulated in a well-posed integral equation, see [4], [5].
The solvability of the frequency domain problem via a combined single-and doublelayer approach for mildly rough surfaces has been shown in [4]. The authors first use Fourier arguments on a flat surface Γ 0 . It is shown that the application of a single-layer potential on a flat surface with kernel corresponds to multiplication with its Fourier transform Note that this is nonzero for h sufficiently small on k ∈ R 2 and that the function at k = κ is bounded. The operator is boundedly invertible on the set of functions whose Fourier transform multiplied with |k| is square integrable, which corresponds to the solvability of the single-layer potential equation For the reconstruction of U s or u s , respectively, on some test surface Γ * (which we choose as a compactly perturbed plane as shown in Figure 2) we need to study the solvability of the single-layer boundary equation. This is carried out in the following result. Theorem 1. Consider the single-layer potential S with kernel G h on a rough surface Γ * . Then for h sufficiently small the operator is boundedly invertible as operator L 2 (Γ * ) → H 1 (Γ * ).
Proof. The theorem is a direct consequence of the results of Lemma 3.3 and Theorem 3.4 in [5], where the density of the single-layer potential on a rough surface is explicitly estimated.

A Time-Dependent Probe Method
Here, we first describe the field reconstruction problem in frequency domain, i.e. for measurements of the total field u(·, κ) = u i (·, κ) + u s (·, κ) for a single fixed wavenumber κ > 0 on a finite measurements plane Γ h,A with for a constant A > 0 our task is to reconstruct u s (x, κ) in Ω. Note that we use the a-priori information h > f + to assure that there are no intersections between Γ h and the unknown surface Γ.
Problem 4 (The Field Reconstruction Problem, frequency domain). Suppose we know the incident field u i (·, κ) = Φ(·, z, κ) and the scattered field u s (·, κ) on the plane Γ h,A for a fixed wavenumber κ with (κ) ≥ 0. We assume that the total field u satisfies the Dirichlet boundary condition u = 0 on the unknown surface Γ. Then, we try to find the total field and the surface Γ such that u is the solution of the direct problem (2) and (u − u i )(·, κ) coincides with u s (x, κ) for all x ∈ Γ h,A .
Let Λ f be a surface given via (1). We define the projection operator for B > 0 by where x = (x, f (x)). In particular, the restrictions of the measurements u s on the plane Γ h,A can be seen as the image of the projection P A u s ∈ P A (L 2 (Γ h )) of the scattered field u s . 8 For the solution of the inverse problem we proceed as follows. We make the ansatz (30) u s (·, κ) = v(·, κ) − Φ(·, z , κ) and seek the projection P A v of the remainder v in a single layer potential, i.e.

Evaluate
Numerically we will employ a finite section version of the above approach, i.e. we study a single-layer potential defined on some finite part Γ * C of Γ * . In the next section we will show convergence of this finite section method to the above solution for the infinite surface Γ * .
We are now prepared to introduce the Time-Domain Probe Method, which relies on time-measurements of the scattered field U s (·, t) on the measurement plane Γ h,A . To this end we assume to know the measurements of the scattered field in time, that is Definition 2. For a point x ∈ Ω we define the first hitting time with respect to the incident field U i by which is called a vertical needle. Then the first hitting parameter λ * is defined as We use a grid of size h λ > 0 for the discretization of the vertical probing, i.e. we employ (39) λ ξ := h λ · ξ, ξ = 0, 1, 2, ...

Algorithm 2 (Time-Domain Probe Method).
Let U s (x, t) be given for all x ∈ Γ h,A and t ∈ R. To identify points x ∈ Γ or to calculate a reconstruction f rec for the surface height function f on Q, respectively, we choose a constant h λ and = 2h λ (where we assume wave speed c = 1) to carry out the following steps: 1. For every x ∈ Q: 2. we successively investigate λ = λ ξ for ξ = 0, 1, 2, .. given by (39): 3. with x λ ξ defined by (37) we reconstruct U s (x λ ξ , t) in the small interval t ∈ (T, T + ) after the first hitting time T = T (x λ ξ ) by the methods described in Algorithm 1.
To achieve a numerical algorithm we will employ a grid of points in the rectangle Q = (a 1 , b 1 ) × (a 2 , b 2 ). With n 1 or n 2 points in x 1 or x 2 direction, respectively, we use the notation We denote our horizontal grid by Q n 1 ,n 2 . For every point x ∈ Q n 1 ,n 2 we investigate the points x λ ξ for ξ = 1, 2, 3, ... until µ(λ ξ ) = 1. For the discrete version we need to make sure that we identify points which are close to the unknown boundary Γ. To this end we need to investigate appropriately chosen intervals (T, T + ) depending on our discretization for simulations. Here, we have used a fixed = 0.2 which was chosen by trial and error.
We will show convergence of the above algorithm in the next section.

Convergence of the Time-Domain Probe Method
The goal of this chapter to prove convergence of the Time-Domain Probe Method for the reconstruction of impenetrable rough surfaces.
As before we write U s α instead of U s to indicate the dependence of the reconstructed scattered field U s in the time-domain on the regularisation parameter α > 0. From now on U s is the true scattered field in the time-domain. We start with the investigation of convergence of U s α to U s . Figure 2: We show the setting for the reconstruction of the scattered field U s (x, t). The surface Γ * is supposed to be in the domain Ω, i.e. above the unknown surface Γ. The measurement surface is Γ A . The figure shows a point x in which U s (x, t) is reconstructed close to or on the boundary of the test surface Γ * .
Theorem 3. Consider the setting shown in Figure 2. Let U i be an incident pulse and U s C,α the reconstructed field using the single-layer approach applied with the truncated surface If the test surface Γ * is in Ω there exists a function C 0 (α) with such that for all C(α) ≥ C 0 (α) we obtain the convergence for every x above or on Γ * .
Proof. For C = ∞ we can choose h > 0 such that the single-layer equation Sϕ = u s from Γ * onto Γ h,A has a unique solution in L 2 (Γ * ). In this case we obtain convergence of the reconstruction for a fixed frequency κ from the standard properties of the Tikhonov regularisation [23]. Denote the Tikhonov operator for Γ * by R α , i.e. R α := (αI + S * P A S) −1 S * P A , and the Tikhonov operator arising from S C := SP C based on Γ * C by R C α . In particular R C α is given by R C α := (αI + P C S * P A SP C ) −1 P C S * P A . Let ϕ α be the regularized density for Γ * , for which we know convergence ϕ α → ϕ towards the true solution ϕ on Γ * for α → 0. Then, we have pointwise convergence This now shows that we have ϕ C(α) α → ϕ when C(α) is chosen appropriately. Finally, we remark that we obtain reconstructions with error estimates uniformly for compact intervals of frequencies with a norm of the inverse operator bounded by Cκ with some constant C, compare [5]. Since the fields are assumed to be in C 3 (R) with respect to time they decay proportional to κ −3 in frequency on the boundary and proportional to κ −2 in the domain Ω. Thus, applying the inverse Fourier transform we obtain convergence for the time-dependent fields as well.
We base our arguments on the following basic properties of the wave equation [32].  Proof. Consider the setting as visualized in Figure 3. The statement is a simple consequence of the fact that for a spherical pulse the total field U (x, t) is zero in space-time range of influence of (x, T (x)). Thus, for t < T (x) the scattered field U s (x, t) can never by positive in x.
Theorem 6 (Convergence of Time-Domain Probe Method). The Time-Domain Probe Method as described in Algorithm 2 provides a complete reconstruction of the surface Γ above the rectangle Q in the sense that for h λ → 0 we have convergence Proof. Assume that x λ ∈ Ω above Q. Then according to Lemma 5 the scattered field U s (x λ , t) is zero for t < T (x) and since U i (x, t) is zero for t < T (x) the same is true for the total field. Further, we want to show that U s (x, t) is zero even in a small neighbourhood (T (x) − , T (x) + )) of T (x). Consider Figure 1, where the evolution of the field U s (x, t) is visualized. For any point x ∈ ∂B R with x ∈ Ω the influence arising from the incident field U i needs some time t > T (x) to reach T (x). This is due to the triangle inequality which states that (52) |z 0 −x| + |x − x| > |z 0 − x| for any pointx ∈ Γ and travel times are calculated by division of the distance by the wave speed c (where T (x) = |z 0 − x|/c). This proves that µ(λ) = 0 for x λ ∈ Ω. Now, consider a point x λ ∈ Γ. Then we know that the scattered field satisfies U s (x, t) = −U i (x, t) and by definition of the first hitting time T (x) we know that in some interval (T (x), T (x) + ) we have |U i (x, t)| > 0. This proves that µ(λ) = 1 for this case. Finally, from both cases and the setup of the needle search which is starting with λ = 0 we obtain convergence of the Time-Domain Probe Method.

Numerical Study of the Time-Domain Probe Method
Our goal here is to show examples for the numerical realization for the Probe Method. We discretize the frequency interval K = [−κ max , κ max ] using a step size h κ = κ max /N with 2N frequency points. Our frequency grid it thus given by with time-stepsize h t . The location of the point source for our simulation is (−3, 0, 10) or (0, 0, 13), respectively, and we choose σ = 15 or σ = 4, respectively, for the Gaussian pulse density. We use forward code developed by Heinemeyer, Lindner and Potthast [19] to calculate the scattered field in the time-domain U s (x, t n ), n = −N, ..., N − 1, for a selection of grid points x in the measurement patch Γ h,A at height h = 10, where we used A = 5. For our numerical implementation we have realized a simplified version of the field reconstruction by a fixed choice of the auxiliary surface Γ * below the unknown surface for the reconstruction of the scattered field u s . We first calculate the frequency components by FFT from the time-domain measurements (55) u s (x, κ n ) = F(U s (x, ·))(κ n ) , n = −N, ..., N − 1, x ∈ Γ h,A .
Then, we reconstruct the scattered field u s (x, κ n ) for x ∈ Q ⊂ R 3 where Q is a threedimensional grid. For every fixed x ∈ Q we evaluate The second part of the numerical implementation is the probing procedure. We chose some fixed constant = 0.2. For every point x in Q we evaluate the first hitting time T (x) and investigate U s (x, t) for t ∈ J := (T (x), T (x) + ). We set µ(x) = 1 for those points for which |U s (x, t)| > ρ in J, where ρ is some numerical threshold which we use to discriminate U s (x, t) = 0 and |U s (x, t)| > 0.
We present simulations where the unknown surface consists of a hill and a valley. We use frequencies from 0 to 6 with a stepsize h = 0.15. In Figure 4 we presents a visualization of the incident pulse in (a). The pictures (b), (c) and (d) show time slices of the scattered field at three different times t 1 , ..., t 3 . This confirms the arguments demonstrated in Figure  1.
The points which are identified as boundary points by the Time-Domain Probe Method are visualized in Figure 5. Figure (a) shows a slice plot along x 2 = 0. In (b) we show a horizontal view onto the reconstructions which proves that the location and height of the hill and valley are correctly found. Finally, a complete reconstruction of the height function is visualized in (c). The numerical results confirm the feasibility of the Time-Domain Probe Method.
The reconstructions here are comparable to frequency domain reconstructions for example by the point source method (c.f. [30], [31]) when we know and use the Dirichlet boundary condition. However, here we do not need to use this condition, so we are in a different setting where standard frequency domain algorithms do not work. A comparison to methods like the range test or the no-response test (compare [30]) needs to be part of future research.