Synthetic aperture imaging is used in many applications such as ultrasonic non-destructive testing, mine detection, surveillance, and radar imaging. The main idea behind synthetic aperture imaging is that a single transmitter/receiver is used to probe an unknown region by emitting known pulses into the medium and recording the time-dependent responses as it moves along a given path. Fourier transforming these time-dependent measurements yields their corresponding frequency responses. In this work we focus our attention on the synthetic aperture radar imaging problem. However, the methodology used here can be directly applied to other related problems.
Several imaging methods have been proposed in the literature for imaging with SAR data. The traditional SAR image is formed by evaluating the data at each measurement location at the travel time that it takes for the waves to propagate from the platform location to a point in the imaging region on the ground and back. The resolution of this image increases with the synthetic aperture and the system bandwidth . When the phases of the waves are recorded with high accuracy, SAR imaging produces high-resolution images of the reflectivity on the ground. It is well known however that SAR imaging is quite sensitive to noise in the phase. Such noise may result from uncertainty in the platform motion and/or scattering by randomly inhomogeneous media. For SAR imaging with noise in the phase, we refer to  for the application of coherent interferometry (CINT) to SAR imaging and to the more recent work in  on a high-resolution interferometric method for imaging through scattering media.
for imaging effectively direction and frequency dependent reflectivities using the multiple measurement vector (MMV) framework. As in other applications, using sparsity-constrained optimization methods significantly increases the resolution of the SAR image. However, the computational cost of optimization is significantly higher than that for sampling methods such as SAR or CINT, which simply consist of evaluating an imaging functional at each grid point on a mesh of the imaging region.
MUSIC (multiple signal classication) is another sampling method that has been widely used in several imaging applications [8, 9, 12, 16]. To explain the main idea of MUSIC, let us consider the single-frequency array imaging problem. For this problem the data is a matrix, called the array response matrix. The -th element of the array response matrix corresponds to the data received at the -th array element when the
-th element is an emitter. The singular value decomposition (SVD) of the array response matrix is used to determine the signal and noise subspaces of the data. Next, a model for the illumination vectoris introduced, with denoting a point in the imaging region. The illumination vector is the vector of measurements received along the receiving array due to a source at . If a target is located at , then is in the signal subspace of the array response matrix. Thus, the projection of onto the noise subspace is zero or very small. In MUSIC, one forms an image by evaluating one over this norm of the projection of onto the noise subspace. The peaks appearing in the MUSIC image give the locations of the targets with high resolution. Although MUSIC effectively and efficiently produces high-resolution images, it does not apply directly to SAR imaging data.
In this paper we introduce a modification and generalization of MUSIC for SAR imaging. This imaging method modifies SAR data by using the Prony method  to rearrange frequency-dependent data at one measurement location as a matrix. Then we form a block-diagonal matrix with the set of Prony matrices from all spatial locations on the flight path. An image of the reflectivity on the ground is then formed using a signal subspace method applied to this block-diagonal matrix. This signal subspace method is a generalization of MUSIC that projects the illuminating vector for each point in the imaging region on both the noise and signal subspaces . The noise subspace provides high spatial resolution and the signal subspace provides quantitative information about the targets. The result of combining these two subspaces is a high-resolution quantitative imaging method. The relative balance between the noise and signal subspaces depends on the noise level in the data which is controlled through a user-defined regularization parameter, .
There are two main results in this paper. The first main result is the resolution analysis for this modified and generalized MUSIC method that shows an enhancement in resolution compared to classical SAR imaging by a factor . Namely we obtain a cross-range resolution of and a range resolution of . Here denotes the speed of the waves, denotes the bandwidth, denotes the synthetic array aperture, denotes the distance from the center of the flight path to the center of the imaging region, and denotes the range offset (see Fig. 1). The second main result is the stability analysis of the method to random perturbations of the travel times. This analysis shows that the method provides stable reconstructions when is chosen to satisfy with
denoting the maximum variance of the random perturbations of the travel times. Our numerical simulations are in agreement with these theoretical findings. Moreover, they show that the proposed method provides statistically stable results with signal-to-noise ratios comparable to CINT, but with much better resolution.
The remainder of this paper is as follows. In Section 2 we give a brief description of synthetic aperture radar imaging and define the measurements. In Section 3 we describe the Prony method that we use to rearrange the frequency data and show why it is appropropriate for signal subspace imaging. We define the two imaging functionals that we use for quantitative signal subspace imaging in Section 4. In Section 5 we give a resolution analysis for the imaging method. We consider this imaging method when the travel times have random perturbations in Section 6 and give results for the expected value and statistical stability of the image formed using this method. We show numerical results that support our theory in Section 7. Section 8 contains our conclusions.
2. Synthetic aperture radar imaging
In synthetic aperture radar (SAR) imaging, a single transmitter/receiver is used to collect the scattered electromagnetic field over a synthetic aperture that is created by a moving platform [6, 7, 14]. The moving platform is used to create a suite of experiments in which pulses are emitted and resulting echoes are recorded by the transmitter/receiver at several locations along the flight path. Let denote the broadband pulse emitted and let denote the data recorded. Here, the measurements depend on the slow time that parameterizes the flight path of the platform, , and the fast time in which the round-trip travel time between the platform and the imaging scene on the ground is measured. In SAR imaging, one seeks to recover the reflectivity of an imaging scene from these measurements.
Although SAR uses a single transmit/receive element, high-resolution images of the probed scene can be obtained because the data are coherently processed over a large synthetic aperture created by the moving platform. As illustrated in Fig. 1, the platform is moving along a trajectory probing the imaging scene by sending a pulse and collecting the corresponding echoes. We call range the direction that is obtained by projecting on the imaging plane the vector that connects the center of the imaging region to the central platform location. Cross-range is the direction that is orthogonal to the range. Denoting the size of the synthetic aperture by and the available bandwidth by , the typical resolution of the imaging system is in range and in cross-range. Here is the speed of light and the wavelength corresponding to the central frequency while denotes the distance between the platform and the imaging region and the offset in range.
We use the start-stop approximation, which is typically done in SAR imaging. This approximation assumes that the change in displacement between the targets and the platform is negligibly small compared to the travel time it takes for the pulse emitted to propagate to the imaging scene and return as echoes. This approximation is valid in radar since the speed of light is orders of magnitude larger than the speed of the targets and the platform.
Using this start-stop approximation, we the consider the measurements only at discrete values of , corresponding to for . Next, we suppose that is digitally sampled at values of . Consequently, these data have a discrete Fourier transform denoted by evaluated at frequencies denoted by for . This choice of samples is to make the notation in Section 3 simpler. With these assumptions, we find that our measurement data is given by the matrix whose columns are
3. Rearranging frequency data
The data matrix is not suitable for direct application of signal subspace methods. Therefore, we introduce a rearrangement of the data based on the Prony method  which, for the -th column of , yields the following matrix,
In this rearrangement, the first column is the truncation of to its first entries. Subsequent columns are sequential upward shifts of truncated to its first entries.
To see why this rearrangement is suitable for signal subspace imaging, consider the Born approximation for a single point target. Let denote the reflectivity of the point target and denote its position. According to the Born approximation, the scattered field at frequency is the spherical wave,
with denoting the wave speed and denoting the field incident on the point target. Suppose that the signal emitted at position on the flight path at frequency is a spherical wave with unit amplitude. For that case, the measurement corresponding to the scattered field evaluated at is given by
It follows that
from which we find that
Next, suppose that the frequencies are sampled according to for with a fixed constant. For that case, we can rewrite (6) as with ,
Here, we have written the reflectivity as and included in and arbitrarily in (7).
Suppose there are non-interacting point targets in the region with reflectivities at positions for . It follows that
Here, , and and are defined just like and in (7), but evaluated on instead. This expression for is a sum of outer products, each of which corresponds to an individual point target. This outer product representation for indicates that signal subspace methods may be effectively used on these matrices for imaging.
4. Quantitative signal subspace method
To combine the matrices formed using the Prony method, we consider the block diagonal matrix,
Suppose we compute the singular value decomposition for each block:for . According to (8), for point targets, the rank of each block will be . Assuming that the signal-to-noise ratio (SNR) is sufficiently high, the first singular values residing in the diagonal entries of will be significantly larger than the others. Those first singular values correspond to the signal subspace. The remaining singular values correspond to the noise subspace. Since we can separate the first singular values, we are able to compute the pseudo-inverse,
with denoting a user-defined parameter.
For search point in the imaging region, we introduce the illumination block-vector
We also consider another imaging functional. For that imaging functional, we introduce the complimentary illumination block-vector,
We form images through evaluation of and over an imaging region. We show below that the image formed using is useful for determining the location and magnitude of reflectivities for point targets and the image formed using is useful for determining the complex reflectivities for point targets.
5. Resolution analysis
To study the performance of imaging using (12) and (14), we consider one point target located in a planar imaging region at position with complex reflectivity . We use a coordinate system in which the origin lies at the center of the planar imaging region. The flight path of the platform is linear and parallel to the -axis. It is offset from the origin along the -axis by range and along the -axis by height . Thus, spatial positions of the measurements are for with and denoting the aperture. Let denote the distance from the center of the flight path to the origin, and let denote the so-called look angle with and . We assume that is the largest length scale in this problem. A sketch of this linear flight path over a planar imaging region is shown in Fig. 2.
To establish estimates for the resolution of image of one point target produced through evaluation ofover an imaging region, let denote the target location, denote the search location, and denote the system bandwidth centered at frequency with corresponding wavenumber . With these quantities defined, we prove the following theorem.
Theorem 5.1 (Resolution estimates for a linear flight path).
Assuming that the SNR is sufficiently high that we can distinguish the singular values corresponding to the signal subspace from those corresponding to the noise subspace, with given in (12) attains a maximum of on , and in the asymptotic limit, , , , and , this image has a cross-range resolution of and a range resolution of .
For a single point target, we have for with and and given in (7). Consequently, the column space of is and is the projection onto subspace orthogonal to . Using
we find that
where we have introduced the quantity,
with denoting the difference in travel times for the search and target locations.
Evaluating (15) on , we find that , so . Because and with both functions evaluating to only at , this result corresponds the maximum value that attains.
with denoting the fraction of the bandwith about the central frequency. Substituting these frequencies into (16) and computing the sum, we find
from which it follows that
In the expression above, we have resubstituted . Assuming we are in a small neighborhood about the target location, we expand the expression above about and obtain
Additionally, we find that
where and . Using these approximations, we find that
Next, we use
For the cross-range resolution, we evaluate (18) on and find that
Using , we find that
The full-width/half-maximum (FWHM) in cross-range satisfies . Substituting the approximation above into this definition, solving for , and expanding that result about , we find
For the range resolution, we evaluate (18) on and find that
It follows that
The full-width/half-maximum (FWHM) in range satisfies . Substituting the approximation above into this definition, resubstituting , solving for , and expanding that result about , we find
This completes the proof. ∎
Theorem 5.1 states that images of a point target formed through evaluation of with given in (12) will form an image that is peaked at the location of the target with magnitude equal to . Because the user-defined parameter can be made arbitrarily small, this imaging method will yield high-resolution images provided that there is sufficient signal that the non-trivial singular values provide accurate quantitative data.
In general, the reflectivity of a point target is complex. To recover the complex reflectivity, we make use of the following theorem.
Theorem 5.2 (Recovery of the complex reflectivity).
For a point target located at with complex reflectivity , when the SNR is sufficient high that we can distinguish the signal subspace from the noise subspace, with given in (14).
Through direct evaluation of given in (14) on , we find . It follows that . ∎
Although Theorem 5.2 states that evaluating yields the complex reflectivity, it is not generally useful for determining the location of the target because this function does not exhibit localized behavior that indicates the region about the target location. For this reason, we propose the following two-stage imaging method.
6. Travel time uncertainty
We now consider the effect of uncertainty in the travel times on images formed through evaluation of with given in (12). Uncertainty in travel times can arise from sampling clock jitter, deviations from the assumed flight path, and random fluctuations in the propagating medium among other practical issues. It is therefore important to understand to what extent images formed using the method described above are useful under uncertain conditions.
To model travel time uncertainty, we use
with denoting the difference in travel times for a homogeneous medium and the vector, denoting a multivariate distribution with and for . Let . Using this model for travel time uncertainty, we prove the following theorem.
Theorem 6.1 (Travel time uncertainty).
Assuming that the SNR is sufficiently high that we can distinguish the singular values corresponding to the signal subspace from those corresponding to the noise subspace, the image formed through evaluation of in a neighborhood about with given in (12) and using (19) with has an expected value whose leading behavior is the result for the homogeneous medium plus a term that is , and has a variance that is .
Since we consider a neighborhood about , we start with (17) and write
with . Substituting (19) yields
Based on our resolution estimates, we introduce the stretched variables for , and obtain
denote the probability density function for. The expected value of the image is then
Assuming that , we expand about and find
denoting the normalized image formed in the homogeneous medium. Substituting this expansion into the integral above for the expected value of the image and using for , we find that
Next, by using the expansion
we determine that
Theorem 6.1 states that when , the leading behavior of the expectation of the image with random perturbations to the travel time is exactly the same as the image in the homogeneous medium. The recovery of the magnitude of the reflectivity , and the resolution estimates of Theorem 5.1 are different by a term that is . Because the variance of the image is , we determine that this image formed is statistically stable.
An immediate consequence of Theorem 6.1 is given in the following corollary.
Corollary 6.1.1 (Resolution with travel time uncertainty).
When is known or can be reliably estimated, one can set the value of so that and Theorem 6.1 will hold.
Setting in this way connects the resolution of the image with the variance of the random perturbations to the travel time.
7. Numerical results
To validate the theoretical results from above, we use numerical simulations to generate data for various scattering scenes. The following values for the parameters are based on the GOTCHA data set . In particular, we have set and , so that . The synthetic aperture created by the linear flight path is . The central frequency is and the bandwidth is . Using , we find that the central wavelength is . The imaging region is at the ground level . We use frequencies so that , and spatial measurements.
7.1. Single point target
We first consider imaging a single point target located at with complex reflectivity on the planar imaging region. Figure 3 shows the image formed through evaluation of with given in (12) with . Measurement noise was added so that the signal-to-noise ratio (SNR) is . The left plot of Fig. 3 shows the color contour plot of the image in a region about the target location. The center plot of Fig. 3 shows the image on as a function of (cross-range), and the right plot shows the image on as a function of (range). These results shown in Fig. 3 show that the image attains its maximum value of corresponding to at the correct target location. The image attains a high resolution due to choice of . Because and , we expect from the resolution estimates given in Theorem 5.1 that the range resolution should be better than the cross-range resolution. This difference in resolution can be observed by noting the values of in the center plot compared to the values of in the right plot of Fig. 3.
In Fig. 4 we show numerically computed FWHM values of (cross-range resolution) and (range resolution) for a single point target when varying (left plot) and (right plot). The blue “” symbols are the computed values of and the red “” symbols are the computed values of both found by numerically determining the FWHM. The solid blue and red curves are the least-squares linear fit through the and data, respectively. For the results shown in the left plot of Fig. 4, all parameters are set to the same values used for Fig. 3, except that , so there is no noise. For the right plot of Fig. 4, we have varied the value of , but all other parameter values are the same as those used for Fig. 3. In these results, we find that for all values of and which is due to the fact that .
The results for cross-range and range resolutions with respect to given in the left plot of Fig. 4 clearly show an behavior which is plotted as a dashed-black curve in the left plot of Fig. 4. The computed least-squares fits are and which numerically validate this behavior.
The results for cross-range and range resolutions with respect to given in the right plot of Fig. 4 clearly show an behavior which is plotted as a dashed-black curve. The least-squares fits are and which numerically validate the behavior.
The behaviors of computed image resolution with respect to and are shown in Fig. 5. For these results, all parameter values are the same as those used for Fig. 3 except that , so there is no noise and is varied in the left plot and is varied in the right plot. The computed range and cross-range FWHM values, and , respectively, are plotted just as in Fig. 4 including the corresponding least-squares fit to lines.
The results for cross-range resolution with respect to shown in the left plot of Fig. 5 clearly show an behavior, which is plotted as a dashed-black curve. The computed least-squares fit to a line is which numerically validates the behavior. In contrast, the range resolution does not vary significantly with . The computed least-squares fit to a line is which quantifies the weak dependence that range resolution has on aperture.
The results for range resolution with respect to shown in the right plot of Fig. 5 clearly show an behavior, which is plotted as a dashed-black curve. The computed least-squares fit to a line is which numerically validates this behavior. In contrast, the cross-range resolution shows a much weaker dependence on . The computed least-squares fit to a line is which quantifies this weak dependence on .
We now show results from evaluating with given in (14). These results use the same parameter values as those used in Fig. 3. When plotting , there is no local behavior to indicate the location of the target. For this reason these images do not provide useful information about the location of targets. However, when we evaluate in a region near the target location, we are able to recover the complex reflectivity. In Fig. 6 we show the real and imaginary parts of in the left plot and of in the right plot. In both plots the actual value is plotted as a black “” symbol. These plots show that when the location of the point target is known, evaluating at the recovered target location provides a method for recovering the complex reflectivity. At the target location, we find which demonstrates a very high accuracy in recovering the complex reflectivity. Provided that the target location is reasonably accurate, the user-defined parameter can be used to regularize these results to enable stable recovery of the complex reflectivity.
In both Theorems 5.1 and 5.2, it is assumed that the SNR is sufficiently high that one can separate the signal subspace from the noise subspace. To investigate the effect of SNR on the recovery of the complex reflectivity, we evaluate for different SNR values and compute the relative error, . These relative error results are shown with as a solid blue curve, as a dashed red curve, and as a dot-dashed yellow curve in Fig. 7. The results in Fig. 7 show that sufficiently high SNR is needed to achieve a high accuracy. Additionally, we observe that larger values achieve higher accuracy for any fixed SNR. This higher accuracy occurs because regularizes thereby stabilizing the recovery of the complex reflectivity. The role of SNR on the resolution becomes more of an issue when imaging multiple targets which we discuss below.
7.2. Multiple point targets
We now consider multiple point targets in the imaging region. We set the origin of the coordinate system to lie at the center of a