An acoustic Lippmann-Schwinger inversion method: applications and comparison with the linear sampling method

We introduce an imaging method based on solving the Lippmann-Schwinger equation of acoustic scattering theory. We compare and contrast the proposed Lippmann-Schwinger inversion with the well-established linear sampling method using numerical examples. We demonstrate that the two imaging methods are physically grounded in different but related wave propagation problems: Lippmann-Schwinger inversion seeks to reconstruct the space and time dependence of a scatterer based on the observed scattered field in a performed physical experiment, whereas the linear sampling method seeks to focus wave fields in a simulated virtual experiment by estimating the space and time dependence of an inverse source function that cancels the effects of the scatterer at a specified focusing point. In both cases, the medium in which the waves propagate is the same; however, neither method requires prior knowledge or assumptions on the physical properties of the unknown scatterer—only knowledge of the background medium is needed. We demonstrate that the linear sampling method is preferable to Lippmann-Schwinger inversion for target-oriented imaging applications, as Lippmann-Schwinger inversion gives nonphysical results when the chosen imaging domain does not contain the scatterer.


Introduction
Imaging is a qualitative inverse scattering problem that seeks to localize the boundaries between media of different physical properties. In this paper, we restrict our attention to imaging acoustic media, in which scalar pressure fields propagate and the scattering is caused by variations in the bulk modulus and mass density of the media. Imaging the boundaries of such variations is challenging due to imprecise (or altogether absent) knowledge of the medium. If the acoustic impedance contrast is sufficiently strong, much of the wave energy that is sent into a medium can become localized amongst the inhomogeneities, where it is scattered and redirected multiple times, before finally being recorded at a receiver. Such strong multiple scattering makes it difficult to interpret where the observed scattered fields originated.
To construct accurate images, it is necessary for imaging algorithms to properly account for the multiple scattering undertaken by a wave field as it propagates through a medium. Of fundamental importance to scattering theory is the Lippmann-Scwhinger equation (e.g. [1]), which explains not only primary (or single) scattered waves, but all multiply scattered waves as well. The Lippmann-Schwinger equation provides an exact representation of the scattered field in terms of a weighted superposition of the impulse response of the background medium over the region containing the scatterer. The weights of the superposition are determined by the physical properties of the scatterer (e.g. the variations in acoustic impedance), and also on the total field incident to the scatterer. Naturally, the Lippmann-Schwinger equation forms the basis of many quantitative inversion schemes that seek to recover some of the physical properties of a scatterer (e.g. [2][3][4]).
The Lippmann-Schwinger equation has been investigated for imaging purposes as well. However, previous attempts at using the Lippmann-Schwinger equation for imaging have often relied on a weak-scattering approximation (e.g. [5]) or a re-normalization of the Lippmann-Schwinger equation to obtain an absolutely convergent series which can then be inverted (e.g. [6,7]). In this paper, we introduce a noniterative technique for inverting the Lippmann-Schwinger equation without making any weak-scattering approximations. The obtained solution approximates the spatiotemporal dependence of the scatterer needed to reconstruct the observed scattered field.
The linear sampling method is another imaging technique that has garnered considerable attention [8][9][10][11][12][13]. To demarcate the boundary of a scatterer, the method relies on a characteristic blowup behavior of the solution to an ill-posed integral equation. The method works by attempting to estimate the spatiotemporal dependence of an inverse source that cancels the effects of the scatterer at a given point in the medium. Such a process has a physical solution with finite energy only when the given point lies inside the scatterer. Clearly, the inverse source cannot exist for any point outside the scatterer, and numerically this 'law of causality' manifests as a blowup in the norm of the solution [13]. Consequently, an image of the scatterer can be obtained by noting where the norm of the solution becomes arbitrarily large. More recently, it has been shown that quantitative information about the scatterer can be extracted from the linear sampling method [14][15][16][17][18].
The linear sampling method has been interpreted as a focusing technique [13,19]. As we show in this paper, the solution to the ill-posed integral equation in the linear sampling method corresponds to a focusing function, which is a type of distributed source for focusing wave fields in a medium. When the specified point lies inside the scatterer, the focusing function contains the inverse time dependence of the scatterer needed to focus wave fields onto the specified point.
Our paper is organized as follows. We formulate the direct acoustic scattering problem to establish the relevant physical quantities of interest in section 2. Lippmann-Schwinger inversion is introduced in section 3 and compared with the linear sampling method in section 4. Numerical experiments are given in section 5 and the conclusions follow in section 6.

Formulation of the direct acoustic scattering problem
We are interested in the problem of reconstructing the shape of a scatterer from measurements of the scattered pressure field. We assume the scattering is due to localized inhomogeneities, which represent variations in the bulk modulus κ and mass density ρ of the acoustic medium. We assume the inhomogeneities have total compact support Ì  D 3 , where we denote the boundary of D by ∂D and the closure of D by Specifically, we assume the acoustic medium is determined by the functions where κ 0 , ρ 0 denote the background bulk modulus and mass density and κ s , ρ s denote perturbations in the bulk modulus and mass density inside the scatterer, respectively. In this paper, we further assume the density is continuous across the boundary of the scatterer, so that r  0 s as  ¶D x from inside the scatterer. Variations in the bulk modulus and mass density will in general cause the pressure field to propagate at different speeds, since the acoustic velocity varies as k r = c . For a point source located at Î  x s 3 , let p denote the total pressure field satisfying the acoustic wave equation where δ is the Dirac delta distribution and z Î  C 2 ( )is a time-dependent function that describes the shape of the wave. With the background velocity given by k r = c 0 0 0 , we define the refractive indexn, the velocity contrastm, and the density ratioq as Using the definitions given in (2), the total pressure field can be decomposed into an unperturbed wave p 0 and a scattered wave p s such that p=p 0 +p s is the unique solution to (1). It follows that if the unperturbed pressure field p 0 satisfies then the scattered pressure field p s satisfies where χ is the contrast source function given by Expression (5) states that the contrast source function χ is a combination of both monopole sources (due to velocity variations) and dipole sources (due to density variations). Note that the time dependence of the contrast source is determined by the total pressure field p at the scatterer. As seen from the definitions given in (2), the velocity contrast m and the gradient of the density ratio ∇q are in general nonzero only inside the scatterer. It follows that the contrast source function χ has compact support D. In deriving expression (5), we have assumed that the density is smoothly varying across the boundary of the scatterer. If the density is discontinuous across the boundary of the scatterer, the contrast source function has an additional term proportional to the jump discontinuity in the density [20]. We do not address such complexities in this paper.
Let G 0 denote the unperturbed Green function of the background medium satisfying  (8) is the Lippmann-Schwinger equation for the scattered pressure field p s . If the contrast source (5) is known, we could use equation (8) to exactly compute the scattered pressure field observed at any point x and at any time t. In practice, however, the contrast source function (5) is unknown since we do not know the velocity or density variations of the scatterer, nor the total pressure field inside the medium. On the other hand, suppose we know (or can estimate) the background medium without the scatterer so that we can compute the unperturbed Green function G 0 . Then, provided we can separate the scattered pressure field p s from the recorded total field p, equation (8) provides a means to solve for the unknown contrast source function χ without requiring any prior knowledge or weak-scattering approximations of the contrast source.

Lippmann-Schwinger inversion
In an imaging experiment, we record the total pressure field p at receiver locations x r , which are typically restricted to an acquisition surface we denote by Γ r . Similarly, the sources used to generate the pressure field are restricted to points x s of an acquisition surface we denote by Γ s (possibly equal to Γ r ). Both surfaces Γ r and Γ s are assumed to be disjoint from D . By evaluating equation (8) on the receiver surface Γ r , we can relate the recorded scattered field to the unknown contrast source function.
The inversion of equation (8) for the contrast source function is well known to be ill-posed due to nonuniqueness of the solution (e.g. [21][22][23]). This ill-posedness holds in general for the inverse source and inverse medium problems [22]. As a simple example, consider a wave incident to a planar interface at normal incidence. Certain combinations of the mass density and velocity exist such that the contrast in acoustic impedance is zero, resulting in a reflectionless interface. The inverse problem of determining the interface from the (identically zero) scattered field does not have a unique solution: we cannot distinguish the trivial solution (i.e. no interface) from any of the specific combinations of velocity and density that yield zero impedance contrast. As another example, Zhang et al [24] construct an 'invisible' medium, where single-reflected waves from deeper layers of the medium exactly cancel multiply reflected waves from shallower layers. Since the reflected waves vanish, one cannot reconstruct the contrast across the interfaces from the data. Moreover, the data we collect in practice is discrete and finite, whereas the contrast source is a continuous function of space and time with infinitely many degrees of freedom. Hence, the inversion of equation (8) is also ill-conditioned. Consequently, a regularization scheme is needed to obtain an approximate solution that satisfies a set of desired constraints (in section 5, we use Tikhonov regularization [25] to solve the inverse source problem). Additionally, the unperturbed Green function satisfying equation (6) may not be a known function, since the background medium may be arbitrarily heterogeneous. In such cases, numerical methods are needed to approximate the Green function. Furthermore, to properly account for the spatial distribution of the contrast source function χ, it is necessary to discretize the volume integral in equation (8) using a quadrature rule. This may prove challenging and computationally expensive when attempting to accommodate the arbitrary structure of the background medium.
To address these issues, let Î  S C 2 ( )be a time-dependent function possessing the same frequency band as the scattered pressure field p s . For an arbitrary point Î  z 3 , we define the test function which represents a band-limited impulse response of the background medium due to a monopole point source at location z, and can be effectively modeled using numerical methods. The choice of the time function S need only be made in consideration of the frequency band of the data, and does not necessarily need to be equal to the generating pulse function ζ (which may be unknown). Next, let Ω denote a subset of the background medium containing the scatterer, so that Ì W D . Our imaging procedure is then based on approximating the spatiotemporal dependence of the contrast source function χ over the domain Ω. We do this by regularizing the solution of the modified Lippmann-Schwinger equation in a least-squares sense. By using expression (9) as a proxy for the unperturbed Green function G 0 in equation (8), we can avoid the cost of performing a quadrature rule for the volume integral in equation (10) by discretizing the domain Ω into a sufficiently dense grid, say a few grid points per the dominant wavelength of the test function Ψ.
Such an approach effectively exploits the fact that a superposition of the band-limited test functions will achieve sufficient smoothness in the obtained solution.
Due to the compact support of the contrast source function, we expect the solution χ to be nonzero only for those points zäD. It follows that an image of the scatterer can be obtained by plotting the energy of the solution to equation (10) over the domain Ω and noting where the energy is greater than zero: We numerically implement our proposed imaging technique as follows. Given a model of the background medium, we discretize the imaging domain Ω into N z distinct points using a regularly sampled grid that covers the scatterer. Let N r denote the number of receivers, N t the number of time samples, and N s the number of sources. For each grid point z n , n=1, K, N z , we discretize the test function observed at the receivers as and Δt is the time sampling interval. For a general heterogeneous background medium, numerical methods are needed to compute these test functions. In section 5, we consider the special case of a constant background medium for which the test functions assume a simple, analytic expression. In any case, the test functions can be efficiently computed using source-receiver reciprocity whenever the number of grid points N z is expected to be greater than the number of receivers N r . Similarly, the scattered pressure field is discretized as Figure 1 shows a schematic in which we construct a regularly sampled grid that covers an unknown scatterer.
Next, we discretize the modified Lippmann-Schwinger equation (10) into a sum over the grid points z n . Since equation (10) is convolutional in time, the problem is efficiently solved in the frequency domain. Here, we use a circumflexˆto denote frequency-domain quantities and w to denote a vector of angular frequencies. For each source j=1, K, N s , we take the scattered pressure field j p , , (· · ) and solve for a regularized solution to the discretized Lippmann-Schwinger equation in a least-squares sense. An image is obtained for each source as where the 2-norm is computed over all angular frequencies w. The final image produced from Lippmann-Schwinger inversion is defined as the root-mean-square of the normalized images I j : It follows that the final image satisfies 0I LSI (z n )1 for all points z n in the imaging domain Ω, where values close to 1 indicate points which are likely inside the scatterer and values close to 0 indicate otherwise.

Comparison with the linear sampling method
The linear sampling method exploits the use of the source position x s ä Γ s as a free parameter in the physical experiment to design wave fields of a simulated virtual experiment [14]. By multiplying equations (3) and (4) by a function j Î G´ L s 2 ( )and integrating over G´ s , it follows by linearity that the functions are also solutions to equations (3) and (4), respectively. The functions v j and w j are the 'unperturbed' and 'scattered' fields of a virtual experiment corresponding to the unperturbed and scattered fields p 0 and p s of the physical experiment, respectively. The dependence of the functions v j and w j on the choice of j is emphasized in the notation.
In the linear sampling method, we force the scattered field w j of the virtual experiment to coincide with a band-limited impulse response Ψ z that radiates from a fixed point z in the imaging domain Ω, as observed on the receiver surface Γ r . This amounts to evaluating expression (14) on the receiver surface Γ r to form the ill-posed near-field equation Since the scattered field of the virtual experiment is made to radiate from a single point z ä Ω, the incident field j v z of the virtual experiment must focus onto this point. To construct this focusing wave field, we regularize the solution of the near-field equation (15) for the function j z in a least-squares sense. We call the solution j z a focusing function, which is parameterized by the focusing point z ä Ω. The focusing function j z extracts the properties of the contrast source χ at the prescribed point z ä Ω from the information encoded in the scattered field p s . Of course, this is physically possible only when the prescribed point lies inside the scatterer. It can be shown that the focusing function j is proportional to the inverse of the contrast source function χ [13]. In particular, the inverse of χ exists and is bounded in the L 2 norm only for points zäD (where χ is strictly nonzero). Otherwise, for points Î  D z 3 ⧹ , the inverse of χ does not exist and the solution j z becomes unbounded (due to a division by zero since χ is zero outside the scatterer). This is the so-called blowup behavior of solutions to the near-field equation.
Note that the impulse response appearing on the right-hand side of equation (15) is that due to a monopole point source, whereas the contrast source function (5) for a more general acoustic scatterer is comprised of both monopole and dipole sources. As demonstrated by Haddar et al [12], the linear sampling method requires a combination of both monopole and dipole test functions to accurately reconstruct obstacles of mixed monopole and dipole source distributions. Here, a dipole test function is defined as the gradient of expression (9) with respect to the variable x r , taken in the direction of a user-specified unit vector Î  d 2 , where  2 denotes a unit sphere embedded in  3 . Of course, the boundary of the scatterer is unknown, so the true orientation of a dipole source at a given point z ä Ω is not known a priori. Consequently, a search procedure must be undertaken to find an optimal direction d * for the dipole test function [12]. The computational cost for such a procedure can be significant in all but the simplest cases (e.g. a scatterer with known spherical symmetry [26]), as it requires the solution of equation (15) for each direction Î  d 2 of the dipole test function at each point z ä Ω. In this paper, we apply the linear sampling method for monopole test functions only.
In the frequency domain, the discretization of the near-field equation (15) leads to the linear system

Due to the blowup behavior of the solution ĵ , it is convenient to define an image as
As was the case for Lippmann-Schwinger inversion, it follows that the image obtained from the linear sampling method satisfies 0I LSM (z n )1 for all points z n in the imaging domain. By comparing the discretized equations (11) and (16), we see that Lippmann-Scwhinger inversion is reciprocal to the linear sampling method in the sense that the roles of the scattered pressure fields p and the test functions Y have been interchanged. In equation (11), we take a linear combination of the test functions over the imaging domain to solve for the contrast source function that best reconstructs the observed scattered field at the receivers. In equation (16), we take a linear combination of the observed scattered fields to solve for the inverse source function that best reconstructs a prescribed impulse response. Thus, the matrix operations in equations (11) and (16) perform the same actions but using opposite kernels. From a computational perspective, this means that a generic convolutional operator can be defined once, and the test functions Y or the scattered pressure fields p can be substituted as the kernel of the operator to switch back and forth between Lippmann-Schwinger inversion and the linear sampling method, respectively.
An important difference between Lippmann-Schwinger inversion and the linear sampling method lies in the domain of definition of their solutions. The contrast source function χ obtained from Lippmann-Schwinger inversion is a global solution defined over the entire imaging domain Ω. Consequently, the image obtained from Lippmann-Schwinger inversion depends on every point in the imaging domain. Conversely, the focusing function j z obtained from the linear sampling method is defined over the acquisition surface Γ s and parameterized by a single image point z ä Ω. Thus, each image point obtained from the linear sampling method depends only on the blowup behavior of the solution j z , and is independent of every other point in the imaging domain.
This difference between the two methods has profound implications should the chosen imaging domain not contain the scatterer. In particular, Lippmann-Schwinger inversion would fail to give a physically plausible source function χ, since the least-squares solution of equation (11) forces a nonzero solution over a given imaging domain. Consequently, the resulting image cannot always be trusted. This issue is of significant practical importance since the location of the scatterer is often not known in advance. In contrast, the image obtained from the linear sampling method would remain unaffected by the choice of the imaging domain, since the method relies on the blowup behavior of the solutions j z to localize the scattering target. We illustrate this important difference in the next section.

Numerical experiments
In what follows, we demonstrate Lippmann-Schwinger inversion as a source reconstruction method for the performed physical experiment and the linear sampling method as a focusing technique for the corresponding virtual experiment. Special attention is paid to the influence of the chosen imaging domain on the obtained images for each method. Finally, we compare and contrast the two methods in their ability to image morphologically random scatterers.  For simplicity, we model the scattering experiments in a two-dimensional acoustic medium. The background medium is taken to be homogeneous, characterized by a constant wave speed c 0 =2 m s −1 and a constant density ρ 0 =1 kg m −3 . In the following experiments, we interrogate the medium using 24 sources and receivers which are placed along a circle surrounding the scatterers. The time dependence of each source is given by the Ricker pulse where the dominant frequency ν will be specified before each example. The test function is computed for each point z n in the background medium as z when kΔtx i −z n /c 0 and is zero otherwise.
To generate the scattered pressure field p s at the receivers, we simulate the total pressure field p in the presence of the scatterers by finite-difference modeling and subtract from it the unperturbed pressure field computed using equation (18) with the grid points z n replaced by the source points x s .
We use the LSMR algorithm [27] and Tikhonov regularization [25] to solve equations (11) and (16) in a leastsquares sense. The value of the regularization parameter α is indicated in each example.

Demonstrating the physical bases of the methods
To test the physical interpretation of Lippmann-Schwinger inversion as a source reconstruction technique and the linear sampling as a focusing technique, we consider the velocity model shown in figure 2. The model is intentionally simple, consisting of a single box-shaped anomaly embedded within the homogeneous background medium. In this experiment, we interrogate the scatterer using the Ricker pulse with a dominant frequency ν=25 Hz. The dominant wavelength λ of the Ricker pulse is shown in figure 2 for scale. We apply Lippmann-Schwinger inversion and the linear sampling method to two different imaging domains: one domain that contains the scatterer and another domain that does not (Ω 1 and Ω 2 in figure 2, respectively). Both imaging domains are discretized into a 15×15 point regularly sampled grid.
The left column of figure 3 shows snapshots of the modeled scattered pressure field corresponding to the source shown by the red dot. The boundary of the scatterer is indicated by the white dashed line. The right column of figure 3 shows the same snapshots of the reconstructed scattered pressure field obtained from Lippmann-Schwinger inversion applied to the imaging domain Ω 1 (shown by the yellow dashed line). In this case, the obtained contrast source function χ correctly replicates the spatiotemporal behavior of the scatterer to reconstruct the scattered pressure field observed at the receivers. The anomalous waveforms seen in the reconstructed wave field are likely due to coarse receiver sampling along the acquisition circle and could be reduced with denser acquisition. Figure 4 shows the same snapshots of the modeled scattered pressure field but now alongside the reconstructed scattered pressure field obtained from Lippmann-Schwinger inversion applied to the imaging domain Ω 2 . In this case, the imaging domain does not contain the scatterer and the obtained contrast source function χ poorly reconstructs the scattered pressure field observed at the receivers. Interestingly, the contrast source still attempts to reconstruct the back-scattered waves observed at the receivers at the top-left quadrant of the acquisition circle. The imprint of the coarse receiver sampling is again obvious in the reconstructed waveforms. Figure 5 shows snapshots of the corresponding virtual experiments obtained from the linear sampling method. The boundary of the scatterer is again shown by the white dashed line. The left column shows the total wave field before, at, and after the time of focusing for the case when the image point (shown as the yellow dot) lies inside the scatterer. In this case, the obtained focusing function j correctly synchronizes the sources (shown as blue dots) to focus the incident field onto the prescribed image point. Conversely, the right column shows the same snapshots of the total wave field for the case when the image point does not lie inside the scatterer. Here, the obtained focusing function j is nonphysical, since it could not extract the properties of the contrast source χ at the given image point. Consequently, focusing is not achieved. Figure 6 shows the obtained images from Lippmann-Schwinger inversion (left column) and the linear sampling method (right column) using definitions (12) and (17), respectively. The top row shows the reconstructions for the imaging domain Ω 1 , which contains the scatterer. Clearly, both methods accurately reconstruct the shape of the box anomaly. The bottom row shows the reconstructions for the imaging domain Ω 2 , which does not contain the scatterer. Here, Lippmann-Schwinger inversion incorrectly suggests there is a scatterer inside the imaging domain, whereas the linear sampling method correctly indicates there is no scatterer. The images in the bottom row are plotted on the same scale as the corresponding images in the top row.

Imaging random inclusions
We compare Lippmann-Schwinger inversion and the linear sampling method in their ability to reconstruct morphologically random inclusions ( figure 7). In this experiment, the random inclusions represent perturbations in both the model wave speed and mass density, with c=2.3 m s −1 and ρ=1.3 kg m −3 inside the inclusions. As discussed in section 2, the jump discontinuities in the density across the boundaries of the random inclusions will contribute an additional term to the contrast source function which is not accounted for in Lippmann-Schwinger inversion. Similarly, because the impulse response appearing on the right-hand side of the near-field equation is due to a monopole source, the linear sampling method will fail to capture the directional dependence of the density discontinuities. Still, it is instructive to compare how well these methods perform when the physics of the scattering problem is not fully captured in the imaging algorithms.
We interrogate the acoustic medium using staggered sources and receivers along the acquisition circle as depicted in figure 7. Since the sources and receivers are not co-located, we use reciprocity to add 'fictitious' data to the linear system of equations and enhance our coverage of the imaging domain. In this experiment, we interrogate the scatterers using the Ricker pulse with a dominant frequency ν=20 Hz. The dominant wavelength λ of the Ricker pulse is shown in figure 7 for scale. We construct an imaging domain covering the region from 0.7 m to 1.3 m along the x-axis and from 0.2 m to 0.8 m along the y-axis with a 31×31 point regularly sampled grid.  We test the imaging methods under three different scenarios. In the first scenario, the background medium is known exactly and the scattered wave is measured noise-free. An example shot gather for this scenario is shown in figure 8(a). In the second scenario, we test the ability of the methods to reconstruct the random inclusions in the presence of noisy data and an erroneous background model. Here, we incorrectly assume the background velocity is c 0 =2.3 m s −1 , which is a 15% increase from the true background velocity. We add random noise to the scattered pressure field that is uniformly distributed within the same frequency band as the data. We set the signal-to-noise ratio (SNR) to 0.3, which is defined as where P denotes average power (i.e. the mean-square amplitude). Here, P signal denotes the average signal power taken over 24 traces per shot gather. A common shot gather of the scattered pressure field with noise is shown in figure 8(b).
In the third scenario, we generate the scattered pressure field p s using only the density model shown in right panel of figure 7, while the acoustic velocity is held fixed at c 0 =2 m s −1 everywhere in the model. Since neither Lippmann-Schwinger inversion nor the linear sampling method are formulated to account for the jump discontinuities in the density, both methods are expected to fail. Figure 8(c) shows a typical shot gather corresponding to this limiting case. Figure 9 shows the images obtained from Lippmann-Schwinger inversion (left column) and the linear sampling method (right column) for each of the three scenarios. The top row shows the images obtained for the case of the exact background model and noise-free data. Both methods clearly reconstruct the shapes of the random inclusions. This suggests that imaging with these methods is still feasible despite jump discontinuities in the density so long as the velocity contrast is the dominant term in the contrast source function (5). The middle row shows the images obtained for the case of an erroneous background model and noisy data. For both methods, we observe a dilation in the size and location of the random inclusions due to the erroneous background velocity. The imprint of the noisy signal is evident in the gritty texture of the reconstructions compared with the images in the top row. The bottom row shows the images obtained for the case in which the scatterers represent pure density contrasts. As expected, both methods have failed to identify any coherent structure of the random inclusions.

Discusion and Conclusions
We have presented an imaging method based on solving the Lippmann-Schwinger equation of acoustic scattering theory. Compared with the linear sampling method, both techniques formulate the imaging problem without using any prior knowledge or weak-scattering approximations of the unknown scatterer. We have validated the physical interpretations of Lippmann-Schwinger inversion as a source reconstruction method for the performed physical experiment and the linear sampling method as a focusing technique for the corresponding virtual experiments. Our numerical experiments show that the two imaging methods give comparable reconstructions when the imaging domain fully contains the scatterer. Even in the presence of noise and an erroneous background model, both imaging techniques give reasonable reconstructions. However, only the linear sampling method gives physically consistent results when the imaging domain does not contain the scatterer.
As discussed by Martin [20], many practical applications of acoustic scattering involve jump discontinuities in the mass density across the boundary of a scatterer. Unfortunately, these jump discontinuities manifest as an additional surface integral in the Lippmann-Schwinger equation over the unknown boundary of the scatterer. This poses a significant challenge for Lippmann-Schwinger inversion in reconstructing such jump discontinuities, especially if the scatterer is primarily caused by strong and sudden density variations. According to expression (5), the contrast source function is sensitive to the illumination direction provided by the known source position. It may be possible to glean more detailed information about the density properties of a scatterer by analyzing the gradients of the contrast source with respect to known source positions. As discussed by Haddar et al [12], the sensitivity of the linear sampling method to dipole source distributions (such as density variations) can be improved by modeling the impulse response due to a dipole source. How to optimally choose the orientation for the dipole test function remains an open problem, however. Additionally, the procedure described by Haddar et al [12] requires solving the near-field equation for monopole and dipole sources separately; it would be of practical interest if this could be done concurrently. Further work on acoustic inverse scattering methods should consider these suggestions for handling jump discontinuities in the mass density.