A Multilevel Sampling Algorithm for Locating Inhomogeneous Media
aa r X i v : . [ m a t h . NA ] D ec A Multilevel Sampling Algorithmfor Locating Inhomogeneous Media
Keji Liu ∗ Jun Zou † July 8, 2018
Abstract
In the reconstruction process of unknown multiple scattering objects in inversemedium scattering problems, the first important step is to effectively locate someapproximate domains that contain all inhomogeneous media. Without such aneffective step, one may have to take a much larger computational domain than ac-tually needed in the reconstruction of all scattering objects, thus resulting in a hugeadditional computational efforts. In this work we propose a simple and efficientmultilevel reconstruction algorithm to help locate an accurate position and shapeof each inhomogeneous medium. Then other existing effective but computation-ally more demanding reconstruction algorithms may be applied in these initiallylocated computational domains to achieve more accurate shapes of the scatter andthe contrast values over each medium domain. The new algorithm exhibits severalstrengths: robustness against noise, requiring less incidences, fast convergence, flex-ibility to deal with scatterers of special shapes, and advantages in computationalcomplexity.
Key Words . Inverse medium scattering, initialization, multilevel reconstruction.
MSC classifications . 35R30, 65N20, 78A46.
This works investigates the numerical identification of inhomogeneous medium scat-terers by scattered fields. The inverse scattering problem can find wide applicationsin medicine, geophysics, biological studies. A large variety of numerical reconstruction ∗ Department of Mathematics, The Chinese University of Hong Kong, Shatin, Hong Kong.( [email protected] ) † Department of Mathematics, The Chinese University of Hong Kong, Shatin, Hong Kong. The workof this author was substantially supported by Hong Kong RGC grants (Projects 405110 and 404611).( [email protected] ) Consider an inverse scattering problem where the scatterer Ω, possibly consisting ofseveral separated disjoint components, is located in a homogeneous background medium R d ( d = 2 , u incj ( x ), j = 1 , , · · · , N i . For each plane waveincidence, the scattered field u sca ( x sq ) is measured by the receivers at locations x s , · · · , x sN s ; see Figure 1 for the incidences and receivers located on a circle S. χ =0 χ≠ χ≠ D receiversS u sca plane−waveincidences χ =0u inc Figure 1:
Geometrical model of the scattering problem.
The inverse scattering problem is to determine the contrast function or index of refraction, χ ( x ) for any point x varying in the scatterer Ω, given a set of N i N s scattering data u sca .The contrast χ has a very important property, namely it vanishes outside the scatteringobstacles. For each incident field, the total field u j satisfies the Helmholtz equation [6]:∆ u j ( x ) + k ( χ ( x ) + 1) u j ( x ) = 0 , x ∈ R d , (2.1)where k is the wavenumber of the homogeneous background media. The total field u j canbe represented by the intergral equation [6], u j ( x ) = u incj ( x ) + k Z Ω g ( x , x ′ ) χ ( x ′ ) u j ( x ′ ) dv ( x ′ ) , (2.2)where g ( x , x ′ ) is the Green’s function of the homogeneous background medium: g ( x , x ′ ) = ( i H (1)0 ( k | x − x ′ | ) for d = 2 , e ik | x − x ′| π | x − x ′ | for d = 3where H (1)0 is the zero-order Hankel function of first kind. We note that the total field u may stand for the acoustic pressure in an acoustic scattering problem, or for the electric3eld vector in an electromagnetic scattering, or for the particle-velocity vector in anelastodynamic scattering. The scattered field is measured on the boundary of a domainS, which is sitting outside the scatterer Ω. As the scatterer Ω is unknown, we shallintroduce a sampling domain D that completely cover the scatterer Ω. As the contrastfunction χ vanishes outside Ω, with the help of (2.2) we can write the scattered field as u scaj ( x ) = u j ( x ) − u incj ( x ) = k Z D g ( x , x ′ ) χ ( x ′ ) u j ( x ′ ) dv ( x ′ ) , x ∈ S . (2.3)For sake of convenience, we shall often introduce the contrast source function w j ( x ) = χ ( x ) u j ( x ) , x ∈ D . (2.4)Then we can write (2.2) and (2.3) in the following more compact forms w j ( x ) = χ ( x ) u incj ( x ) + χ ( x )( G D w j )( x ) , x ∈ D , (2.5)and u scaj ( x ) = ( G S w j )( x ) , x ∈ S , (2.6)where G D and G S are two integral operators given by( G D w )( x ) = k Z D g ( x , x ′ ) w ( x ′ ) dv ( x ′ ) ∀ x ∈ D , ( G S w )( x ) = k Z D g ( x , x ′ ) w ( x ′ ) dv ( x ′ ) ∀ x ∈ S . Equations (2.5) and (2.6) will be two fundamental equations for our proposed multilevelinitialization algorithm.
We can easily see that the support of the contrast source function w = χu describesthe exact locations and geometries of all the inhomogeneous media, which generate ascattered field u sca . The aim of this work is to propose a fast and less expensive algo-rithm that can help locate all the inhomogeneous media and provide good initial guessesfor some computationally more demanding iterative algorithms to find more accurateapproximations of the contrast function χ .Our algorithm will rely on the approximate contrast source obtained by backpropa-gation. Backprogation is widely used in inverse medium scatterings, see [17] [10] and thereferences therein. In this section we shall give a rigorous mathematical explanation ofthe approximate contrast source by backpropagation. Let ( · , · ) L ( S ) and ( · , · ) L ( D ) be thescalar products respectively in L ( S ) and L ( D ), and G ∗ s : L ( S ) → L ( D ) be the adjointof operator G s : L ( D ) → L ( S ). G ∗ s is called the backpropagation operator and given by( G ∗ S w )( x ) = k Z S g ( x , x ′ ) w ( x ′ ) ds ( x ′ ) ∀ x ∈ D .
4e shall need the following backpropagation subspace of L ( D ), V b = span { G ∗ s u sca } , which is formed by all the fields generated by the backpropagation G ∗ s on the scattereddata u sca . It follows from (2.6) that u sca ( x ) = ( G S w )( x ) , x ∈ S . (3.1)The backpropagation is to seek a best approximate solution w b to the equation (3.1) inthe backpropagation subspace V b , namely || u sca − G s w b || L ( S ) = min v b ∈ V b || u sca − G s v b || L ( S ) . (3.2)It is easy to see that the solution w b to (3.2) solves the variational system:( u sca − G s w b , G s v b ) L ( S ) = 0 ∀ v b ∈ V b , (3.3)or equivalently, ( G s w b , G s v b ) L ( S ) = ( G ∗ s u sca , v b ) L ( D ) ∀ v b ∈ V b . (3.4)As w b , v b ∈ V b , we can write w b = λ G ∗ s u sca , v b = µ G ∗ s u sca (3.5)for some constants λ, µ . Substituting the two expressions into (3.4) we obtain λ = || G ∗ s u sca || L ( D ) || G s G ∗ s u sca || L ( S ) , (3.6)which gives the approximate contrast source by backpropagation: w b = || G ∗ s u sca || L ( D ) || G s G ∗ s u sca || L ( S ) G ∗ s u sca . (3.7) In this section we propose a fast multilevel algorithm to find the locations and geome-tries of all the inhomogeneous media, which are described by the contrast function χ in(2.1). The algorithm proceeds iteratively, and carries out two important steps at eachiteration based on the two fundamental equations (2.5) and (2.6), namely the state andfield equations. In the first step, we apply the backpropagation technique to compute anapproximate contrast source w j corresponding to each incident u incj ( j = 1 , , · · · , N i ). Itfollows from (3.7) that this approximation is given by w j = || G ∗ S u scaj || L ( D ) || G S G ∗ S u scaj || L ( S ) G ∗ S u scaj , j = 1 , , · · · , N i . (4.1)5ith these approximate contributions w j of the exact contrast source w corresponding toeach incident u incj , we approximate the contrast χ pointwise by minimizing the residualequation corresponding to the state equation (2.5), namely min χ ( x ) ∈ R N i X j =1 (cid:12)(cid:12)(cid:12)(cid:16) χu incj − w j + χG D w j (cid:17) ( x ) (cid:12)(cid:12)(cid:12) , (4.2) which yields an explicit formula to compute an approximate contrast value χ ( x ) at everypoint x ∈ D when an approximate contrast source w j is available: χ ( x ) = P N i j =1 w j ( x ) (cid:0) u incj + G D w j )( x ) P N i j =1 | ( u incj + G D w j )( x ) | , (4.3)where the overbear denotes the complex conjugate.Clearly both (4.1) and (4.3) are rather crude mostly, and may provide rather poorapproximations for the exact contrast source w and contrast profile χ . But as it willbe seen, when we combine these two poor approximations in a novel manner with somemultilevel technique, it generates a very efficient and robust algorithm for locating anaccurate position and shape of each inhomogeneous medium.Our basic idea is also based on another simple observation. We know that the exactcontrast function χ ( x ) vanishes outside the scatterer Ω, so its support provides the lo-cation and shape of the scatterer Ω, which is formed by all the inhomogeneous media.This observation, along with the previous two explicit evaluation formulae (4.1) and (4.3)and a novel multilevel technique, forms the foundation of our new multilevel samplingalgorithm.For the description of the algorithm, we first introduce two new concepts, the smallestdistance and the first gap interval with index M . For a given finite positive sequence, { a , a , · · · , a m } , its smallest distance is the positive smallest one among all the distancesbetween two neighboring elements, namely dist ( a i , a i +1 ), i = 1 , , · · · , m −
1. Among allthese m − j such that 2 ≤ j ≤ m − dist ( a j , a j +1 ) is M times larger than the smallest distance of the sequence { a , a , · · · , a j } ,then [ a j , a j +1 ] is called a gap interval . The first such interval is called the first gap interval .Now we are ready to state our new algorithm. Multilevel Sampling Algorithm.
1. Choose a sampling domain D that contains the scatterer Ω;Select a uniform (coarse) mesh on D , consisting of rectangular (2D) or cubic (3D)elements; write the mesh as D ;Select a tolerance ε and an index M ; set an initial cut-off value c := 0 and k := 1.2. Compute an approximate value of the contrast χ k ( x ) at each grid point x ∈ D k ,using the formulae (4.1) and (4.3). Then do the following:6.1 Order all the values of χ k ( x ) satisfying χ k ( x ) ≥ c k − into a non-decreasingsequence;Find the first gap interval of the sequence with index M ;Choose the right endpoint of this first gap interval as the next cut-off value c k .2.2 If χ k ( x ) ≥ c k at a grid point x , select all the grid points of the elements whichshare x as one of its vertices;Remove all the grid points in D k , which are not selected;Update D k by all those selected grid points.3. If | c k − c k − | ≤ ε , go to Step 4;otherwise refine the mesh D k to get D k +1 ; set k := k + 1 and go to Step 2.4. Output all grid points in D k for the supports of all inhomogeneous media;Output the approximate contrast values χ k ( x ) at all grid points of D k .We can easily see that the Multilevel Sampling Algorithm does not involve any opti-mization process or matrix inversions, and its major cost is to update the contrast valuesusing the explicit formulae (4.1) and (4.3) at each iteration, and the computational sam-pling domain D k shrinks as the iteration goes. So the algorithm is rather simple andless expensive. In addition, as the cut-off values are basically to distinguish the homo-geneous background medium where χ ( x ) vanishes and the inhomogeneous media where χ ( x ) should be essentially different from 0, so our cut-off values are rather easy to chooseand insensitive to the size and physical features of scatterers. In fact, the cut-off valuecan start simply with zero, then it is updated automatically with the iteration. As weshall see from numerical examples in the next section, the algorithm works well with fewincidents, even with one; and it is self-adaptive, namely it can recover some elements thathave been removed at the previous iterations due to the computational errors. In termsof these aspects the new Multilevel Sampling Algorithm outperforms the popular linearsampling methods [15], including the improved multilevel variant [11]. Discretization . We end this section with a brief discussion about numerical dis-cretization of the integrations involved in the above multilevel algorithm. We illustrateonly the discretization of (2.5) and (2.6), as all other integrations involved in the algorithmcan be approximated similarly. To do so, we divide the domain D into smaller rectangularor cubic elements, whose centers are denoted as x , x , . . . , x M . Using the coupled-dipolemethod (CDM) or discrete dipole approximation (DDA) [3, 9], we can discretize (2.5) by w j ( x m ) = χ ( x m ) u incj ( x m ) + k χ ( x m ) X n = m A n g ( x m , x n ) w j ( x n ) , m = 1 , , . . . , M, (4.4)where A n is the area or volume of the n -th element. Similarly, we can discretize equation(2.6) at every point x ∈ S by u scaj ( x ) = k M X n =1 A n g ( x , x n ) w j ( x n ) for j = 1 , , · · · , N i . (4.5)7 Numerical simulations
In this section we present several examples to test the effectiveness of our multilevelsampling algorithm. We first list the parameters that are used in our numerical simu-lations: the wave number k and wave length λ are taken to be k = 2 π and λ = 1; thenumber of incidences and receivers are set to be N i = 6 and N s = 30 respectively, bothequally distributed on the circle of radius 5 λ ; the index M of the first gap interval andthe tolerance parameter ε are chosen to be 100 and 10 − respectively. In all the numericalsimulations, random noises are added to the exact scattering data in the following form: u scaj ( x ) := u scaj ( x )[1 + ξ ( r ,j ( x ) + ir ,j ( x ))] , j = 1 , , · · · , N i where r ,j ( x ) and r ,j ( x ) are two random numbers varying between -1 and 1, and ξ corresponds to the level of the noise, which is usually taken to be 10% unless specifiedotherwise. All the programs in our experiments are written in MATLAB and run on a2.83GHz PC with 4GB memory. Example 1 . This example shows a scatterer Ω consisting of two squares of side length0 . λ , located respectively at ( − . λ, − . λ ) and (0 . λ, . λ ), with their contrast valuesbeing 1 and 2 respectively; see the two red squares in Figure 2(a). We take the samplingdomain D = [ − . λ, . λ ] × [ − . λ, . λ ], which is quite large compared to the scattererΩ, with its area of 64 times of the area of one scatterer component. More importantly, wesee that these two small objects are quite close to each other. The mesh refinement duringthe multilevel algorithm is carried out based on the rule h k = 0 . λ/ k , where k is the k -th refinement, and h and h k are respectively the mesh sizes of the initial mesh and themesh after the k th refinement. The numerical reconstructions are shown in Figure 2(b) -2(d) respectively for the 1st, 3rd and 5th iteration. One can observe from the figuresthat the algorithm converges very fast and provides accurate locations of the two mediumcomponents in only 5 iterations. Moreover, we can see an important advantage of thealgorithm, i.e., it can separate the disjoint medium components quickly. Example 2 . This example is the same as Example 1, except that the contrast valuesof two medium components are now variable functions, namely χ ( x, y ) = sin π (10 | x | − . sin π (10 | y | − . . The numerical reconstructions are shown in Figure 3(b)- 3(d) for the 1st, 3rd and 6thiteration. Again, we observe from the figures that the algorithm converges fast, providesvery satisfactory locations of the two medium components in only 6 iterations, and it canseparate the disjoint medium components quickly.
Example 3 . This example considers a scatterer Ω of a thin annulus with the innerand outer radii being 0 . λ and 0 . λ respectively and centered at the origin. The contrast8 D (a) (b) D D (c) (d)Figure 2: (a) The initial (coarse) mesh on the sampling domain; (b)-(d) Reconstructionsat the 1st, 3rd and 5th iteration. value χ ( x ) is 2 inside the thin annulus. The sampling domain D is taken to be a squareof side length with 5 . λ , as shown in Figure 4(a).It is easy to see the sampling domain D has an area of about 62 times as large asthe one of the annulus, and the annulus has a very thin thickness, i.e., 0 . λ . The meshrefinement during the multilevel algorithm is carried out based on the rule h k = 0 . λ/ k ,where k indicates the k -th refinement and h k is the mesh size after the k th refinement. Thenumerical reconstructions are shown in Figure 4(b)- 4(d) for the 1st, 3rd and 4th iteration.Same as for the previous two examples, the reconstructions are quite satisfactory and theaccurate locations of the scatterer can be achieved. χ There are many numerical reconstruction methods available in the literature for thecontrast profile function χ , that are robust and more accurate than our new multilevelmethod. But these methods are usually more complicated technically and much more9 D (a) (b) D D (c) (d)Figure 3: (a) The first coarse mesh on the sampling region; (b)-(d) Reconstructions at the1st, 3rd and 6th iteration. demanding computationally, as they mostly involve nonlinear optimizations and matrixinversions. Without a reasonably good initial location for each inhomogeneous medium,we may have to take a much larger sampling domain than required so these methods canbe extremely time consuming, especially in three dimensions. Using the newly proposedmultilevel algorithm in Section 4, we can first locate a much smaller sampling domainthan usual (or the one we originally selected) in a numerical reconstruction for the con-trast χ . Then we can apply any existing reconstruction algorithms for more accuratereconstructions, starting with an initial sampling domain provided by the new multilevelalgorithm. This may save us a great fraction of the entire computational costs.Next we show numerical tests by combining our new multilevel method with thepopular extended contrast source inversion (ECSI) method [17]. Firstly, we considerthe same scatterer Ω and the set-ups as in Example 1 of Section 5.1; see Figure 5(a).Then we apply the ECSI method [17] with mesh size h = 0 . λ to the reconstructeddomain (cf. Figure 2(d)) by the multilevel algorithm. The reconstruction result is shownin Figure 5(b). Clearly, one can see that both the location and the values of the contrast10 D (a) (b) D D (c) (d)Figure 4: (a) The first coarse mesh on the sampling region; (b)-(d) Reconstructions atthe 1st, 3rd and 4th iteration. are well reconstructed.Secondly, we consider the same scatterer Ω and the set-ups as in Example 3 of Section5.1; see Figure (5)(c). We apply the ECSI method [17] with mesh size h = 0 . λ tothe reconstructed domain (cf. Figure 4(d)) obtained by the multilevel algorithm. Thereconstruction result is shown in Figure 5(d). Clearly, both of the location and the valuesof the contrast are well reconstructed. Example 4 . This example tests a three-dimensional scatterer Ω consisting of twosmall cubic components:Ω = [ − . λ, − . λ ] , Ω = [0 . λ, . λ ] . The two squares are quite close to each other, both with constant contrast values 2; seeFigure 6(a). We take the sampling domain to be D = [ − . λ, . λ ] , which is about 500times of the volume of Ω or Ω . 11 −1 −0.5 0 0.5 1−1−0.8−0.6−0.4−0.200.20.40.60.81 00.20.40.60.811.21.41.61.82 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−1−0.8−0.6−0.4−0.200.20.40.60.81 00.20.40.60.811.21.41.61.82 (a) (b) −2 −1 0 1 2−2.5−2−1.5−1−0.500.511.522.5 00.20.40.60.811.21.41.61.82 −2 −1 0 1 2−2.5−2−1.5−1−0.500.511.522.5 00.20.40.60.811.21.41.61.82 (c) (d)Figure 5: (a) True scatterer in Example 1; (b) Reconstruction by combining the multilevelmethod with ECSI for Example 1; (c) True scatterer in Example 3; (d) Reconstruction bycombining the multilevel method with ECSI for Example 3 We take an initial mesh size of h = 0 . λ in the multilevel algorithm. The mesh refine-ment during the multilevel algorithm is carried out based on the rule h k = 0 . λ/ k , where k is the k -th refinement, and h k is the mesh size after the k th refinement. The numericalreconstructions are shown in Figure 7. Same as for the previous two-dimensional exam-ples, the reconstructions are quite satisfactory and the accurate locations of the scatterercan be achieved, and two inhomogeneous medium objects can be quickly separated. Example 5 . In this test we consider a torus scatterer (see Figure 6(b)), with acontrast value 2. The torus has the following representation, (cid:18) R − p x + y (cid:19) + z = r , where r = 0 . λ and R = 0 . λ ( R is the radius from the center of the hole to the centerof the torus tube, r is the radius of the tube). The sampling doman is taken to be D = [ − . λ, . λ ] , which is about 170 times of the volume of the torus.12 D (a) (b)Figure 6: Scatterers imbedded in a large sampling domain: two cubic components closeto each other (a); a torus (b) D D (a) (b)
D D (c) (d)Figure 7:
Numerical reconstructions by the first 4 iterations D (a) (b) D D (c) (d)Figure 8:
Numerical reconstructions by the first 3 iterations; (d) is the same as (c), butviewed in a different angle.
We take an initial mesh size of h = 0 . λ in the multilevel algorithm. The meshrefinement during the multilevel algorithm is carried out based on the rule h k = 0 . λ/ k ,where k is the k -th refinement, and h k is the mesh size after the k th refinement. Thenumerical reconstructions are shown in Figure 8. Same as for the previous two-dimensionalexamples, the reconstructions are quite satisfactory and the accurate locations of thescatterer can be achieved. This work proposes a multilevel sampling algorithm which helps locate an initial com-putational domain for the numerical reconstruction of inhomogeneous media in inversemedium scatterings. The algorithm is an iterative process which starts with a large sam-pling domain, and reduces the size of the domain iteratively based on the cut-off values,which are computed adaptively by using the updated contrast source strengths and con-trast values at each iteration. The iterative algorithm can be viewed actually as a direct14ethod, since it involves only matrix-vector operations and does not need any optimiza-tion process or to solve any large-scale ill-posed linear systems. The algorithm works withvery few incident fields and its cut-off values are easy to compute and insensitive to thesizes and shapes of the scatterers, as well as the noise in the data. This is a clear ad-vantage of the algorithm over some popular existing sampling methods such as the linearsampling type methods, where the cut-off values are sensitive to the noise and difficultto choose, and the number of incident fields can not be small. In addition, the multilevelalgorithm converges fast and can easily separate multiple disjoint scattering components,often with just a few iterations to find a satisfactory initial location of each object. An-other nice feature of the new algorithm is that it is self-adaptive, that is, it can remedythe possible errors from the previous levels at each current level. With an effective initiallocation of each object, we may then apply any existing efficient but computationallymore demanding methods for further refinement of the estimated shape of each scatteringobject as well as for recovery of the contrast profiles of different media.
References [1] A. Abubakar and P. M. van den Berg, The contrast source inversion method forlocation and shape reconstructions, Inverse Problems (2002), pp. 495-510.[2] G. Bao and P. Li, Inverse medium scattering for the Helmholtz equation at fixedfrequency, Inverse Problems (2005), pp.1621-1641.[3] K. Belkebir, P. C. Chaumet and A. Sentenac, Superresolution in total internal re-flection tomography, J. Opt. Soc. Amer. A (2005), pp. 1889-1897.[4] X. Chen, Application of signal-subspace and optimization methods in reconstructingextended scatterers, J. Opt. Soc. Amer. A (2009), pp. 1022-1026.[5] X. Chen, Subspace-based optimization method for solving inverse-scattering prob-lems, IEEE Trans. Geosci. Remote Sensing (2010), pp. 42-49.[6] D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory,2nd ed., Springer Verlag, Berlin, 1998.[7] K. Ito, B. Jin and J. Zou, A direct sampling method to an inverse medium scatteringproblem, Inverse Problems (2012), 025003 (11pp).[8] A. Kirsch, The MUSIC-algorithm and the factorization method in inverse scatteringtheory for inhomogeneous media, Inverse Problems (2002), pp. 1025-1040.[9] A. Lakhtakia, Strong and weak forms of the method of momoents and the coupleddipole method for scattering of time-harmonic electromagnetics fields, Int. J. ModernPhys. C (1992), pp. 583-603. 1510] B. Levy and C. Esmersoy, Variable background Born inversion by wavefield back-propagation, SIAM J. Appl. Math. (1988), pp. 952-972.[11] J. Li, H. Liu and J. Zou, Multilevel linear sampling method for inverse scatteringproblems, SIAM J. Sci. Comp. (2008), pp. 1228-1250.[12] J. Li, H. Liu and J. Zou, Strengthened linear sampling method with a reference ball,SIAM J. Sci. Comp. (2009), pp. 4013-4040.[13] K. Liu, Y. Xu and J. Zou, A parallel radial bisection algorithm for inverse scatteringproblems, Inv. Prob. Sci. Eng. (2012), pp. 1-13.[14] E. A. Marengo, F. K. Gruber and F. Simonetti, Time-reversal MUSIC imaging ofextended targets, IEEE Trans. Image Proc. (2007), pp. 1967-1984.[15] R. Potthast, A survey on sampling and probe methods for inverse problems, InverseProblems (2006), pp. R1-R47.[16] P. M. van den Berg and R. E. Kleinman, A contrast source inversion method, InverseProblems (1997), pp. 1607-1620.[17] P. M. van den Berg, A. L. van Broekhoven and A. Abubakar, Extended contrastsource inversion, Inverse Problems15