Data recovery: from limited-aperture to full-aperture

Data recovery: from limited-aperture to full-aperture

Abstract

The inverse scattering problems have been popular for the past thirty years. While very successful in many cases, progress has lagged when only limited-aperture measurement is available. In this paper, we perform some elementary study to recover data that can not be measured directly. To be precise, we aim at recovering the full-aperture far field data from limited-aperture measurement. Due to the reciprocity relation, the multi-static response matrix (MSR) has a symmetric structure. Using the Green’s formula and single layer potential, we propose two schemes to recover full-aperture MSR. The recovered data is tested by a recently proposed direct sampling method and the factorization method. The numerical results show the possibility to recover, at least partially, the missing data and consequently improve the reconstruction of the scatterer.

Keywords: inverse scattering, multi-static response matrix, limited-aperture, data recovery.

AMS subject classifications: 35P25, 35Q30, 45Q05, 78A46

1 Introduction

The inverse scattering theory has been a fast-developing area for the past thirty years. The aim is to detect and identify the unknown objects using acoustic, electromagnetic, or elastic waves. Many methods have been proposed, e.g., iterative methods, decomposition methods, the linear sampling method, the factorization method and direct sampling methods [4, 5, 6, 7, 9, 11, 13, 14, 15, 18, 19, 23, 24]. Most of the above algorithms use full-aperture data, i.e., data of all the observation directions due to all incident directions. However, in many cases of practical interest, it is not possible to measure the full-aperture data, e.g., underground mineral prospection, mine location in the battlefield, and anti-submarine detection. Consequently, only limited-aperture data over a range of angles are available.

Various reconstruction algorithms using limited-aperture data have been developed [1, 3, 10, 12, 14, 17, 20, 21, 25]. Although uniqueness of the inverse problems can be proved in some cases [8], the quality of the reconstructions are not satisfactory. Indeed, limited-aperture data put forward a severe challenge for all the existing numerical methods. A typical feature is that the “shadow region” is elongated in down range [17]. Physically, the information from the “shadow region” is very weak, especially for high frequency waves [20]. For two-dimensional problems, the numerical experiments of the decomposition methods in [12, 25] indicate that satisfactory reconstructions need an aperture not smaller than 180 degrees.

Other than developing methods using limited-aperture data, we take an alternative approach to recover the data that can not be measured directly. As a consequence, methods using full-aperture data can be employed. We take the acoustic scattering by time-harmonic plane waves as the model problem. The measurement data are only available for limited-aperture observation angles but for all incident directions. The goal is to recover data for all observation angles. The case to recover full-aperture data from limited-aperture observation angles due to limited-aperture incident directions will be considered in future.

For scattering problems, it is well-known that the full-aperture data can be uniquely determined by the limited-aperture data. However, because of the severely ill-posed nature of the analytic continuation, it is in general not possible to recover full-aperture data using techniques such as extrapolation [2]. We take a different way by seeking an analytic function in a suitable space based on the PDE theory governing the scattering problem. More precisely, we look for kernels of layer potentials that generate the measured data approximately by regularization. Then these kernels are used to obtain the full-aperture data.

This paper is organized as follows. In Section 2, we briefly introduce the scattering problem of interests and the multi-static response (MSR) matrix, which is the far field pattern. Due to the reciprocity relation of the far field pattern, the MSR has a symmetry property, which can be used to recover partial missing data. In Section 3.1, we propose a technique using the Green’s formula to recover the full MSR. Another recovery technique based on the single layer potential is proposed in subsection 3.2. Combining these techniques and the symmetry property, a novel algorithm is proposed to recover the full-aperture MSR. In Section 4, numerical examples are presented to demonstrate the performance of the data recover techniques. The recovered data are tested using a direct sampling method and the factorization method. We draw some conclusions and discuss future works in Section 5.

2 The Multi-static Response Matrix

Let be the wave number of a time harmonic wave and be a bounded domain with Lipschitz-boundary such that the exterior is connected. Let the incident field be given by

 ui(x) = ui(x;d)=eikx⋅d,x∈Rn, (2.1)

where , denotes the direction of the plane wave.

The scattering problem for an inhomogeneous medium is to find the total field such that

 Δu+k2(1+q)u=0in Rn, (2.2) limr:=|x|→∞rn−12(∂us∂r−ikus)=0, (2.3)

where such that its imaginary part and in . The Sommerfeld radiation condition (2.3) holds uniformly with respect to all directions . If the scatterer is impenetrable, the direct scattering problem is to find the total field such that

 Δu+k2u=0in Rn∖¯¯¯¯Ω, (2.4) B(u)=0on ∂Ω, (2.5) limr:=|x|→∞rn−12(∂us∂r−ikus)=0, (2.6)

where denotes one of the following three boundary conditions:

 (1)B(u):=uon ∂Ω;(2)B(u):=∂u∂νon ∂Ω;(3)B(u):=∂u∂ν+λu% on ∂Ω

corresponding, respectively, to the cases when the scatterer is sound-soft, sound-hard, and of the impedance type. Here, is the unit outward normal to and is the (complex valued) impedance function such that almost everywhere on . Uniqueness of the scattering problems (2.2)–(2.3) and (2.4)–(2.6) can be shown with the help of Green’s theorem, Rellich’s lemma and unique continuation principle, see e.g., [8]. The proof of existence can be done by variational approaches (cf. [8, 22] for the Dirichlet boundary condition and [4] for other boundary conditions) or by integral equation methods (cf. [8]).

Radiating solutions of the Helmholtz equation have the following asymptotic behavior [14, 18]:

 us(x;d)=eiπ4√8kπ(e−iπ4√k2π)n−2eikrrn−12{u∞(^x;d)+O(1r)}as r:=|x|→∞ (2.7)

uniformly with respect to all directions . The complex valued function defined on the unit sphere is known as the scattering amplitude or far-field pattern with denoting the observation direction.

It is well known that the scatterer can be uniquely determined by the far field pattern for all [8]. Due to analyticity, for is uniquely determined by for if has a nonempty interior. Unfortunately, it is practically impossible to obtain on from on using the analytic continuation (see Atkinson [2]).

In this paper, we consider the discrete version of , i.e., the multi-static response (MSR) matrix in . Let ,

 di:=(cosθi,sinθi),i=1,2,…,2m,

and

 ^xj:=(cosθj,sinθj),j=1,2,…,2m.

The multi-static response (MSR) matrix is defined as

 Ffull:=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝u∞1,1u∞1,2⋯u∞1,2mu∞2,1u∞2,2⋯u∞2,2m⋮⋮⋱⋮u∞2m,1u∞2m,2⋯u∞2m,2m⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (2.8)

where for corresponding to observation directions and incident directions .

Assume that the far field pattern can only be measured in a limited-aperture. In particular, the measured data are the first columns of

 F(l)limit:=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝u∞1,1u∞1,2⋯u∞1,lu∞2,1u∞2,2⋯u∞2,l⋮⋮⋱⋮u∞2m,1u∞2m,2⋯u∞2m,l⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,1≤l<2m. (2.9)

The inverse problem considered in this paper is to firstly recover from , and then reconstruct the scatterer from the recovered . Note that is NOT symmetric, i.e., . Here and throughout the paper we use the superscript to denote the transpose of a matrix. We can partition the -by- MSR matrix into a -by- block matrix

 Ffull=(F11F12F21F22), (2.10)

where . The following theorem is a consequence of the reciprocity relation.

and .

Proof.

Recall that the far field pattern is the same if the direction of the incident field and the observation direction are interchanged [8], i.e.,

 u∞(^x,d)=u∞(−d,−^x),for all% ^x,d∈S1. (2.11)

For all , using the reciprocity relation (2.11), we have

 u∞i,j = u∞(^xj;di) = u∞(−di;−^xj) = u∞(−(cosθi,sinθi);−(cosθj,sinθj)) = u∞((cos(θi+π),sin(θi+π));(cos(θj+π),sin(θj+π))) = u∞((cosθi+m,sinθi+m);(cosθj+m,sinθj+m)) = u∞j+m,i+m,1≤i,j≤m.

Thus, we have .

Similarly, For all , using the reciprocity relation (2.11) again, we have

 u∞i,j+m = u∞(^xj+m;di) = u∞(−di;−^xj+m) = u∞(−(cosθi,sinθi);−(cosθj+m,sinθj+m)) = u∞((cos(θi+π),sin(θi+π));(cos(θj+m+π),sin(θj+m+π))) = u∞((cosθi+m,sinθi+m);(cosθj,sinθj)) = u∞j,i+m,1≤i,j≤m.

Thus, we have . The equality can be treated analogously. ∎

As a direct consequence of Theorem 2.1, the following data

 ˜F(l)limit:=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝u∞m+1,l+1u∞m+1,l+2⋯u∞m+1,2mu∞m+2,l+1u∞m+2,l+2⋯u∞m+2,2m⋮⋮⋱⋮u∞m+l,l+1u∞m+l,l+2⋯u∞m+l,2m⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,1≤l<2m, (2.12)

can be obtained directly from . Here, we have set if .

Remark 2.2.

If we set

 ˆFfull:=(F12F11F22F21),

then is symmetric, i.e., . The result also holds for phaseless MSR matrix.

3 Data Recover Schemes

In this section, we propose two methods to recover from . The first one is based on the Green’s formula. The second one is based on the single layer potential.

3.1 Method of Green’s Formula

Let be a bounded domain with connected complement such that and the boundary is of class . Let denote the unit normal vector to the boundary directed into the exterior of . The fundamental solution of the Helmholtz equation is given by

 Φ(x,y):=i4H(1)0(k|x−y|), (3.13)

where is the Hankel function of the first kind of order zero. The scattered field is a radiating solution to the Helmholtz equation in such that the Green’s formula holds [8]

 us(x;d)=∫∂B{us(y;d)∂Φ(x,y)∂ν(y)−∂us(y;d)∂ν(y)Φ(x,y)}ds(y),x∈R2∖¯¯¯¯B. (3.14)

Letting tend to the boundary and using the jump relations, it can be shown that

 (ϕ,ψ):=(us,∂us∂ν)∣∣∂B∈H1/2(∂B)×H−1/2(∂B)

solves the following boundary integral equations

 ϕ(x)=2∫∂B{ϕ(y)∂Φ(x,y)∂ν(y)−ψ(y)Φ(x,y)}ds(y),x∈∂B, (3.15) ψ(x)=2∂∂ν(x)∫∂B{ϕ(y)∂Φ(x,y)∂ν(y)−ψ(y)Φ(x,y)}ds(y),x∈∂B. (3.16)

For later use, we define the space

 W:={(ϕ,ψ)∈H1/2(∂B)×H−1/2(∂B):(ϕ,ψ) is a solution to(???)−(???).}.

Using the Green’s formula (3.14), has the following form (cf. [14])

 u∞(^x;d)=∫∂B{us(y;d)∂e−ik^x⋅y∂ν(y)−∂us∂ν(y;d)e−ik^x⋅y}ds(y),^x∈S1. (3.17)

Note that the Cauchy data is independent of the variable . From (3.17) the far field pattern can be computed in any direction if the Cauchy data is known.

Let be the measurement surface, which is an open subset of the unit sphere with nonempty interior (open relative to ). If we already know the far field pattern in , then it is natural to approximate the Cauchy data by solving the following integral equation

 F(ϕ(⋅;d),ψ(⋅;d))(^x)=u∞(^x;d),^x∈S10, (3.18)

where is defined by

 (3.19)
Theorem 3.1.

The operator is compact, injective with dense range in .

Proof.

The operator is certainly compact since its kernel is analytic in both variables.

Let satisfy in . By analyticity, we have

 ∫∂B{ϕ(y;d)∂e−ik^x⋅y∂ν(y)−ψ(y;d)e−ik^x⋅y}ds(y)=0,^x∈S1. (3.20)

Note that the left hand side of (3.20) is actually the far field pattern of the scattered field given by

 ws(x):=∫∂B{ϕ(y;d)∂Φ(x,y)∂ν(y)−ψ(y;d)Φ(x,y)}ds(y),x∈R2∖¯¯¯¯B.

By Rellich’s lemma and (3.20), vanishes in . Now, jump relations yield

 0=12ϕ(x)+∫∂B{ϕ(y)∂Φ(x,y)∂ν(y)−ψ(y)Φ(x,y)}ds(y),x∈∂B, 0=12ψ(x)+∂∂ν(x)∫∂B{ϕ(y)∂Φ(x,y)∂ν(y)−ψ(y)Φ(x,y)}ds(y),x∈∂B.

Recall that , which implies that is also a solution of (3.15)-(3.16). Hence, and is injective.

We consider the adjoint of and show that it is injective as well which proves the denseness of the range of . For all we extend by zero in to obtain . Recall the Herglotz wave function of the form

 vh(y):=∫S1eiky⋅^xh(^x)ds(^x),y∈R2.

Then we obtain that the adjoint operator is given by

 F∗h=(∂vh∂ν,−vh)∣∣∂B. (3.21)

Interchanging the order of integration, we have

 (F(ϕ,ψ),h)L2(S10) = ∫S10∫∂B{ϕ(y;d)∂e−ik^x⋅y∂ν(y)−ψ(y;d)e−ik^x⋅y}ds(y)¯¯¯¯¯¯¯¯¯¯¯h(^x)ds(^x) = ∫∂B⎧⎨⎩ϕ(y;d)∂¯¯¯¯¯¯¯¯¯¯¯¯vh(y)∂ν(y)−ψ(y;d)¯¯¯¯¯¯¯¯¯¯¯¯vh(y)⎫⎬⎭ds(y) = ⟨(ϕ,ψ),(∂vh∂ν,−vh)⟩ = ⟨(ϕ,ψ),F∗h⟩,

where the last two equalities hold in the sense of dual paring . Here and in the following, denotes the complex conjugate of .

We proceed by showing that the adjoint operator is injective. Let be such that on . Again extending by zero in to obtain . We find that the Cauchy data of the Herglotz wave function vanishes on . Note that the Herglotz wave function is an entire solution of the Helmhlotz equation in . Thus by Holmgren’s uniqueness theorem we deduce that vanishes identically in . This further implies that on [8] and the proof is complete. ∎

Method of Green’s Formula (MGF)

• Given partial data , set by Theorem 2.1.

• Solve for (3.18) using by Tikhonov regularization.

• Use (3.17) to obtain .

3.2 Method of Single Layer Potential

We consider the scattered field in the form of a single-layer potential

 us(x)=∫∂BΦ(x,y)ϕ(y)ds(y),x∈R2∖¯¯¯¯B

with an unknown density . Its far field pattern is given by

 u∞(^x)=∫∂Beik^x⋅yϕ(y)ds(y),^x∈S1. (3.22)

Inspired by this, we introduce the following integral equation of the first kind

 S∞ϕ=u∞, (3.23)

where the far field integral operator is defined by

 (S∞ϕ)(^x):=∫∂Be−ik^x⋅yϕ(y)ds(y). (3.24)

The properties of the far field integral operator have been collected in the following theorem.

Theorem 3.2.

The far field integral operator is injective and has dense range provided is not a Dirichlet eigenvalue for the negative Laplacian in .

We omit the proof since, except for minor adjustments, they literally coincide with those of Theorem 5.19 in [8], where, the 3D case is considered. The requirement on is not essential since we have the freedom to choose .

To recover full-aperture data , one may firstly solve the equation (3.23) by the Tikhonov regularization in and then insert the solution into to obtain the missed data. Again, the integral operator has an analytic kernel and therefore equation (3.23) is severely ill-posed.

Method of Single Layer Potential (MSLP)

• Given partial data , set by Theorem 2.1.

• Solve for (3.23) using by Tikhonov regularization.

• Use (3.22) to obtain .

3.3 From limited-aperture to full-aperture

The direct application of the above two methods does not produce good recovery of the full-aperture data due to the severe ill-posedness. To improve the result, we propose a step by step alternative technique which makes use of the symmetry of . Roughly speaking, we use MGF or MSLP to recover a few data using the known data. Then Theorem 2.1 is used to obtain more data. The process is repeated until is recovered.

Now we ready to introduce the algorithm to recover the full-aperture MSR.

• DR-MSR:

• Step 1. Measure the limited-aperture far field pattern .

• Step 2. Using MGF or MSLP, recover the far field pattern

 Mnew:={^xl+1,^xl+1,⋯^xl+t,^x2m−st−t+1,^x2m−st−t+2,⋯,^x2m−st},

i.e., compute data in new directions close to known data.

• Step 3. Recover in (2.9) using Theorem 2.1.

• Step 4. If is obtained, stop. Otherwise, set and go to Step 2.

Remark 3.3.

The scheme makes no use of the boundary conditions or topological properties of the underlying object . In other words, the full-aperture data is retrieved without any a priori information on .

4 Numerical Examples and Discussions

The numerical examples are divided into two groups. We first present some numerical examples to demonstrate how to use the proposed methods to recover data. The second group of numerical examples are to test the recovered data, which are used by some non-iterative methods to reconstruct the support of the scatterer. The results show that the reconstruction improves significantly, indicating that the recovered data does help in certain cases. Two scatterers are considered (see Fig. 1):

 Kite: x(t) =(cost+0.65cos2t−0.65,1.5sint),0≤t≤2π, (4.1) Peanut: x(t) =√3cos2t+1(cost,sint),0≤t≤2π. (4.2)

The synthetic scattering data are generated by the boundary integral equation method, i.e., , for equidistantly distributed directions in . We first verify Theorem 2.1. Let , and take the kite as an example. The MSR matrix is . The four block matrices are given as follows.

 F11 = ⎛⎜ ⎜ ⎜⎝−2.6282+1.8817i0.1698+0.4158i0.1657−0.2286i1.0722−0.6313i0.0028−0.9694i−2.5830+1.9160i0.3264+0.1581i−0.4424−0.9227i−0.0740+0.7809i0.2839−0.5024i−2.4052+1.5689i0.3264+0.1581i0.1929−0.5886i0.2202+0.4858i0.2839−0.5024i−2.5830+1.9160i⎞⎟ ⎟ ⎟⎠, F22 = ⎛⎜ ⎜ ⎜⎝−2.6282+1.8817i0.0028−0.9694i−0.0740+0.7809i0.1929−0.5886i0.1698+0.4158i−2.5830+1.9160i0.2839−0.5024i0.2202+0.4858i0.1657−0.2286i0.3264+0.1581i−2.4052+1.5689i0.2839−0.5024i1.0722−0.6313i−0.4424−0.9227i0.3264+0.1581i−2.5830+1.9160i⎞⎟ ⎟ ⎟⎠, F12 = ⎛⎜ ⎜ ⎜⎝−0.5250+0.1132i1.0722−0.6313i0.1657−0.2286i0.1698+0.4158i1.0722−0.6313i−0.0050−0.3054i0.1510+0.3285i−0.5603−0.0594i0.1657−0.2286i0.1510+0.3285i−0.3128−0.5104i−0.6526−1.3338i0.1698+0.4158i−0.5603−0.0594i−0.6526−1.3338i−0.0441−0.9080i⎞⎟ ⎟ ⎟⎠, F21 = ⎛⎜ ⎜ ⎜⎝−0.4473−0.3633i0.1929−0.5886i−0.0740+0.7809i0.0028−0.9694i0.1929−0.5886i−0.0441−0.9080i−0.6526−1.3338i−0.5603−0.0594i−0.0740+0.7809i−0.6526−1.3338i−0.3128−0.5104i0.1510+0.3285i0.0028−0.9694i−0.5603−0.0594i0.1510+0.3285i−0.0050−0.3054i⎞⎟ ⎟ ⎟⎠.

It is obvious that Theorem 2.1 holds:

 F11=FT22,F12=FT12 and F21=FT21.

In the rest of the section, we fix and divide uniformly into directions. Assuming that the incident directions cover the full aperture, the observation directions only span a subset of . In particular, let with the observation angle . Then we consider the measurements for three cases:

 (1)ϕ∈(0,π/2),(2)ϕ∈(0,2π/3),and(3)ϕ∈(0,π).

Namely, we have the synthetic data for , , and , respectively. Then is perturbed by random noises

 F(l),δlimit=F(l)limit+δ∥F(l)limit∥R1+R2i∥R1+R2i∥,

where and are two matrixes containing pseudo-random values drawn from a normal distribution with mean zero and standard deviation one. The value of is the noise level, which is taken as .

4.1 Data Recovery

Examples in this subsection are to test the validity of data recover algorithms proposed in Section 3. Given , the goal is to recover , the full-aperture MSR. We take the kite as the scatterer. The artificial domain is chosen to be a disc centered at the origin with radius .

Figures 2-4 show the data reconstructions with different measurement apertures and respectively. We show the recovered data for the incident direction , i.e., the first row of in Figures 2-4. Only few data are well reconstructed. This is reasonable because both methods, MGF and MSLP, involve solving severely ill-posed integral equations. In particular, the symmetry does not apply in this case. Similar results with respective to are shown in Figures 2-4.

For incident angles in , we have obtained nearly exact far field pattern for all observation directions. In Figures 2-4 we show results for two incident directions and . This further verify the symmetric structure of the multi-static response matrix.

4.2 Applications in Sampling Methods

We first test the recovered data by a novel direct sampling method (DSM) proposed in [18], which uses an indicator functional defined as

 I(z):=|ϕ(z;−d)FfullϕT(z;^x)|2, (4.3)

where and . The indicator takes its maximum on or near the boundary of the scatterer. Consequently, the plot of the indicator can be used to reconstruct the scatterer.

This method can be modified to use only limited-aperture data by introducing

 Ilimit(z):=|ϕ(z;−d)F(l)limitϕTlimit(z;^x)|2, (4.4)

where corresponds to the limited-aperture observation directions and is the limited-aperture data given by (2.9).

Denote by and the recovered full-aperture data using MGF and MSLP, respectively. We introduce the following indicator

 I(ii)full(z):=|ϕ(z;−d)F(ii)fullϕT(z;^x)|2,ii=2,3. (4.5)

Alternatively, one can reconstruct the scatterer by , using recovered full-aperture data. As will seen shortly, the quality of the reconstructions indeed improves.

We used a grid of equally spaced sampling points on some rectangle . For each point , we compute the indicator functionals given in (4.4)-(4.5).

The typical feature for limited-aperture problems is that the concave part cannot be reconstructed if the the observation angles do not cover the concave part of the obstacle. A common criterion for judging the quality of a reconstruction method is whether the concave part of the obstacle can be successfully recovered. The resulting reconstructions by using the indicator functional with limited-aperture far field patterns are shown in Figures 5, , and for different observation apertures. Clearly, the quality improves with the increase of observation apertures. We also observe that the illuminated part is well constructed, but the shadow region is elongated down range.

As shown in the second and third columns of Figure 5, the reconstructions are indeed improved by using the data recover techniques. In particularly, the two wings of the kite appear and the shadow region is reconstructed very well. Considering the severe ill-posedness of the data reconstruction of an analytic function and the relative noise level , the target reconstructions given in Figure 5 are satisfactory. Similar results are shown in Figures 6 for the peanut.