Inverse Scattering from a Sound-Hard Crack via Two-Step Method

Abstract

We present a two-step method for recovering an unknown sound-hard crack in ℝ² from the measured far-field pattern. This method, based on a two-by-two system of nonlinear integral equations, splits the reconstruction into two consecutive steps which consists of a forward and an inverse problems. In this spirit, only the latter needs to be regularized.

1. Introduction

In this paper we consider an acoustic scattering problem from a sound-hard crack. This problem is modeled by an exterior boundary value problem governed by Helmholtz equation for an open arc with Neumann boundary conditions on both sides of the arc. Our major concern is the inverse problem which aim is to reconstruct the crack from some measurements. This kind of problem is of fundamental importance, for example, in material investigation, nondestructive testing, or in seismic exploration.

The inverse scattering problem for an open arc was first investigated by Kress [1]. In his article, integral equation method was used to solve both the direct and inverse problems for a sound-soft crack. The scattering problem in the unbounded domain was then converted into a boundary integral equation. Mönch [2] extended this approach to a Neumann problem. Their integral-equation-based method was to solve the so-called far-field equation. 𝐹(Γ)=𝑢_∞, (1.1) for the unknown crack Γ from the measured far-field data. Because of the nonlinearity and the compactness of the far-field operator, both linearization and regularization are needed to keep the solving of this equation in the framework of linear regularization theory. One of the major drawbacks of the regularized Newton’s method used in [1, 2] is that at every iteration, a direct problem with different boundary data must be solved for the Fréchet derivative of this far-field operator which is an essential part of this method.

The nonlinear integral equations method, proposed by Kress and Rundell [3], can be seen as a remedy for the Newton’s method. They proposed a two-by-two system of nonlinear integral equations which was shown to be equivalent to the original inverse problem. The Fréchet derivatives were obtained simply through the solving of this system. This method was extended to other boundary conditions and cracks in [4, 5] and to obstacle scattering in [6].

Another group of method which is not iterative is the decomposition method (see [7–9]). The idea is to split the nonlinearity and the ill-posedness of the original inverse problem. The scattered field 𝑢^𝑠 is computed in the first linear ill-posed step from the far-field pattern 𝑢_∞. In the second step, which is nonlinear, the unknown crack Γ is reconstructed as the location where the boundary condition fits in the least square sense. This method is thus carried out without computation of the direct problem. One of its disadvantages is that the reconstruction is less accurate than that of the Newton’s method.

In [10], a method consisting of two steps was proposed for a sound-soft crack. It is a mixture of the above mentioned methods. In this paper, we’d like to extend this method to the case of a sound-hard crack and to the case of limited aperture. The plan of the paper is as follows. For the sake of completeness and also for the introduction of notations, we will briefly summarize the main results of the direct problem in Section 2. In Section 3 we will consider the inverse scattering problem in its equivalence form of a two-by-two system of integral equations. In Section 4 we will discuss an iteration scheme for numerical computation of the inverse problem. This will be followed by some numerical examples in the final section.

2. Direct Neumann Problem

Given a regular nonintersecting 𝐶³-smooth open arc Γ⊂ℝ² which can be represented as Γ=𝑧(𝑠)∶𝑠∈[−1,1],𝑧∈𝐶³[−1,1],𝑧^′(𝑠)≠0,∀𝑠∈[−1,1]. (2.1) The two end points of the crack are denoted by 𝑧^∗₋₁, 𝑧^∗₁ respectively. The left hand side and the right-hand side of the crack are written by Γ₊ and Γ₋, respectively. The unit outward normal to Γ₊ is denoted by 𝜈. Further we set Γ₀∶=Γ⧵{𝑧^∗₋₁,𝑧^∗₁}.

The direct scattering problem for a sound-hard crack that we are considering is as follows.

Problem 1 (DP). Given an incident plane wave 𝑢^𝑖(𝑥,𝑑)∶=𝑒^{𝑖𝑘⟨𝑥,𝑑⟩} with a wave number 𝑘>0 and a unit vector 𝑑 giving the direction of propagation, we find a solution 𝑢^𝑠∈𝐶²(ℝ²⧵Γ)∩𝐶(ℝ²⧵Γ₀) to the Helmholtz equation Δ𝑢^𝑠+𝑘²𝑢^𝑠=0,inℝ²⧵Γ, (2.2) which satisfies the Neumann boundary conditions: 𝜕𝑢^𝑠_±𝜕𝜈=−𝜕𝑢^𝑖𝜕𝜈onΓ₀, (2.3) on both sides of the crack and the Sommerfeld radiation condition: lim_𝑟→∞√𝑟𝜕𝑢^𝑠𝜕𝑟−𝑖𝑘𝑢^𝑠=0,𝑟∶=|𝑥| (2.4) uniformly for all directions 𝑥∶=𝑥/|𝑥|.

In (2.3), the limits 𝜕𝑢^𝑠_±(𝑥)𝜕𝜈∶=lim_ℎ→0𝜈(𝑥),grad𝑢^𝑠(𝑥±ℎ𝜈(𝑥)),𝑥∈Γ₀ (2.5) are required to exist in the sense of locally uniform convergence.

Note that the boundary conditions (2.3) can be reformulated as homogeneous Neumann conditions for the total field 𝑢∶=𝑢^𝑖+𝑢^𝑠, that is, 𝜕𝑢_±/𝜕𝜈=0.

Using boundary integral equations, this direct problem can be solved via the layer approach. We refer to the monograph [11] for the details. In terms of the fundamental solution to the Helmholtz equation in ℝ²Φ𝑖𝑥,𝑦∶=4𝐻₀⁽¹⁾𝑘𝑥−𝑦,𝑥≠𝑦, (2.6) the following theorem, which was proved in [12], ensures the unique solvability of the direct scattering problem.

Theorem 2.1. The direct Neumann Problem 1 has a unique solution given by 𝑢^𝑠(𝑥)=_Γ𝜕Φ𝑥,𝑦𝑦𝜑𝑦𝑦𝜕𝜈𝑑𝑠,𝑥∈ℝ²⧵Γ, (2.7) where 𝜑∈𝐶^1,𝛼,∗(Γ) is the (unique) solution to the following integral equation 𝜕𝜕𝜈(𝑥)_Γ𝜕Φ𝑥,𝑦𝑦𝜑𝑦𝑦𝜕𝜈𝑑𝑠=−𝜕𝑢^𝑖(𝑥)(𝑥).𝜕𝜈 (2.8) The function space 𝐶^1,𝛼,∗(Γ) is defined by 𝐶^1,𝛼,∗(Γ)∶=𝜑∣𝜑𝑧^∗₋₁=𝜑𝑧^∗₁(𝑧(𝑠=0,𝑑𝜑))=𝜑()𝑑𝑠arccos𝑠√1−𝑠²,𝜑∈𝐶^0,𝛼[0,𝜋], (2.9) for 0<𝛼<1.

At this place, we introduce the far-field pattern 𝑢_∞ of the scattered field 𝑢^𝑠. The far-field pattern describes the behavior of the scattered wave at the infinity 𝑢^𝑠(𝑥)=𝑒^{𝑖𝑘|𝑥|}|𝑥|𝑢_∞()1𝑥+𝑂|𝑥||𝑥|⟶∞ (2.10) uniformly for all directions 𝑥∈Ω∶={𝑥∈ℝ²∣|𝑥|=1}. The one-to-one correspondence between radiating wave and its far-field pattern is established by the Rellich’s lemma. In the case of a sound-hard crack, the far-field pattern is given by 𝑢_∞(𝑥)=𝜌_∞Γ𝜈𝑦,𝑥𝑒^{−𝑖𝑘⟨}𝑥_,𝑦⟩𝜑𝑦𝑑𝑠𝑦, (2.11) with the density function 𝜑 from Theorem 2.1 and the constant 𝜌_∞=(𝑘/8𝜋)𝑒^{−𝑖(𝜋/4)}. We note that from the viewpoint of the inverse scattering problems, the far-field pattern is of particular interest since usually one has no access to the close neighborhood of the crack.

It is convenient for our further treatment to rewrite our integral equation (2.8) in operator form. For this purpose, we introduce the following integral operators: 𝑆𝜑(𝑥)∶=_ΓΦ𝜑𝑦𝑦(𝑥)𝜕𝑥,𝑦𝑑𝑠,𝑥∈Γ,𝑇𝜑∶=𝜕𝜈(𝑥)_Γ𝜕Φ𝑥,𝑦𝑦𝜑𝑦𝑦𝐹𝜕𝜈𝑑𝑠,𝑥∈Γ,_∞𝜑(𝑥)∶=𝜌_∞Γ𝜈𝑦,𝑥𝑒^{−𝑖𝑘⟨}𝑥_,𝑦⟩𝜑𝑦𝑑𝑠𝑦,𝑥∈Ω. (2.12) The operator 𝐹_∞(𝜑), which is termed the far-field operator in the literatures, plays an important role in the scattering theory. In terms of these operators, the solving of the direct scattering problem amounts to the solving of the following system of two operator equations: 𝑇𝜑=−𝜕𝑢^𝑖,𝑢𝜕𝜈_∞=𝐹_∞𝜑. (2.13) In words, given the crack, the unknown density function 𝜑 is firstly solved from the first equation. Inserting 𝜑 into the very definition of the far-field operator, one obtains the far-field pattern 𝑢_∞ from the second equation of the system (2.13) straight forward. This reasoning forms the basic idea of our numerical scheme for the inverse problem. At this place, we want to reduce the hypersingularity of the operator 𝑇. Using Maue’s identity, the hypersingular operator 𝑇 is splitted into two parts (Theorem 7.29 in [11]): 𝜕𝑇𝜑=𝑆𝜕𝜗𝜕𝜑𝜕𝜗+𝑘²𝜈,𝑆𝜑𝜈, (2.14) where 𝜗 is the unit tangent vector.

3. Inverse Scattering Problem

After introducing the notations in the last section, we consider the following inverse problem.

Problem 2 (IP). Determine the crack Γ if the far-field pattern 𝑢_∞ is known for one incident wave.

For the uniqueness of this inverse problem, that is, for the identifiability of the arc, we refer to [2] where infinitely many incident waves are used. Although it is not yet proven, it is widely believed that the scatterer can be uniquely determined from one incident wave. Since our numerical method is based on Newton’s method which has the advantage of dealing with only one incident wave well, this leads to our problem setting. For the general solution theory of the inverse scattering problem, we refer to the monograph [13].

Motivated by the method of nonlinear integral equations in [3], a system of nonlinear integral equations which solution is the solution to the inverse problem was presented in [5] for a Neumann crack. In [5], the following system was proposed: 𝐵Γ,𝜑=−𝜕𝑢^𝑖𝜕𝜈_Γ,𝐹_∞Γ,𝜑=𝑢_∞, (3.1) where the operators 𝐵 and 𝐹_∞ are no more than the operators 𝑇 and 𝐹_∞ in (2.13), respectively. The additional parameter Γ refers to the dependence on the crack. In [5, Theorem 2], the equivalence of the inverse problem and the solving of the system (3.1) was proven. The basic idea of this system is that both the direct and inverse problem share the same integral equations, with slightly different interpretations. This makes the inverse problem really inverse to the original physical problem.

In [5], (3.1) was treated as a coupled system with two nonlinear ill-posed equations. This means that the two equations in (3.1) have to be linearized and regularized at the same time, and at every iterative step. The Fréchet derivatives for every operator both w.r.t. the unknown boundary and the unknown density have to be calculated at every step. Besides, two regularization parameters have to be selected for the scheme.

Start the inverse problem with the same system (3.1), we suggest a different scheme in the next section in details. Before formulating our algorithm, we need to parametrize the system.

As in Section 2, we use the parametrization Γ={𝑧(𝑡)∶𝑡∈[−1,1]} for the sound-hard crack. To incorporate the square root singularities of the solution 𝑢^𝑠 at the crack tips, we use the cosine substitution 𝑡=cos𝜏, as suggested in [14]. Keeping the notations in [5], we write 𝛾(𝜏)∶=𝑧(cos𝜏) for the crack Γ and 𝜓(𝜏)=𝜑(cos𝜏) for the unknown density.

The system of integral equation (3.1), using Maue’s identity, is now parametrized as follows: 𝐴₀𝛾,𝜓^′−𝐴₁𝛾𝐴𝛾,𝜓=𝑎,_∞𝛾,𝜓=𝑢_∞, (3.2)(3.3) where 𝐴₀𝛾,𝜓^′∶=^𝜋₀𝐾₀(𝜏,𝜎)𝜓^′(𝜎)𝐴d𝜎,₁𝛾,𝜓∶=^𝜋₀𝐾₁()𝜓(𝜎)𝐴𝜏,𝜎sin𝜎d𝜎,_∞𝛾,𝜓∶=𝜌_∞^𝜋₀𝐾_∞()𝜓(𝜎)𝑎𝛾⟨𝑛(𝜏)⟩𝑢𝑥,𝜎sin𝜎d𝜎,∶=−2𝑖𝑘,𝑑^𝑖𝛾(𝜏),𝑑, (3.4) with the kernels 𝐾₀(𝜏,𝜎)∶=𝑖𝑘4𝐻^(1)′₀𝑘𝛾(𝜏)(𝜎𝛾(𝜏)(𝜎)(𝜏−𝛾)−𝛾,𝛾′)𝛾(𝜏)(𝜎,𝐾−𝛾)₁(𝜏,𝜎)∶=𝑖𝑘²4𝛾^′(𝜏),𝛾′(𝜎)𝐻₀⁽¹⁾𝑘𝛾(𝜏)(𝜎,𝐾−𝛾)_∞(𝑥,𝜎)∶=⟨𝑛(𝜎),𝑥⟩𝑒^{−𝑖𝑘⟨𝛾(𝜎),}𝑥_⟩ (3.5) and the parametrized outward normal 𝑛=(𝜈∘𝛾)⋅|𝛾′|.

4. A Two-Step Method

Now we want to introduce our two-step algorithm. The parametrized systems (3.2), (3.3), or the original system (3.1) can be interpreted differently as in [5]. The first equation of the system can be seen as the equation for the direct problem. To be more precise, given the crack Γ, which is equivalent to fixing the variable 𝛾, (3.2) is solved for 𝜓. Hence this equation is well posed according to Theorem 2.1. Moreover, it is a linear integral equation in 𝜓. Another advantage is numerically as it can be just solved by the direct solver. With this 𝜓 at hand, we can solve (3.3) for 𝛾. This is a nonlinear ill-posed equation which we will linearize and regularize. Based on Newton’s method, the Fréchet derivative of 𝐴_∞ w.r.t. 𝛾 is needed for the linearization. The Fréchet derivative of the integral operator is simply given by the Fréchet derivative of its kernel (cf. [15]). For brevity, we set 𝐚^⟂=(𝑎₂,−𝑎₁)^𝑡 if 𝐚=(𝑎₁,𝑎₂)^𝑡. We have 𝐴^′_∞𝛾,𝜓;𝑞∶=𝜌_∞^𝜋₀𝐾^′_∞𝛾,𝜓;𝑞𝜓(𝜎)sin𝜎d𝜎, (4.1) where 𝐾^′_∞𝛾,𝜓;𝑞=𝑞′(𝜎)^⟂,𝑥−𝑖𝑘𝑞(𝜎),𝑥⟨𝑛(𝜎),𝑥⟩𝑒^{−𝑖𝑘⟨𝛾(𝜎),}𝑥_⟩. (4.2)

The two-step scheme reads 𝐴₀𝛾,𝜓^′−𝐴₁𝛾𝐴𝛾,𝜓=𝑎,_′∞𝛾,𝜓𝑞+𝐴_∞𝛾,𝜓=𝑢_∞. (4.3)(4.4) Because of the compactness of the operator 𝐴^′_∞, (4.4) still has to be regularized. To accomplish this, we apply the Tikhonov regularization, that is, instead of (4.4), the equation 𝛼𝐼+𝐴_∞^′∗𝛾,𝜓𝐴^′_∞𝛾,𝜓𝑞=𝐴_∞^′∗𝛾,𝜓𝑢_∞−𝐴_∞𝛾,𝜓 (4.5) has to be solved with a regularization parameter 𝛼>0.

At this point, we want to discuss the uniqueness and the solvability of the system (4.3), (4.5). From the solution theory of the direct problem in Section 2, (4.3) is uniquely solvable. However, the Fréchet derivative 𝐴^′_∞ is not injective as pointed out in [1]. The null space of this operator is given by 𝑁𝐴^′_∞=𝑞∶𝜈(𝜏)⋅𝑞(𝜏)=0,𝜏∈[0,𝜋]. (4.6)

This reflects the fact that different parametrizations of the arc leading to the same set of points turn out to have the same far-field pattern. We can avoid this ambiguity by limiting our solution space to the set of arcs representable as the graph of a function as suggested in [1]. Once this restriction is made, the operator 𝐴^′_∞ is then an injective linear compact operator. As a consequence of the regularization theory, (4.5) is uniquely solvable, see [11]. Again, (4.3) is as a direct problem uniquely solvable as stated in Theorem 2.1. Thus, we can have the following theorem.

Theorem 4.1. If the pairs (𝜓,𝑞) and (𝜓,𝑞) are solutions to the system (4.3), (4.5), then 𝜓=𝜓 and 𝑞=𝑞.

Numerically, we solve the systems (4.3), (4.5) in two steps: At the start, given an initial guess 𝛾₀ for the unknown crack Γ. The first equation is treated as the corresponding direct problem and is solved for the density 𝜓. Then the second equation updates the crack by means of the regularized Newton’s method.

The algorithm for our method can be formulated as follows:(i)given an initial guess 𝛾₀ for the unknown crack;(ii)iterative steps, for 𝑛=0,1,2,…. (1)Step 1: solve 𝐴₀𝛾_𝑛,𝜓^′−𝐴₁𝛾_𝑛,𝜓=𝑎𝛾_𝑛, (4.7) for 𝜓.(2)Step 2: solve 𝛼𝐼+𝐴_∞^′∗𝛾_𝑛,𝜓𝐴^′_∞𝛾_𝑛,𝜓𝑞_𝑛=𝐴_∞^′∗𝛾_𝑛,𝜓𝑢_∞−𝐴_∞𝛾_𝑛,𝜓, (4.8) for the update 𝑞_𝑛 of 𝛾_𝑛.(iii)𝛾_𝑛+1=𝛾_𝑛+𝑞_𝑛. (4.9)At this place we comment that from the numerical point of view our scheme is very attractive. As a variant of the nonlinear integral equations method, the Fréchet derivative of the integral operator is computed directly by solving the integral equation (4.8). Besides, the derivative is very simple as compared to those in [5] since only the far-field operator which has a smooth kernel is differentiated. Numerically it can be easily solved via the rectangular rule, for example. Another advantage of this method is that it splits the problem into two smaller parts and thus makes the computation cheaper. However, being a Newton’s method, the convergence of our numerical scheme is still open.

5. Numerical Results

In this section we will demonstrate the applicability of our method via some examples. We reconstruct the unknown crack from the knowledge of the far-field pattern at a number of points resulted from just one incident wave. For the direct problem, the forward solver is applied with 63 collocation points. In order to avoid committing an inverse crime, the number of collocation points used in the inverse problem is chosen to be different from that of the forward solver. The point is, in the inverse problem, (4.3) is solved with the same forward solver as for the direct problem. We therefore choose 31 collocation points in the inverse algorithm. In all our examples, we take the incident wave coming from the direction 𝑑=(1/√2)(−1,1). The far-field pattern is measured at 15 points evenly distributed on the unit circle.

The basis functions for the parametrization are taken from the space: 𝑉_𝑚=span𝑇₀,𝑇₁,…,𝑇_𝑚−1, (5.1) where 𝑇^′_𝑘𝑠 are the kth Chebyshev monomials 𝑇_𝑘(𝑥)=cos(𝑘cos⁻¹𝑥). The selection of the Chebyshev polynomials is based on the fact that they can take care of the square root singularity in their own right. Thus the updates for the crack can be represented in the form: 𝑞(𝑥)=^𝑚−1_𝑘=0𝑏¹_𝑘𝑇_𝑘(𝑥),^𝑚−1_𝑘=0𝑏²_𝑘𝑇_𝑘(𝑥), (5.2) with the unknown coefficients 𝑏¹_𝑘, 𝑏²_𝑘, 𝑘=0,…,𝑚−1 which have to be determined from (4.5). We test our scheme for two different wave numbers, 𝑘=1,3. In the case where 𝑘=1, 𝑚 is taken to be 4. For 𝑘=3, we choose 𝑚=5. The starting curve, that is, the initial guess for the regularized Newton’s method, is taken to be the straight line 𝑦=0 in all examples. According to the discrepancy principle, the stopping criterion for the iterative scheme is given by the relative error: 𝑢_∞−𝑢_∞,𝑛2𝑢_∞2≤𝜖, (5.3) which is taken to be 0.003 in case of exact data and 0.03 in the case where 3% data error are present. In all our figures, the dotted line (blue) represents the initial guess. We denote by the dashed line (red) the true solution and by the solid line (black) the reconstruction.

Example 5.1. For the first example, we take the arc 𝑥Γ=𝑥,exp³2,𝑥∈[−1,1], (5.4) which is representable of a graph of a function. The numerical results are given in the Figures 1, 2 for 𝑘=1 and 𝑘=3, respectively.
From the numerical results, we see that the reconstructions for exact data are very good. The reconstructions in the case where random are present, the reconstructions are not bad.

(a) Exact data, 𝑁=47

(b) 3% error, 𝑁=31

Figure 1 

𝑘=1.

(a) Exact data, 𝑁=62

(b) 3% error, 𝑁=44

Figure 2 

𝑘=3.

Example 5.2. To demonstrate the benefit of our method, we choose a curve which does not belong to our solution space. To this end, we take the arc Γ=2sin3𝜋84𝑥+3,−sin3𝜋44𝑥+3,𝑥∈[−1,1], (5.5) which is not representable as a function.
The results are shown in Figures 3 and 4.

(a) Exact data, 𝑁=37

(b) 3% error, 𝑁=43

Figure 3 

𝑘=1.

(a) Exact data, 𝑁=24

(b) 3% error, 𝑁=25

Figure 4 

𝑘=3.

We see that the reconstructions are very good, even in the case of erroneous data.

Example 5.3. In this final example, we choose the same curve as in the last example for the case of limited aperture. We assume that the data are only measurable within a certain range apart from the incident angle. Figure 5 shows the results for uniformly distributed far-field measurements within 90 degrees from both sides of the incident wave. We see that in this case, the reconstructions are very accurate both for the exact and erroneous data. In Figure 6, the results are shown for measurements taken within 45 degrees from both sides of the incident wave. We see that the reconstructions are convincible for only 15 data points.

(a) Exact data, 𝑁=59

(b) 3% error, 𝑁=34

Figure 5 

𝑑±𝜋/2.

(a) Exact data, 𝑁=60

(b) 3% error, 𝑁=50

Figure 6 

𝑑±𝜋/4.

We conclude this paper with some remarks. First of all, our examples above show the feasibility of the proposed numerical method. On one hand, being a Newton-type method, our method is conceptually simple and numerically more accurate than the traditional decomposition method. On the other hand, being a variant of the nonlinear integral equations method, the derivative is directly computable from the algorithm itself which makes the method more easily accessible than the classical Newton’s method. The splitting of the problem into two smaller parts makes our method more competitive as the computational cost concerns. Finally we want to point out that our method can be carried over to other type of boundary conditions and also to other type of scattering problems.

Acknowledgment

This work is partially supported by the NSC Grant NSC-100-2115-M-006-003-MY2.