Iterative and range test methods for an inverse source problem for acoustic waves

We propose two methods for solving an inverse source problem for time-harmonic acoustic waves. Based on the reciprocity gap principle a nonlinear equation is presented for the locations and intensities of the point sources that can be solved via Newton iterations. To provide an initial guess for this iteration we suggest a range test algorithm for approximating the source locations. We give a mathematical foundation for the range test and exhibit its feasibility in connection with the iteration method by some numerical examples.


Introduction
The scattering of time harmonic acoustic point sources at n source points with locations s j and intensities c j for j = 1, . . . , n at a sound-sound soft obstacle D ⊂ R 2 can be modeled by the solution of subject to the Dirichlet boundary condition u = 0 on ∂D and the Sommerfeld radiation condition at infinity. Here, δ s stands for the Dirac delta distribution. This paper is concerned with the inverse problem to recover the location and the intensity of the sources from knowledge of the normal derivative of u on the boundary ∂D where we assume the boundary ∂D to be known. We propose an iterative method based on the reciprocity gap principle, that is, on Green's integral theorem. Our approach is motivated by a series of papers using similar ideas for the iterative solution of inverse problems for the shape of the boundary that was initiated by Kress and Rundell [9] for an inverse boundary value problem for the Laplace equation. For a survey of the extension of this approach to the Helmholtz equation we refer to Ivanyshyn et al [6]. For the case of the inverse source problem the reciprocity gap principle leads to a nonlinear equation for the location and the intensities that can be solved iteratively via simultaneous linearization with respect to all unknowns. Our approach modifies the method considered in [14] through the use of point sources on the boundary ∂D rather than point sources in the exterior R 2 \D of the scatterer. This results in more accurate reconstructions as indicated through our numerical examples. The description of the iterative scheme is provided in section 3.
In order to provide an initial guess for the iterative scheme we propose a range test algorithm for finding approximations for the source locations. For a closed curve containing the closureD of the scatterer in its interior we construct a compact linear operator A : L 2 ( ) → L 2 (∂D) such that the range of A can be used as an indicator whether the source locations are contained in the annulus between and ∂D or not. By choosing different shapes for the curve and numerically deciding via Tikhonov regularization on the solvability of the ill-posed operator equation Aϕ = f for appropriate right-hand sides f , depending on the given normal derivative of u, it is possible to approximate the source locations. We will give a detailed analysis of this algorithm in section 4 including the construction of the operator A and establishing injectivity and dense range as prerequisites for using Tiknonov regularization in the range test. In section 5, numerical examples illustrate the feasibility of the range test and its combination with the iterative scheme. In principle, the approach can of course be extended to three dimensions.
Range test and probe methods have been more recently suggested and developed in inverse scattering for gaining information on the location and shape of scatterers by Luke and Potthast [11] and by Potthast et al [13]. For a survey we refer to Potthast [12]. However, to our knowledge, this type of methods has not yet been employed for inverse source problems.
In practical remote sensing, faraway sources radiate fields that, within measurement precision, are nearly those radiated by point sources. Hence, the inverse source problem occurs for example in astronomy. If instead of the exterior problem we consider the interior problem with unknown source locations and intensities within D, then the inverse source problem has applications, for example, in electro-and magneto-encephalography. Algorithms like the MUSIC, in principle, also address an inverse source problem but with a different set of data (see, for example, Kirsch [7]). For the related work on the type of inverse problem that we are addressing we refer to [1,4,10].

The inverse source problem
We proceed with a more specified description of the inverse source problem under consideration. Let D ⊂ R 2 be a simply connected bounded domain with a C 2 boundary ∂D and outward unit normal ν. Denote by the fundamental solution to the Helmholtz equation with positive wave number k given in terms of the Hankel function of the first kind and of order zero. Consider the direct scattering problem for the sound-soft obstacle D with the incident field u i generated by a source distribution with source points s j in R 2 \D and intensities c j ∈ C\{0} for j = 1, . . . , n. The scattered field u s has to satisfy the Helmholtz equation the sound-soft boundary condition and the Sommerfeld radiation condition uniformly for all directions. In the distributional sense, the total field u := u i + u s satisfies the inhomogeneous Helmholtz equation with the inhomogeneity c j δ s j (2.6) given in terms of the Dirac delta distribution δ s .
The inverse source problem we are interested in is to determine the source locations and intensities, including their number n, from knowledge of the normal derivative of the total field. We note that this inverse problem is clearly nonlinear, since the total field u and therefore its normal derivative g depend nonlinearly on the source locations s j . Moreover, the problem is also ill-posed, since due to the well-posedness of the direct scattering problem for the sound-soft scatterer D, given any ε > 0, we can add a source term with intensity one located faraway from the boundary ∂D to f in (2.6) such that the difference of the normal derivatives corresponding to the original source f and the perturbed source in the L 2 norm is smaller than ε. In this sense, the solution to the inverse problem does not depend continuously on the data.
To be more precise, we note that while considering noisy data there are some classical situations where it is virtually impossible to distinguish between two different sources aligned far from the measurement set. Assume that we have two sources at the points s 0 = Rd and s ε = (R + ε)d, where d is any direction, R > 0 is large and ε relatively small. While measuring these incident waves on a boundary relatively near the origin, for some intensities the difference may in practice be indistinguishable. From the expansions choosing the intensity c ε = c 0 e −ikε 1 + ε/R for the total incident field we just have and this may completely be overwritten by the noise on the measurements. We cite the uniqueness result for the inverse problem as considered in section 3.2 of [14] with a sketch of its proof.
Before we proceed with describing our solution algorithm, we note that in the case of the Laplace equation there exists an interesting simplification that connects the exterior inverse source problem to an interior inverse source problem. Denote by the fundamental solution for the Laplace equation and consider a disk D of radius 1 centered at the origin. For any source point s ∈ R 2 \D we then have Green's function for the Dirichlet problem given by with the reflected point y * = |y| −2 y ∈ D. Hence, for a source point s ∈ R 2 \D and corresponding s * = |s| −2 s ∈ D, we simultaneously have together with the boundary condition Therefore the boundary data produced by a set of external source points s j ∈ R 2 \D can be emulated using a corresponding set of internal source points s * j ∈ D. Hence, it suffices to solve the internal inverse source problem (see [1,4,10]) to recover s * j ∈ D and then by a simple transformation obtain the associated outer sources s j ∈ R 2 \D. This simplification is no longer possible for the Helmholtz equation. Nevertheless, it should be noted that this procedure may be used to produce an initial guess for the source locations when the wave number is small enough (see [14]).

Iterative solution
In the unpublished master thesis [14], supervised by one of us, an iterative solution method is proposed based on the reciprocity gap principle. Applying Green's theorem in view of the boundary conditions (2.3) and (2.7) it follows that the total field u satisfies for all solutions v ∈ H 1 loc (R 2 \D) to the Helmholtz equation satisfying the radiation condition. The method presented in [14] applies the reciprocity gap (3.1) for v = (x, ·) with source locations x on some auxiliary closed curve surrounding ∂D to obtain a set of nonlinear equations for the source locations s j and the intensities c j .
We suggest to slightly modify this approach and use as test functions v in (3.1) the fundamental solution with source points on the known boundary ∂D instead of source points away from the boundary. As to be expected this improves on the accuracy of the source reconstructions. Putting the source points on the boundary leads to Green's formula which can be seen as a particular case of Huygen's principle, that is, of the representation of the total field in theorem 3.12 in [2] by taking the trace on the boundary ∂D. We use (3.2), or more precisely a discretized version of (3.2), to solve the inverse source problem iteratively by linearizing simultaneously with respect to the intensities and the locations. Given a current approximation c 1 , c 2 , . . . , c n for the intensities and s 1 , s 2 , . . . , s n for the locations, we determine updates c 1 + γ 1 , c 2 + γ 2 , . . . , c n + γ n and s 1 + ξ 1 , s 2 + ξ 2 , . . . , s n + ξ n from the linearized equation This equation is linear with respect to the γ j and ξ j . It needs to be solved in a least squares sense after collocating it at a sufficient number of collocation points x m ∈ ∂D, m = 1, . . . , M.
For stability, a penalty term on the unknown locations and intensities has to be added, that is, Tikhonov regularization has to be employed to the solution of (3.3). Alternatively, one could reformulate the nonlinear problem (3.2) directly as a nonlinear least squares problem and use the Levenberg-Marquardt algorithm, requiring basically the same derivatives. For the integral on the right-hand side, i.e., the single-layer potential with density g the logarithmic quadrature formulae as described in section 3.5 of [2] can be employed. Of course, it would be desirable to have a general convergence result available and a criterion on the accuracy of the initial guess ensuring convergence. However convergence of Newton-type iterations for nonlinear ill-posed problems, in general, and for inverse scattering problems, in particular, is a very deep and challenging issue. So far it has not been clarified whether the general results on regularized Newton iterations on the solution of ill-posed nonlinear equations in a Hilbert space setting (among others see [3,5]) are applicable to inverse scattering problems. Hence, since we only want to present an initial analysis, we refrain from pursuing the convergence issue.
Instead of further investigations on convergence, we concentrate our attention on the more immediate issue of how to actually obtain a reasonable initial guess. We note that in order to start the iterations only an initial guess is required for the locations. From this, the initial guess for the intensities can be obtained via solving (3.2) in a least squares sense. Note that (3.2) is linear with respect to the intensities. An algorithm for providing an initial guess for the source locations is the subject of the following section.
Before we proceed we wish to point out that (3.2), in principle, can also be used for the case of incomplete data, i.e., for the case when g is known only on an open subset ⊂ ∂D. In this case, we just consider the missing part g| ∂D\ as an additional unknown and split the integral over ∂D accordingly. Of course, it is to be expected that this will effect the degree of ill-posedness.

Initial guess via a range test
Let B be a simply connected domain with C 2 boundary with outward normal ν such that D ⊂ B. We will design a compact linear operator A : L 2 ( ) → L 2 (∂D) such that the range of A depends on whether the source locations are contained in the annulus G := B\D or not.
To this end, we denote by N : L 2 (∂D) → H 1 (∂D) the Neumann-to-Dirichlet operator for the exterior domain R 2 \D subject to the Sommerfeld radiation condition. Further we denote by u 0 the solution to the exterior Neumann problem with boundary condition ∂u 0 ∂ν = g on ∂D. (4.1) In particular, we have u 0 | ∂D = Ng. Further we introduce the compact operators V , W : L 2 ( ) → L 2 (∂D) by the single-layer potential and its normal derivative If there exist densities ϕ ∈ L 2 ( ) and ψ ∈ L 2 (∂D) such that the total field u can be represented in the form then u is regular in G and therefore none of the source points s j , j = 1, . . . , n, lies within G.
Conversely, if none of the source points is contained in G, then u − u 0 is regular in G and can be represented as a single-layer potential with density in L 2 (∂G) provided k 2 is not a Dirichlet eigenvalue of the negative Laplacian neither for G nor for D. The latter assumption ensures bijectivity of the single-layer potential operator from and Then, in view of u = 0 and ∂ ν u = ∂ ν u 0 on ∂D by Holmgren's theorem the representation (4.2), that is, the composition u = u 0 + v + w in G, is equivalent to With the aid of the Neumann-to-Dirichlet operator N, that is, we can eliminate w from these two equations and arrive at the ill-posed equation Conversely, for a solution ϕ ∈ L 2 ( ) we define v by (4.3). Then we can represent the solution w of the exterior Neumann problem with boundary condition (4.5) as a single-layer potential with density ψ ∈ L 2 (∂) since the single-layer operator from L 2 (∂D) to H 1 (∂D) is bijective as a consequence of our assumption that k 2 is not a Dirichlet eigenvalue for D.
Summarizing we have proven the following theorem. We note that instead of (4.6) we could also use the equivalent equation with the operator where the inverse N −1 is given by the Dirichlet-to-Neumann operator.
The idea now is to check the solvability of (4.6), for example, via Tikhonov regularization. If A is injective and has dense range, the regularized solution of (4.6) converges as the regularization parameter tends to zero if and only if (4.6) is solvable (see section 4 in [2]). Numerically we can test for solvability by performing the Tikhonov regularization for a couple of reasonably small regularization parameters. If the regularized solution remains bounded while the regularization parameter is decreased then the equation is considered as solvable. Otherwise it is considered as unsolvable.
Here, the assumptions on injectivity and dense range are essential to ensure that the Tikhonov regularization scheme converges to the correct solution rather than to the solution given by the Moore-Penrose inverse. In this sense, we provide the following two theorems.  Proof. It suffices to show that the adjoint operator A * : L 2 (∂D) → L 2 ( ) of A is injective. Clearly, we have with the adjoints V * , W * : L 2 (∂D) → L 2 ( ) of V and W , and the adjoint N * : The adjoints of the potential operators are given by To characterize the adjoint N * of the Neumann-to-Dirichelet operator N, given any ψ 1 , ψ 2 ∈ L 2 (∂D), we define w 1 and w 2 as the radiating solutions to the Helmholtz equation in R 2 \D satisfying the Neumann boundary conditions ∂ ν w 1 = ψ 1 and ∂ ν w 2 = ψ 2 on ∂D, respectively. By Green's theorem we conclude Consequently, N * is given by For an explicit form of N, we introduce the operators S : L 2 (∂D) → H 1 (∂D) and Observing that k is not an interior eigenvalue for D, a single-layer representation of the solution to the exterior Neumann problem leads to the decomposition for the Neumann-to-Dirichelet operator (cf [2]). Here, I stands for the identity operator. We further introduce the operator T : H 1 (∂D) → L 2 (∂D) as the normal derivative of the double-layer potential by and note the identity (see section 3 in [2] or section 8 in [8]). From (4.11) and (4.12) we conclude that where we made use of (4.10). From this we observe that A * ψ = 0 implies v = w on . By the uniqueness for the exterior Dirichlet problem and the analyticity of v and w in R 2 \D this implies that v = w in R 2 \D, whence ∂v ∂ν = ∂w ∂ν on ∂D (4.14) in the sense of the jump relations in the L 2 sense (see section 3 in [2]). The latter also imply that ∂v ∂ν = (−I + K )ψ where we made use of (4.13). In view of the last two equations, from (4.14) we now can conclude that ψ = 0 and the proof is complete. Now theorem 4.1 suggests the following procedure to find the position of the unknown sources. For simplicity assume that ∂D is close to a circle. Extensions to other shapes are obvious. We choose for a sequence n of concentric circles with radius r n := r 0 + δn for n = 0, 1, 2, . . . with some appropriately selected increment δ > 0 and r 0 such that 0 is close to ∂D. We increase the radius until equation (4.6) becomes unsolvable for some radius r n in the sense that the L 2 norm of two different Tikhonov solutions to (4.6) for a couple of regularization parameters changes drastically. Then we know that one or more sources must be located in the annulus between the two circles n−1 and n . Therefore we modify the previous circle n−1 intending it locally through an exterior bump of diameter roughly equal to δ and rotate the modified curve n−1 equidistantly (see figure 1). If for some of these rotated bumpy contours n−1 equation (4.6) is unsolvable we choose a circle containing the corresponding bump. We then consider equation (4.6) for the union of all these circles. If the equation is solvable, the method for finding an initial guess stops and we choose the center of the circles as initial guesses. Otherwise, we go back to increasing the concentric circles, starting with n , but now we exclude from G the disks around the sources that we found already. For example, if we have, let say, m sources located on a ray extending from the center of D (and if they are not spaced too close to each other) we can determine their location by alternating m times between increasing the radius of the circles and rotating the bumpy curve.
Note that in this sense, in principle, our approach can be interpreted as a range test algorithm for finding the sources. In a postprocessing procedure, the accuracy of the reconstruction of the source locations could be improved by testing the solvability of (4.6) with the circles individually decreased or shifted around.

Numerical examples
In our examples we confine ourselves to the case where D is a disk of radius 1 centered at the origin and to the wave number k = 1. For the numerical discretization of the right-hand side in (3.2) and of the Neumann-to-Dirichlet operator (4.11) we used the logarithmic quadrature rules in [2] with 80 quadrature points over the boundary ∂D. The synthetic data were generated by using the combined single-and double-layer approach for solving the direct problem. Note that no inverse crime can be committed since the algorithm for the inverse problem does not contain the integral equations used for the solution of the direct problem.
We will present numerical results both for the range test method for finding an initial guess and for the iterative method.

Range test method for initial guess
First we describe the numerical implementation of the range test method as suggested at the end of section 4. We start by choosing as a concentric circle with radius 2 and then increase the radius by 1 until the equation becomes unsolvable. To test solvability, we solve equation (4.6) by Tikhonov regularization with regularization parameters α 1 = 10 −8 and α 2 = 10 −12 . We consider equation (4.6) to be unsolvable if the relative L 2 error between the two solutions for the different regularization parameters is larger than 50%. For the discretization of the integrals over we used the trapezoidal rule with 120 quadrature points.
After we obtain a circle for which unsolvability occurs we consider bumpy curves parametrized by where r is the largest radius of the concentric circle for which equation (4.6) was solvable, 2n is the number of bumpy curves used for fixed radius r and θ is an angle describing the location of the bump. The number n = round 5r 4 seems to be a good empirical choice along with 2n equidistant angles θ j = πj n , j = 1, 2, . . . , 2n, in order to cover most of the annulus between the circles of radii r and r + 1. For each θ i for which equation (4.6) was not solvable we picked the circle with radius 5 4 centered at r + 3 10 (cos θ i , sin θ i ) to exclude the corresponding source. 5.1.1. First example: one source. As a first example we considered a single source located at s 1 = (4, 0) with intensity 1, that is, the incident field is given by u i (x) = (x, (4,0)). We applied the range test algorithm with the configuration as described above and the result is plotted in figure 2. The disk that was obtained indeed covers the unknown source and the center of the circle can be used as an initial guess for the iterative method.
Applying the range test algorithm, again the disks cover the unknown sources as presented in figure 3.
However, we should note some restrictions of our approach. The circles considered to exclude the sources do not cover all the plane and therefore the method (as applied here) might not work in some cases. A more complete coverage of the plane by disjoint sets with well-behaved boundaries would then be needed and we leave this for future research. We also might have the situation where two or more sources lie within a circle, that is, our method does not actually separate all the sources. However, this could be taken into account by taking several initial guesses inside each circle for the iterative method. We also point out that our method has some difficulties in finding sources away from the obstacle, which is due to the ill-posed nature of the problem as pointed out above. Finally, we admit that as it is obvious from its construction the range test algorithm is very sensitive to noisy data.

Iterative method
We now illustrate the application of the iterative method as presented in section 3 using as initial guess the approximations for the source locations from the previous subsection. As a stopping criterion for the iterations we used an adaptation of Mozorov's discrepancy principle, for the mth iteration with the noise level n l . In this way, the regularization parameter decreases with the number of iterations at a slower rate for higher noise level.

First example: one source.
Again we considered a single source located at (4, 0) with intensity 1. We first worked with exact data using as initial guess the center of the circle obtained in section 5.1.1. The result as presented in figure 4 is extremely accurate with an error of order 10 −8 for the source location. We also applied the method with 1% and 5% noise on the data, using as initial guess two sources located ats 1 = (1, 1) ands 2 = (−1, −1), both with intensity 1, that is, with an incorrect number of sources in the initial guess. The results in figures 5 and 6 show that the sources 1 approximates the existing source, while the intensity of the second source tends to zero. The errors are of the same order as the noise level, as to be expected from Mozorov's principle.

Second example: three sources.
We go back to the example in subsection 5.1.2, that is, we want to recover the locations and intensities of the sources generating the incident field u i (x) = (x, (4, 0)) + 3 (x, (−3, 1)) − 2 (x, (−2, 4)).    Again we start by applying the method with exact data using as initial guess the centers of the circles obtained in subsection 5.1.2. The result is presented in figure 7, an extremely accurate reconstruction is exhibited again. Analogously to the previous example we applied the method with 1% and 5% noise on the data, using as initial guess four sources located at  all with intensity 1. In a similar way, the results in figures 8 and 9 show that three of the sources approximate the existing sources, while the intensity of the remaining source is close to zero.

Conclusion
We have proposed and analyzed a combination of a range test algorithm and an iterative scheme based on the reciprocity gap principle for the inverse source problem. Numerical examples provide evidence for the feasibility of this approach. Further research is required to make the method more efficient and to extend it to three dimensions. In principle, the method can also be applied to other boundary conditions than the sound-soft scatterer, to interior problems and to other differential equations.