Compressive Inverse Scattering II. Multi-Shot SISO Measurements with Born Scatterers
Inverse scattering methods capable of compressive imaging are proposed and analyzed. The methods employ randomly and repeatedly (multiple-shot) the single-input-single-output (SISO) measurements in which the probe frequencies, the incident and the sampling directions are related in a precise way and are capable of recovering exactly scatterers of sufficiently low sparsity.
For point targets, various sampling techniques are proposed to transform the scattering matrix into the random Fourier matrix. Two schemes are particularly interesting: The first one employs multiple frequencies with the sampling angle always in the back-scattering direction resembling the synthetic aperture (SA) imaging; the second employs only single frequency with the sampling angle in the (nearly) forward scattering direction in the high frequency limit, resembling the setting of X-ray tomography.
The results for point targets are then extended to the case of localized extended targets by interpolating from grid points. In particular, an explicit error bound is derived for the piece-wise constant interpolation which is shown to be a practical way of discretizing localized extended targets and enabling the compressed sensing techniques.
For distributed extended targets, the Littlewood-Paley basis is used in analysis. A specially designed sampling scheme then transforms the scattering matrix into a block-diagonal matrix with each block being the random Fourier matrix corresponding to one of the multiple dyadic scales of the extended target. In other words by the Littlewood-Paley basis and the proposed sampling scheme the different dyadic scales of the target are decoupled and therefore can be reconstructed scale-by-scale by the proposed method. Moreover, with probes of any single frequency the coefficients in the Littlewood-Paley expansion for scales up to can be exactly recovered.
Consider the scattering of the incident plane wave
by the variable refractive index where is the incident direction. The scattered field satisfies the Lippmann-Schwinger equation 
where is the Green function of the operator . We assume that the wave speed is unity and hence the frequency equals the wavenumber .
The scattered field has the far-field asymptotic
where is the scattering amplitude. In inverse scattering theory, the scattering amplitude is the measurement data determined by the formula 
The main objective of inverse scattering then is to reconstruct the medium inhomogeneities from the knowledge of the scattering amplitude. In Part I  and this paper the target to be imaged consists of a finite number of point scatterers. And the main techniques for reconstruction are from theory of compressed sensing. In  we analyze the one, but high, frequency imaging method with the single-input-multiple-output (SIMO) and multiple-input-multiple-output (MIMO) measurements in which for every incident plane wave the scattering amplitude is sampled at multiple directions independent of the incident wave.
In this paper the focus is on the multi-shot single-input-single-output (SISO) measurement in which for every randomly selected incident plane wave the scattering amplitude is sampled at only one direction correlated with the incident wave.
Our motivation for this alternative imaging method is practical as well as theoretical. On the theoretical aspect, the analysis of the high-frequency SIMO/MIMO method employs the coherence theory of compressed sensing which deals with only random targets under a suitable sparsity constraint. On the other hand, the multi-shot SISO method proposed in this paper is amenable to the restricted isometry theory of the random Fourier matrix which guarantees reconstruction for all targets under the weakest known sparsity constraint. On the practical aspect, the present method can achieve a comparable performance with a much lower frequency or bandwidth (Figure 7). Moreover, the case of extended targets can be treated by either interpolating from grid points or using the wavelet basis in this approach. The main drawback, though, of the present approach, in comparison to that of , is that the multiple scattering effect is not accounted for.
In Sections 2, we discuss the case of point scatterers and propose several sampling schemes to transform the scattering matrix into the random Fourier matrix which is amenable to the compressed sensing techniques. One scheme employs multiple frequencies with the sampling angle always in the back-scattering direction resembling the synthetic aperture (SA) imaging; another scheme employs only single frequency with the sampling angle in the (nearly) forward scattering direction in the high frequency limit, resembling the setting of X-ray tomography. We then extend these results to the case of localized extended targets by interpolating from grid points in Section 3. In Section 4 we analyze the case of distributed extended targets using the Littlewood-Paley basis and propose a sampling scheme to block-diagolize the scattering matrix. Each block is in the form of random Fourier matrix and corresponds to one dyadic scale of the target. Hence our method has the capability of imaging the target scale-by-scale by the compressed sensing techniques. Moreover, the coefficients in the Littlewood-Paley expansion for scales up to can be exactly recovered by using probes of any single frequency . We numerically test these sampling methods and compare their success probabilities in Section 5. We conclude and comment on the issue of resolution in Section 6.
2. Point scatterers
For the simplicity of notation, we will focus on two dimensions below.
We consider the medium with point scatterers located in a square lattice
of spacing . The total number of grid points in is a perfect square. Without loss of generality, assume where and . Let be the strength of the scatterers. Let be the locations of the scatterers. Hence .
For the discrete medium the scattering amplitude becomes a finite sum
Unlike  which covers both linear and nonlinear scattering, here we work exclusively under the Born approximation in which the exciting field is replaced by the incident field ; unlike  which deals exclusively with one frequency, here we will work with multiple frequencies.
Let be various incident and sampling directions for the frequencies to be determined later. Define the measurement vector with
The measurement vector is related to the target vector by the sensing matrix as . Let be the polar angles of , respectively. The -entry of is
Note the all entries of have unit modulus. As in  we reconstruct as the solution of the -minimization, called Basis Pursuit (BP):
In the presence of noise/error of size as in
(8) is replaced by the relaxation scheme
A fundamental notion in compressed sensing under which BP yields the unique exact solution is the restrictive isometry property (RIP) due to Candès and Tao . Precisely, let the sparsity of a vector be the number of nonzero components of and define the restricted isometry constant to be the smallest positive number such that the inequality
holds for all of sparsity at most and some constant . For the target vector let denote the best -sparse approximation of in the sense of -norm, i.e.
where denotes the number of nonzero components, called the sparsity, of . Clearly, consists of the largest components of .
For general , one can consider the normalized version of (9)
by which it follows that
We wish to write the -entry of the sensing matrix in the form
where are independently and uniformly distributed in in view of the following theorem.
To construct a sensing matrix of the form (14) we proceed as follows. Write in the polar coordinates as
where is a parameter to be determined later (30). Equivalently we have
This set of equations determines the single-input--single-output- mode of sampling.
and set the frequency to be
Then the entries (7) of the sensing matrix have the form
By the square-symmetry of the problem, it is clear that the relation (19) can be generalized to
On the other hand, the symmetry of the square lattice should not play a significant role and hence we expect the result to be insensitive to any fixed , independent of , as long as (20) holds. Indeed this is confirmed by numerical simulations below (Figure 5).
Let us focus on three specific measurement schemes.
Scheme I. This scheme employs band limited probes, i.e. . This and (20) lead to the constraint:
. In this case the scattering amplitude is always sampled in the back-scattering direction. This resembles the synthetic aperture imaging which has been previously analyzed under the paraxial approximation in . In contrast, the forward scattering direction with almost surely violates the constraint (23).
Scheme II. This scheme employs single frequency probes no less than :
with . The difference between the incident angle and the sampling angle is
which diminishes as . In other words, in the high frequency limit, the sampling angle approaches the incident angle. This resembles the setting of the X-ray tomography.
Scheme III. This scheme employs probes of unlimited frequency band. Let be arbitrary distinct numbers in and let and be determined by (22) and (20), respectively. The possibility of having a small divisor in (20) renders the bandwidth unlimited in principle.
Let be independently and uniformly distributed in and let be the polar coordinates of , i.e.
3. Localized extended targets
In this section we consider extended targets that are localized in space. We represent such targets by interpolating from grid points and thus extend the preceding results for point targets to localized extended targets. We will treat the case of distributed extended targets that are not localized in space in the next section.
Suppose that the target function is continuous and has a compact support. Consider the filtered version of
where for any integrable filter function such that . Clearly, tends to as .
Next we discretize (31) by replacing the integral by the Riemann sum of step size . We obtain
which is as smooth as the interpolation element is.
Since has a compact support, is a finite set. For simplicity let be the square sublattice
of total cardinality . Let . Define the target vector with . Let and be the probe frequencies and directions, respectively, and let be the sampling directions for . Now we write the data vector in the form (9) with the sensing matrix elements
and the error term due to the filtered discretization. For sufficiently small we may assume for a given (see Remark 1 and below).
The crucial observation is that the sensing matrix (33) is identical to (7) with and any isotropic filter function . Therefore Theorem 1 holds verbatim for the case of localized extended targets formulated above.
How small must and be in order to ensure that ? This can be answered roughly as follows. First, by the inequality it suffices to have . Now consider the transformation , defined by
from the space of continuous functions supported on to , cf. (4). By definition
and we have
For the sampling schemes I and II (with ) of Section 2, we can give an explicit bound for the discretization error
where denotes the -norm of functions defined on .
For example, consider the smooth isotropic filter function
with With the choice , the condition (34) can be formulated as
For another example, consider the indicator function on the unit square . The resulting interpolation is the piece-wise constant approximation. Then
Setting again we obtain the condition
which appears to be much more useful than (35).
4. Distributed extended targets
In this section, we consider extended targets represented by square-integrable functions .
To this end we use the Littlewood-Paley basis
Then the following set of functions
forms an orthonormal wavelet basis in . In terms of the Littlewood-Paley basis we have the expansion
Define the sensing matrix elements to be
and let , where are given below.
for some , be the column and row indexes, respectively, of . It is important to keep in mind how and are related to and , respectively, in order to understand the scheme described below.
Let be independent, uniform random variables on and define
Suppose for all . This holds true, for example, when the frequencies satisfy the constraint
Let be the polar coordinates of .
As before, let be the angles of incidence and sampling, respectively, which are chosen according to
By (48) and the definition of it is clear that are zero if . Consequently the sensing matrix is the block-diagonal matrix with each block (indexed by ) in the form of random Fourier matrix
be the target vector. Let
We can then write the measurement vector where . The above observation means that the target structures of different dyadic scales are decoupled and can be determined separately by our approach using compressed sensing techniques.
by (50). On the other hand if then
Using the RIP for the random Fourier matrix of each block, we obtain the following theorems analogous to Theorem 1.
For each , let (15) be satisfied for and and any . Suppose
The condition (56) implies that the wavelengths are no larger than the scales under interrogation, consistent with the classical resolution criterion.
Theorem 3 allows reconstruction with probes of single sufficiently high frequency
The beauty of the theorem, however, lies in the fact that with probes of any frequency the coefficients in the Littlewood-Paley expansion for scales up to can be recovered (exactly in the absence of noise).
5. Numerical results
Greedy algorithms have significantly lower computational complexity than linear programming and have provable performance under various conditions. For example under the condition the Subspace Pursuit (SP) algorithm is guaranteed to exactly recover via a finite number of iterations . We have used SP for reconstruction in all our simulations with the following parameters: . The probability of recovery is calculated by using 1000 independent runs. In Figures 3-7, the vertical axis is for the probability of recovery and the horizontal axis is for the number of point scatterers.
To test Scheme I numerically, we use (25) and
to see if the deviations from (24) have any impact on performance. Their probabilities of recovery are plotted as a function of the sparsity in Figure 3. Clearly, the performance deteriorates rapidly as the difference between the sampling and incident angles decreases. This is due to increasingly more frequent and more severe violation of (23) as a result. In other words, the backward scattering direction is the optimal sampling direction for Scheme I.
Likewise, to test Scheme II numerically, we use (26),
The results for case (a) at are shown in Figure 4 and the results for case (b) with at are displayed in Figure 5. The slight degradation in performance for e.g., results from the violation of (29) (and hence (23)) which affects the performance significantly for small . For large , it has been shown in  that the sampling angles can be chosen independently of the incident angles to maintain a high level of performance (Figure 7). On the other hand, when (29) holds, the performance is essentially independent of both and , Figure 5.
For Scheme III we use equally spaced incident angles , (20) and
As shown in Figure 6, the performance is essentially independent of .
We demonstrate numerically the high-frequency SIMO schemes analyzed in  and compare their performance with that of the multi-shot SISO schemes presented above. A case in point would be to use (26) and (28) but for a fixed incident angle, say (In this case (29) is almost certainly violated). In  it is established that the SIMO schemes with sampling angles, independent of the incident angles, achieves a high performance in reconstructing a sparse target with a sufficiently high-frequency probe wave (i.e. ). The success probabilities of the SIMO schemes for are calculated and plotted in Figure 7. Consistent with the theory , the low frequency case with has the worst performance. Clearly the performance of the SIMO schemes improves with and in the limit approaches that of the multi-shot SISO Scheme II. There is negligible difference between the performances for and both of which follow closely that of the SISO Scheme II with (black dashed line in Figure 7).
Finally we demonstrate the reconstruction with the Littlewood-Paley basis. Since the scattering matrix is block-diagonalized by the proposed sampling scheme (43)-(44) and (46)-(47), we consider targets of single pair of dyadic scales. In this simulation, the reconstruction is carried out for a target on the scales () with nonzero coefficients () among possible modes (, modes in each coordinate) by using samples (). The target and its reconstruction in the domain are color-coded and displayed in Figure 8.
6. Conclusion and discussion
We have proposed, analyzed and numerically tested several multi-shot SISO sampling schemes which transform the scattering matrix into the random Fourier matrix in the case of point and localized extended scatterers and the block-diagonal form of random Fourier matrices in the case of distributed extended targets.
In the case of point scatterers, these sampling schemes are either multi-frequency band limited (I) or single frequency outside band (II). For Scheme I , the sampling direction is the backward direction analogous to synthetic aperture radar (SAR) while for Scheme II in the high frequency limit the optimal sampling direction is in the forward direction, analogous to the X-ray tomography. Both schemes produce nearly the same recovery probability with the resolution given by (30), i.e.
We have extended this approach to the case of localized extended targets by interpolating from grid points. We have formulated the approximate scheme as inversion with noisy data. In particular, we have derived an explicit error bound for the simple piece-wise constant interpolation.
In the case of distributed extended targets, the block-diagonal form of the scattering matrix in the Littlewood-Paley representation means that different dyadic scales of the target are decoupled and can be imaged scale-by-scale separately by our method. Moreover, we can determine the coefficients in the Littlewood-Paley expansion (40) for scales up to by using probes of any single frequency . The disadvantage of the Littlewood-Paley basis is that a localized target has slowly decaying coefficients and hence is not compressible in this basis.
The SIMO schemes in which the scattered field of an incident angle is measured at multiple sampling angles have been studied in [12, 13, 14]. Except for the special case of the periodic scatterers lying on a transverse plane , it is not known if the sensing matrices of the SIMO schemes in general satisfy the RIP or not. In this case, the approach based on the notion of incoherence is taken to analyze the SIMO schemes [12, 14]. This approach is generally more flexible and should be applicable to the SISO schemes considered here.
The main advantage of the SIMO schemes is that in the one-shot setting (one incident field) the inverse scattering problem can be solved exactly without the Born approximation by inverting an auxiliary nonlinear system of equations . We are working to extend the idea to the imaging methods with multi-shot measurements.
On the other hand, the SISO schemes with the RIP tend to have a better performance which can be matched by that of the SIMO schemes only at high frequency as demonstrated in Figure 7. In the high frequency regime the sampling angles can be chosen independently of the incident angles .
- M. Born and E. Wolf, Principles of Optics, 7-th edition, Cambridge University Press, 1999.
- A.M. Bruckstein, D.L. Donoho and M. Elad, “From sparse solutions of systems of equations to sparse modeling of signals,” SIAM Rev. 51 (2009), 34-81.
- E. J. Candès, “The restricted isometry property and its implications for compressed sensing,” Compte Rendus de l’Academie des Sciences, Paris, Serie I. 346 (2008) 589-592.
- E. J. Candès, J. Romberg and T. Tao, “Robust undertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inform. Theory 52 (2006), 489-509.
- E.J. Candès, J. Romberg and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Commun. Pure Appl. Math. 59 (2006), 1207Ð23.
- E. J. Candès and T. Tao, “ Decoding by linear programming,” IEEE Trans. Inform. Theory 51 (2005), 4203Ð4215.
- E. J. Candès and T. Tao, “ Near-optimal signal recovery from random projections: universal encoding strategies?,” IEEE Trans. Inform. Theory 52 (2006), 54-6-5425.
- S.S. Chen, D.L. Donoho and M.A. Saunders, “Atomic decomposition by basis pursuit,” SIAM Rev. 43 (2001), 129-159.
- D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory. 2nd edition, Springer, 1998.
- W. Dai and O. Milenkovic, “Subspace pursuit for compressive sensing: closing the gap between performance and complexity,” arXiv:0803.0811.
- I. Daubechies, Ten Lectures on Wavelets. SIAM, Philadelphia, 1992.
- A. Fannjiang, “Compressive inverse scattering I. High frequency SIMO measurements, ” arXiv: 0906.5405
- A. Fannjiang, “Compressive imaging of subwavelength structures,” SIAM J. Imag. Sci. 2 (2009), 1277-1291.
- A. Fannjiang, P. Yan and Thomas Strohmer, “Compressed remote sensing of sparse objects,” arXiv:0904.3994
- H. Rauhut, “Stability results for random sampling of sparse trigonometric polynomials,” preprint, 2008.
- M. Rudelson and R. Vershynin, “On sparse reconstruction from Fourier and Gaussian measurements,” Comm. Pure Appl. Math. 111 (2008) 1025-1045.
- V. N. Temlyakov, “Greedy approximations,” Acta Num. (2008), pp. 235Ð409.