Synthetic Aperture Radar Imaging and Motion Estimation via Robust Principal Component Analysis

Synthetic Aperture Radar Imaging and Motion Estimation
via Robust Principal Component Analysis

Liliana Borcea222Computational and Applied Mathematics, Rice University, MS 134, Houston, TX 77005-1892. (borcea@caam.rice.edu and tscallaghan@rice.edu), Thomas Callaghan222Computational and Applied Mathematics, Rice University, MS 134, Houston, TX 77005-1892. (borcea@caam.rice.edu and tscallaghan@rice.edu), and George Papanicolaou333Department of Mathematics, Stanford University, Stanford, CA 94305. (papanicolalou@stanford.edu)
Abstract

We consider the problem of synthetic aperture radar (SAR) imaging and motion estimation of complex scenes. By complex we mean scenes with multiple targets, stationary and in motion. We use the usual setup with one moving antenna emitting and receiving signals. We address two challenges: (1) the detection of moving targets in the complex scene and (2) the separation of the echoes from the stationary targets and those from the moving targets. Such separation allows high resolution imaging of the stationary scene and motion estimation with the echoes from the moving targets alone. We show that the robust principal component analysis (PCA) method which decomposes a matrix in two parts, one low rank and one sparse, can be used for motion detection and data separation. The matrix that is decomposed is the pulse and range compressed SAR data indexed by two discrete time variables: the slow time, which parametrizes the location of the antenna, and the fast time, which parametrizes the echoes received between successive emissions from the antenna. We present an analysis of the rank of the data matrix to motivate the use of the robust PCA method. We also show with numerical simulations that successful data separation with robust PCA requires proper data windowing. Results of motion estimation and imaging with the separated data are presented, as well.

1 Introduction

In synthetic aperture radar (SAR), the basic problem [6, 16, 7] is to image the reflectivity supported in a set on the ground surface using measurements obtained with an antenna system mounted on a platform flying above it, as illustrated in Figure 1. The antenna emits periodically a probing signal and records the echos , indexed by the slow time of the SAR platform displacement and the fast time . The slow time parametrizes the location of the platform at the instant it emits the signal, and the fast time parametrizes the echoes received between two consecutive illuminations .

Figure 1: Setup for synthetic aperture imaging.

The echoes are approximately, and up to a multiplicative factor, a superposition of the emitted signals time-delayed by the round-trip travel-time between the platform and the locations of scatterers on the ground

(1.1)

Here is the wave speed, assumed constant and equal to the speed of light.

A SAR image is formed by superposing over a platform trajectory of length (aperture) the data convolved with the time reversed emitted signal (matched filtered), and then backpropagated to points in the imaging domain using travel times. With high bandwidth probing signals and large flight apertures, SAR is capable of generating images with roughly ten centimeter resolution at ranges of ten kilometers away from the platform. Such resolution cannot be achieved for complex scenes because moving targets appear blurred and displaced in the images. Image formation should therefore be done in conjunction with target motion estimation. In fact, in applications such as persistent surveillance SAR, tracking and imaging the moving targets is one of the primary objectives.

Because targets may have complicated motion over lengthy data acquisition trajectories, the targets are tracked over successive small sub-apertures. Each sub-aperture corresponds to a short time interval over which the target is in approximate uniform translational motion. The problem is to estimate this motion for each time interval in order to bring the small aperture images of the moving targets into focus. Then, the images are superposed to form high resolution images over larger apertures.

The existing algorithms for motion estimation fall roughly into two categories: The first is for the usual SAR setup with a single moving antenna, and the target motion is estimated from the phase modulations of the return echoes [2, 9, 8, 1, 21, 24, 13, 15, 18, 20, 10]. These algorithms assume that all the targets are in the same motion, and are sensitive to the presence of strong stationary targets. The second class of methods uses more complex antenna systems [14, 23, 22], with multiple receiver and/or transmitter antennas. They form a collection of images with the echoes measured by each receiver-transmitter pair, and then they use the phase variation of the images with respect to the receiver/transmitter offsets, pixel by pixel, to extract the target velocity.

In this paper we consider imaging and motion estimation of complex scenes with the usual SAR setup, using a single antenna. We address two challenges: (1) the detection of moving targets in the scene and (2) the separation of data in subsets of echoes from the stationary scene and echoes from the moving targets. The stationary scene can be imaged by itself after such separation, and the motion estimation can be carried out on the echoes from the moving targets alone. We propose and analyze a detection and data separation approach based on the robust principle component analysis (robust PCA) method [4]. Robust PCA is designed to decompose a matrix into a low rank one plus a sparse one. The main contribution of this paper is to show with analysis and numerical simulations that by appropriately pre-processing and windowing the SAR data we can decompose it into a low rank part, corresponding to the stationary scene, and a sparse part, corresponding to the moving targets. Our theoretical and numerical study describes the rank of the pre-processed SAR data as a function of the velocity, location, and density of the scatterers. It specifies in particular how slowly a target can move and still be distinguishable from the stationary scene. It also addresses the question of proper windowing of the data for the separation with robust PCA to work.

The paper is organized as follows: We begin in section 2 with a brief description of basic SAR data processing and image formation. There are two processing steps that are key to the data decomposition: pulse compression and range compression. We illustrate with numerical simulations in section 3 that robust PCA can be used for motion detection and data separation if it is complemented with proper data windowing. The analysis is in section 4. Additional numerical results on data separation, motion estimation and imaging are in section 5. We end with a summary in section 6.

2 Basic SAR data processing and image formation

The antenna emits signals that consist of a base-band waveform modulated by a carrier frequency ,

(2.2)

Its Fourier transform is

(2.3)

with supported in the interval , where is the bandwidth. The support of is the same interval with center shifted at , and its mirror image in the negative frequencies.

Imaging relies on accurate estimation of travel times. Recall that the echoes are approximately, up to some amplitude factors, superpositions of the emitted signals delayed by the round trip travel times between the antenna and the targets in the scene. Suppose that

with a function of a dimensionless argument and compact support in the unit interval . Then has the form of a pulse, and the travel times can be estimated with precision , the pulse support. But such pulses are almost never used in SAR because of power limitations at the antenna [16]. The SAR echoes should be above the antenna’s noise level, so the emitted signals should carry large power. However, the instantaneous power at the antenna is limited, while large net power can be delivered by longer signals, such as chirps. The problem is that the received scattered energy is spread out over the long time support of such signals, making it impossible to resolve the time of arrival of different echoes. This is overcome by compressing the long echoes, as if they were created by an incident pulse. This pulse compression amounts to convolving the data with the complex conjugate of the time reversed emitted signal

(2.4)

The result is as if the antenna emitted the signal

which turns out to be a pulse of support , as explained in detail in [16].

Figure 2: Left: Pulse compressed data from a single scatterer at , recorded by a SAR antenna moving at speed m/s along a linear aperture of m, at range km from . Note the hyperbolic curve defined by the peak. Right: Pulse and range-compressed data with respect to point . The abscissa is the shifted fast time . Note that most of the slow time dependence is removed in the range-compressed data.

Another common data pre-processing step is range compression. It amounts to evaluating the pulse compressed data at times offset by the round trip travel time to a reference point ,

(2.5)

Here is the fast time shifted by , so that

In the Fourier domain, the range compression amounts to removing the large phase of . This is obviously advantageous from the computational point of view. But range compression plays a bigger role in our context. It is essential for the robust PCA method to work. To illustrate this, we show in Figure 2 the pulse compressed echo from a single scatterer, for a linear aperture. The scatterer is at , with the first two coordinates defining the location in the imaging plane and the third the elevation, which is always zero. The coordinates in the imaging plane are range and cross-range, with origin at the reference point. The range is the coordinate of along the direction pointing from the SAR platform (at the center of the aperture) to . The cross-range is the coordinate of in the direction orthogonal to the range.

We plot in Figure 2 the amplitude of as a function of and , and note that the location of the peak, defined by equation

is a hyperbola for the linear aperture. The peak lies on some other curve in the plane for other apertures. The amplitude of the range compressed data is shown in the right plot of Figure 2. We note that the dependence on the slow time has been approximately removed by the range compression. Explicitly, it lies on the curve in the plane defined by equation

This curve is close to the vertical axis because and are close to each other,

Consequently, the matrix with entries , sampled at discrete and , appears to be of low rank and it can be handled by the robust PCA method.

We work with the pulse and range compressed data from now on, and to simplify notation, we drop the prime from the shifted fast time . We borrow terminology from the geophysics literature and call the pulse and range compressed echoes data traces. The image is formed by superposing the traces over the aperture, and backpropagating them to the imaging points using travel times

(2.6)

Here are the discrete slow time samples in the interval defining the aperture along the flight track. The sampling is uniform, at intervals , and

Assuming a large , that is a small , we approximate the sum in (2.6) by an integral over the aperture.

When the aperture is very large, the data in (2.6) is weighted by a factor that compensates for geometrical spreading effects over the long flight track, and thus improves the focus of the image. Here we work with small apertures where geometrical spreading plays no role, which is why there are no weights in the imaging function (2.6).

3 Robust PCA for motion detection and SAR data separation

We begin in section 3.1 with a brief discussion of the robust PCA method. Then, we give in section 3.2 a heuristic explanation of why it makes sense to use it for motion detection and data separation. We also illustrate in section 3.3 the difficulties arising in the separation, and the improvements achieved by proper data windowing.

3.1 Robust PCA

The robust PCA method, introduced and analyzed in [4], applies to matrices that are sums of a low rank matrix and a sparse matrix . It solves a convex optimization problem called principle component pursuit:

(3.1)
(3.2)

with

(3.3)

Here is the nuclear norm, i.e. the sum of the singular values of , and is the matrix 1-norm of . The optimization can be done for any matrix, but the point is that if , with low rank and sparse and high rank, then the principle component pursuit recovers exactly and . The analysis in [4] gives sufficient conditions under which the decomposition is exact. These conditions are bounds on the rank of and the number of non-zero entries in the high rank matrix . They are not necessary conditions, meaning that the decomposition can be achieved for a much larger class of matrices than those fulfilling the assumptions of the theorems in [4]. This is already pointed out in [4].

3.2 The structure of the matrix of range compressed SAR data

We illustrate in this section the pulse and range compressed echoes (the traces) from a complex scene. The point of the illustration is to show that typically, the sampled traces from the stationary scene form a low rank matrix, whereas those from moving targets give a high rank but sparse matrix. We refer to section 5 for the description of the setup of the numerical simulations used to produce the results presented here and in the next section.

Figure 3: Desired separation of the range and pulse compressed SAR data. Left: for a sample scene with thirty stationary targets. The vertical traces correspond to the stationary targets, while the sloped trace corresponds to the moving target. Middle: The stationary component of . Right: The moving component of . We plot the absolute value of these matrices only to increase contrast for the figures.

Let be the matrix with entries given by the data traces sampled at the discrete slow and fast times

(3.4)

where

(3.5)

and

(3.6)

The sampling is at uniform intervals and , and and are positive integers satisfying

(3.7)

We assume throughout that and are even and large.

We show in the left plot of Figure 3 the matrix for a scene with a single moving target and thirty stationary scatterers. The ideal separation of this matrix would be

where is given by the echoes from the stationary targets alone and is given by the echo from the moving target. We plot and in the middle and right plot of Figure 3. It appears to the eye that is a matrix with almost parallel columns, so we expect it to be low rank. The matrix is given by the sloped curve in the plane. The range compression does not remove the dependence of the echo from the moving target, as is the case for the stationary ones, and this is why we see this sloped curve. Consequently, has higher rank than . But it is sparse, because we have only one moving target, i.e., one sloped curve. Obviously, the matrix will remain sparse for a small number of moving targets, as well.

Thus, the matrix appears to have the appropriate structure for the robust PCA method to be useful. But we show next that robust PCA alone will not do the job adequately. It must be complemented with proper data windowing.

3.3 Data windowing for separation with robust PCA

We illustrate with two numerical simulations that robust PCA cannot be applied as a black box to the SAR data traces and produce the desired separation of the stationary and moving target parts. However, if it is applied to properly calibrated subsets (or windows) of the data, robust PCA can separate the traces for a large class of complex scenes.

The results of the first simulation are in Figure 4. We have a scene with seven stationary targets at coordinates m, m, m and m in the imaging plane, which is at elevation zero. There is also a moving target that is at location , at slow time corresponding to the center of the aperture. The target moves with velocity m/s in the plane. We refer to section 5 for a detailed description of the setup for the simulations. In the top row of plots in Figure 4 we show the results with robust PCA applied to the matrix of data traces shown on the left. The middle and right plots are the resulting low rank and sparse parts given by the robust PCA algorithm. Due to the fact that the entire matrix is sparse to begin with, much of the stationary target data is captured by the sparse component. Thus, we do not have a good separation. However, the result is very good when we apply the robust PCA algorithm on a smaller fast-time window of the data, as shown in the second row of plots. The number of nonzero entries in the windowed is no longer small relative to its size, and robust PCA separates the moving target trace from the rest.

Figure 4: Robust PCA separation of traces from a scene with seven stationary targets and a moving target. From left to right we show the matrix of traces, the low rank , and sparse components from robust PCA. Top row is robust PCA applied to the entire data, bottom row is robust PCA applied to a windowed part of the data. In each plot the matrices are normalized by the largest value of . The sparse components shown in the right column are plotted in (decibel) dB scale to make the contrast more visible.

In Figure 6 we show the results for a scene with thirty stationary targets and a moving one. The data traces from this scene are displayed in Figure 5. The moving target is as in the previous example. The point of the simulation is to show that when we work in a large fast time window the traces from all the stationary targets may no longer form a low rank matrix. The top row of plots in Figure 6 shows that the sparse component of the matrix , as returned by the robust PCA algorithm, has a large residual part from the stationary targets. The much improved results in the bottom are obtained by applying robust PCA on successive small fast time windows of the data, and then reassembling the separated traces. To give an idea of the size of the windows, the matrix has dimensions by . The matrices are windowed in fast-time and are of size .

Figure 5: SAR scene with thirty stationary targets and a moving one.
Figure 6: Robust PCA separation of traces from a scene with thirty stationary targets and a moving one. The top row of plots shows the results with robust PCA applied to the entire data matrix in Figure 5. The low rank and sparse components are shown in the middle and on the right. The bottom row of plots shows the results with robust PCA applied to successive small time windows of the data. The sparse components shown in the right column are plotted in dB scale to make the contrast more visible.

The two examples given above show that data windowing plays an important role in achieving a successful data separation with robust PCA. The last simulation also illustrates the role of robust PCA in the detection of the moving target. While the trace from this target is faint and difficult to distinguish from the others over most of the fast time interval in Figure 5, it is clearly visible in Figure 6 after the processing of the traces with robust PCA.

4 Analysis of rank of the matrix of SAR data traces

We present here an analysis of the rank of the matrix of SAR data traces for simple scenes with one or two targets. The goal is to understand how the rank depends on the position of the targets and their velocity. We limit the analysis to at most two targets to get a simple structure of the matrix , for which we can calculate the rank almost explicitly. The numerical results presented above and in section 5 show that the data separation works for complex scenes, with many stationary targets.

We begin in section 4.1 with the scaling regime used in the analysis. We illustrate it with the GOTCHA Volumetric SAR data set [5] for X-band surveillance SAR. The model of the matrix is given in section 4.2. We analyze its rank in section 4.3 for a single target, and in section 4.4 for two targets. We end with a brief discussion in section 4.5.

4.1 Scaling regime and illustration with GOTCHA volumetric SAR

The important scales in the problem are: The central frequency , the bandwidth , the typical distance (range) from the SAR platform to the targets, the aperture , the magnitude of the speed of the targets, the speed of the SAR platform and , the diameter of the imaging set . We assume that they are ordered as follows

(4.1)

and

(4.2)

and we let

(4.3)

We also suppose that the target speed is smaller than that of the SAR platform

(4.4)

As an illustration, consider the setup in GOTCHA Volumetric SAR. The central frequency of the probing signal is GHz and the bandwidth is MHz, so assumption (4.1) holds. The SAR platform trajectory is circular, at height km, with radius km and speed km/h or m/s. One circular degree of trajectory is m. The pulse repetition rate is per degree, which means that a pulse is sent every m, and s. A typical distance to a target is km and we consider imaging domains of radius of at most m, so assumptions (4.2) hold. The target speed is km/h or m/s, so it satisfies (4.4).

For a stationary target we obtain from basic resolution theory (for single small aperture imaging) that the range can be estimated with precision cm, and the cross range resolution is m, with one degree aperture and central wavelength cm. The image of a moving target is out of focus unless we estimate its velocity and compensate for the motion in the imaging function.

4.2 Data model

Our model of the data assumes that the scatterers lying on the imaging surface behave like point targets. We neglect any interaction between the scatterers, meaning that we make the single scattering, Born approximation. The pulse and range compressed data is approximated by

(4.5)

in the case of targets at locations

with reflectivity . Here we used the so-called stop-start approximation which neglects the displacement of the targets during the round trip travel time. This is justified in radar because the waves travel at the speed of light that is many orders of magnitude larger than the speed of the targets. We refer to [2] for a derivation of the model (4.5) in the scaling regime described above.

Because of our scaling assumptions (4.2), we can approximate the amplitude factors as

and obtain the simpler model

(4.6)

where we let

(4.7)

We assume henceforth, for simplicity, that the targets are identical

and write

(4.8)

with

(4.9)

The matrix of traces analyzed below is given by discrete samples of (4.9),

(4.10)

with slow times and fast times defined in (3.5) and (3.6). We also take for convenience a compressed pulse given by a Gaussian modulated by a cosine at the central frequency,

(4.11)

4.3 Analysis of rank of the data traces for one target

In the case of one target at location , the entries of matrix (4.10) are given by

(4.12)

with defined by (4.7). The target is moving at speed , so we have that

(4.13)

where we let

be the location of the target at the time , corresponding to the center of the aperture.

We study the rank of , which is equivalent to studying the rank of the symmetric, square matrix , with entries given by

(4.14)

in terms of the function

(4.15)

Here we used the assumption that is small enough to approximate the Riemann sum over by the integral over . Because the traces vanish for , we extended the integral to the whole real line.

We obtain after a calculation given in appendix A that if the aperture is small enough, the matrix has a Toeplitz structure.

Proposition 1

Assuming that the Fresnel number is bounded by

(4.16)

the matrix is approximately Toeplitz, with entries given by

(4.17)

The dimensionless parameter depends linearly on the velocity in the range direction and the cross-range offset of the target with respect to the reference point . It is given by

(4.18)

with the unit vector pointing in the range direction from the center of the aperture

(4.19)

and the orthogonal projection

(4.20)

The unit vector is defined by

(4.21)

It is tangential to the flight track, at the center of the aperture.

Figure 7: The matrix for one stationary target at m. The reference point is . The aperture is m, the range km and the central wavelength is cm. The plot shows that the matrix is essentially constant along the diagonals, it is approximately Toeplitz.

As an illustration, we plot in Figure 7 the matrix for an aperture m, in the GOTCHA regime with central wavelength cm and range km. We have a stationary target at m, and is at the origin. The Fresnel number

(4.22)

is only half the ratio , with m, and yet the matrix is essentially Toeplitz.

Remark: Since

the angle between two rows of the matrix of data traces, indexed by and , is given by

(4.23)

Nearby rows with indices satisfying

are nearly parallel, with a small, rapidly fluctuating angle between them. But rows that are far apart, with indices satisfying

are essentially orthogonal, because the right hand side in (4.23) is approximately zero. When the target is stationary, and its cross-range offset is zero, then , and all the entries in are constant and equal to one. Moreover, all the rows of are parallel to each other, and the matrix has rank one. When we view the data traces we see a vertical line at , as in the left plot of Figure 8. Obviously, the larger is, the closer the indices of the rows that are nearly orthogonal. So we expect the rank of to increase with . This is indeed the case, as we show next. We also illustrate this fact in the right plot of Figure 8, where we show the matrix of data traces for a stationary target that is offset from the reference point . Note that by definition (4.18), is large when the cross-range offset of the target and/or the target speed are large.

Figure 8: Matrix of data traces for a stationary target located at m (left) and at m (right).

4.3.1 Asymptotic characterization of the rank

Because matrix is Toeplitz and large, we can use the asymptotic Szegő theory [17, 3] to estimate its rank. To do so, let us define the sequence with entries

(4.24)

A finite set of this sequence, for indices , defines approximately the diagonals of the matrix ,

(4.25)

Since multiplication of a large Toeplitz matrix with a vector is approximately a convolution, and since convolutions are diagonalized by the Fourier transform, it is not surprising that the spectrum of is defined in terms of the symbol , the coefficients of the Fourier series of (4.24),

(4.26)

With this symbol, we can characterize asymptotically in the limit the rank of , using Szegő’s first limit theorem [3], that gives

(4.27)

Here is the indicator function of the interval and is the number of eigenvalues of that lie in this interval.

We show in appendix B that the symbol is given approximately by

(4.28)

where is defined by

(4.29)

We use this result and (4.27) to obtain an asymptotic estimate of the essential rank, defined by

(4.30)

Here is a small threshold parameter, and is of the order of the largest singular value of . It follows from the Szegő theory [17, 3] that this singular value is given by the maximum of the symbol, which is of the order of . We obtain that for ,

(4.31)

where the last equality follows from direct calculation.

(a) Stationary
(b) Moving
Figure 9: Comparison of the computed and estimated rank of for a single target. Left: Stationary target at various cross-range positions. Right: Moving target with various velocities.

As an illustration, we show with green in Figure 9a the computed rank of the matrix for a single stationary target, versus the cross-range position. Positions with larger cross-range result in larger and thus larger rank, as expected. The asymptotic estimate (4.31) of the rank is shown in blue. It exhibits the same linear growth in cross-range and therefore in , but it is lower than the computed rank. This is because in this simulation is not sufficiently large. We get a good approximation when for , so we can approximate the series (4.26) that defines the symbol by the truncated sum for indices . This is not the case in this simulation, so there is a discrepancy in the estimated rank. However, the result improves when we increase , by increasing the aperture. This is illustrated in Figure 10, where we show the rank normalized by the size of the matrix, and note that it has the predicted asymptote as increases.

The plot in Figure 9b shows the computed rank (in green) and the asymptotic estimate (4.31) (in blue) for a moving target located at at . We plot the rank as a function of the range component of the velocity. The rank increases with the velocity, and therefore with , as expected. Moreover, the asymptotic estimate is very close to the computed one, because in this case the entries in the sequence decay faster with . Finally, notice that even for small velocities, the rank is much larger for the moving target than the stationary target.

Figure 10: Convergence of the rank of normalized by the size . The blue line is the computed value and the green line is the asymptotic estimate.

4.4 Analysis of the rank of the data traces for two targets

In the case of two targets at locations for , the entries of the matrix follow from equations (4.9) and (4.10),

(4.32)

We take for simplicity the case of two stationary targets

Extensions to moving targets with speeds and are straightforward. They amount to redefining the parameters and defined below by adding two linear terms in the target velocity, as in equation (4.18).

The expression of matrix is given in the following proposition. It is obtained with a calculation that is similar to that in appendix A, using the same assumption on the Fresnel number as in Proposition 1.

Proposition 2

Assume that the Fresnel number satisfies the bound (4.16), with since the targets are stationary. The matrix has entries defined by the function

(4.33)

sampled at the discrete slow times. Here we let

(4.34)

and

(4.35)

The structure of the matrix is now more complicated. The first term in (4.33) gives a Toeplitz matrix, as before. The other two, which are due to the interaction between the targets, give matrices that are approximately g-Toeplitz or g-Hankel, depending on the sign of the ratio .

Definition 1

An g-Hankel matrix with shift is defined by a sequence as

The matrix is Hankel when . A g-Toeplitz matrix is defined similarly, by replacing with .

To analyze the spectrum of matrix , we choose the reference point in such a way that

(4.36)

Then, is given by

(4.37)

with Toeplitz matrix

(4.38)

and g-Hankel matrix

(4.39)

Here and are sequences with entries

(4.40)

and

(4.41)

An example of matrix and its Toeplitz and g-Hankel parts is shown in Figure 11.

Figure 11: Components of matrix for two stationary targets located at m and m. On the left we show the Toeplitz part . The next two plots show the g-Hankel part and its transpose . The right plot shows the sum, i.e., the matrix .

If the range offsets are large, so that

the g-Hankel matrix has small entries, and is approximately Toeplitz, as in the single target case. The difference of the travel times between the SAR platform and such targets is larger than the compressed pulse width, and their interaction in (4.33) is negligible. If the range offsets are small, the structure of the matrix is as in equation (4.37), and the estimate of its rank follows from the recent results in [19, 11, 12]. They say that the g-Hankel terms have a negligible effect on the rank in the limit . See appendix C for more details. Thus, in either case, the rank estimate of is given by equation (4.27), in terms of the symbol defined by (4.26), using the sequence with entries (4.40).

Explicitly, the symbol is given by

(4.42)

with defined by

(4.43)

for , and the rank is given by

(4.44)
Figure 12: Data traces from two stationary targets. Left: the targets are at positions m and m. Middle: m and m. Right: m and m.

4.4.1 Illustration

We compare in Figure 13 the computed and estimated rank for two stationary targets. The results on the left are for one target fixed at location m. We vary the location of the other target between m and m. The range separation is m, so that and the g-Hankel matrix has negligible entries. The results on the right are for one target at m and the location of the other varying between m to m. The range separation is m, so that and the g-Hankel matrix is no longer negligible. We see that in spite of being neglible or not, the rank of matrix behaves essentially the same, as predicted by the asymptotic theory. This would be difficult to guess by just looking at the data traces displayed in Figure 12.

The computed and the estimated ranks grow at the same rate with the cross range of the second target, i.e., with . The growth is monotone except in the vicinity of the local minimum corresponding to the targets having exactly the same cross-range. In this special configuration , and therefore . The symbol (4.42) simplifies to

and it exceeds the threshold for in a smaller subset of , than in the general case with . The rank is defined by the size of this set, so it should have a minimum as observed in Figure 12.

Similar to the result in Figure 9 for a single stationary target, there is a discrepancy between the computed and estimated rank, due to the aperture not being large enough. This discrepancy diminishes as we increase and therefore the aperture.

Figure 13: Computed and estimated rank of matrix for two stationary targets. One is fixed at location m and the location of the other varies on the line segment between m and m.

4.5 Discussion

The analysis above shows that the matrix of data traces from a moving target has much higher rank than that from a stationary target. The larger the speed, the higher the rank. The rank of the matrix of traces from a stationary target is smallest, equal to one, when the target is at the same cross-range as the reference point . The rank increases at a linear rate with the cross-range offset from .

Comparing the results for one and two stationary targets, we see that the rank increases. The rank depends strongly on the cross-range offset of the targets. There is a small effect due to the separation of the targets in range, but it becomes negligible in the asymptotic limit .

Although we have not presented an analysis for more than two targets, we observe numerically that the rank increases as we add more and more stationary targets. The implication is that the matrix of traces from a stationary scene with many targets is not in general low rank. This is why the data separation with robust PCA should be done in successive small time windows, with each window containing the traces from only a few stationary targets. These traces give a matrix that is low rank, and thus can be separated from the traces due to moving targets. The simulation results shown in Figure 6 illustrate this point.

5 Numerical simulations

We begin with the setup for the numerical simulations. Then we present three sets of results.

5.1 Setup

We use the GOTCHA Volumetric SAR setup described in section (4.1). The data traces are generated with the model (4.5). In all the simulations but the last one, the point targets are assumed identical, with reflectivity . The images are obtained by computing the function (2.6) at points in the square imaging region of area centered at . The motion estimation results are obtained with the phase space algorithm introduced in [2]. This algorithm requires that we know the location of the target at one instant. We choose it at the center of the aperture, which is why there is no error in the target trajectory at .

The principle component pursuit optimization in the robust PCA is solved with an augmented Lagrangian approach. It requires the computation of the top few singular values and corresponding singular vectors of large and sparse matrices, which we do with the software package PROPACK.

5.2 Simulation 1

The first simulation is for a collection of stationary targets placed randomly in the imaging region, and a single moving target located at m at , and moving in the plane with velocity m/s. The data traces are shown in Figure 5.

The estimated moving target trajectory is shown in the left plot in Figure 14. The blue line corresponds to the true trajectory. The red and green lines are the estimated trajectories with the sparse component of the matrix of traces, as returned by robust PCA with and without windowing. These sparse components are shown in the right plots in Figure 6. The separation with robust PCA is better for the windowed traces, and so is the estimate of the target trajectory. This is more clear in the rigÄht plot of Figure 14, where we show the error of the trajectory.

Figure 14: Estimation of the trajectory of a moving target in a complex scene with 30 stationary scatterers, placed at random in a imaging region. We compare the results obtained with the sparse part of the data traces returned by robust PCA with and without windowing. These sparse parts are shown in the top right and bottom right plots in 6.

In the left plot of Figure 15 we show the image obtained with the data traces, and with exact compensation of the motion of the target. The image is focused at the initial location m of the moving target, as expected. However, the stationary targets are out of focus, and the image appears noisy. The right plot of Figure 15 shows the image obtained with the sparse component of the traces, separated successfully by robust PCA with data windowing. The motion compensation is with the estimated velocity. We note that the artifacts due to the stationary targets are now removed, and the image peaks at the expected location m. There are two ghost peaks, due to the error in the estimated target velocity, but they are much smaller than the peak at the correct location.

Figure 15: Images of a scene with 30 randomly placed stationary scatterers and one moving target with velocity 28 m/s. Left: the image given by the original data with exact compensation of the motion of the target. Right: the image given by the sparse component of the data, separated from the other traces by robust PCA with time windowing. The images are normalized by the largest pixel value and plotted in dB.

5.3 Simulation 2

The second simulation is for a scene with 20 stationary targets and two moving targets. The first moving target is as in the first simulation, The second one is located at m at and moves in the plane with velocity m/s.

The data traces are plotted on the left in Figure 16. The separation with robust PCA is shown in the middle and right plots of Figure 16. Each is normalized by the maximum of , and then plotted on the same color scale. Note how the traces from the two moving targets are separated from the other traces. This simplifies the motion estimation.

Figure 16: Data trace separation with robust PCA for a scene with 20 stationary targets and two moving targets. From left to right we show the matrix of data traces, the low rank , and the sparse parts. In each plot we show absolute values normalized by the largest value of . The sparse component is plotted in dB scale to emphasize the contrast.

5.4 Simulation 3

Our last simulation considers again a scene with 30 stationary targets and a single moving target. This target is like that in simulation one, except that its reflectivity is ten times larger than that of the stationary targets.

The data separation results are in Figure 17. We show in Figure 18 the estimated target trajectory with the original data traces (left plot in Figure 17), and the sparse component (right plot in Figure 17). We obtain as before that the estimation is better after the data separation.

Figure 17: Data trace separation with robust PCA for a scene with 30 stationary targets and a moving one with reflectivity that is ten times stronger than the others. From left to right we show the matrix of data traces, the low rank , and the sparse parts. In each plot we show absolute values normalized by the largest value of . The sparse component is plotted in dB scale to emphasize the contrast.

The left plot in Figure 19 shows the image computed with the original SAR data traces. The image is focused at the stationary targets, but there is a strong artifact (streak), due to the moving target. The middle plot in Figure 19 shows the image computed with the low rank component of the data traces, displayed in the middle in Figure 17. The effect of the moving target is now considerably smaller. The right plot in Figure 19 shows the image obtained with the sparse component of the traces, with motion compensation using the estimated target velocity. There is no artifact due to the stationary targets and the image is focused at the location m, as expected.

Figure 18: Estimation of the trajectory of a moving target in a complex scene with 30 stationary scatterers, placed at random in a imaging region. The reflectivity of the moving target is ten times stronger than that of the stationary ones. Left: estimated target trajectory using the original traces (green) and the separated traces (red). The true trajectory is in blue. Right: errors of the estimated trajectories.
Figure 19: Images of a scene with 30 stationary targets and a moving one with reflectivity that is ten times stronger than the others. Left: image obtained with the original data traces. Middle: image obtained with the low rank component of the data traces, returned by robust PCA. Right: image obtained with the sparse component of the data traces, returned by robust PCA. The motion compensation is with the estimated target velocity.

6 Summary

In this paper we consider the problem of synthetic aperture radar (SAR) imaging of complex scenes consisting of a few moving and many stationary targets. The SAR setup is the usual one with a single antenna mounted on a platform flying above the region to be imaged. With large bandwidth probing signals and long data acquisition trajectories, SAR can produce high resolution images of stationary scenes. However, the presence of moving targets may cause serious degradation of the images. When the targets have moderate speed, they appear displaced and blurred in the images. Fast moving targets create significant artifacts such as prominent streaks.

To bring the images of the moving targets in focus we need to estimate their motion. This is necessarily done with small successive sub-apertures, corresponding to short acquisition times over which the target motion can be approximated by uniform translation. Imaging with motion estimation is difficult for at least the following reasons: First, the echoes from the moving targets may be overwhelmed by those from the stationary scenes, so the targets may be difficult to detect. Second, even if we can detect the presence of a moving target in a complex scene, it is difficult to estimate its motion with the existing algorithms, unless we use multiple receiver or transmitter antennas. The algorithms that work with the usual SAR setup assume that all the targets move the same way and are sensitive to the presence of strong stationary targets. Third, even if we detect and estimate well the target motion, when we compensate for it in the image formation process we may bring the stationary targets out of focus, and thus still get images with significant artifacts.

To address these challenges, we propose a pre-processing step of the SAR data designed to separate the stationary target echoes from those due to the moving targets. The main result of the paper is to show that this can be accomplished with the robust principal component analysis (PCA) algorithm complemented with appropriate data windowing. The robust PCA algorithm decomposes a matrix into a low rank part and a sparse part . In our context, the matrix is given by the pulse and range compressed echoes received at the SAR platform. We show with analysis and numerical simulations that the contribution of the stationary targets to is a low rank matrix when we observe it in a small enough time window. Thus, we may think of it as the component of . The contribution of a few moving targets to is a sparse matrix that has higher rank, depending on the target velocity. Therefore, we expect that robust PCA separates it from the low rank part , due to the stationary targets.

We show with numerical simulations that indeed, robust PCA can accomplish such data separation. But the algorithm cannot be applied as a black box. It must be complemented with proper windowing of the pulse and range compressed SAR data in order to achieve a good separation. We present results for various imaging scenes containing multiple stationary targets, and one or two moving targets that may be stronger or weaker than the stationary ones. For weaker targets, we demonstrate that robust PCA can detect their faint echoes and it separates them from those due to the stationary targets. We also show motion estimation and imaging results with and without the data separation step, in order to demonstrate its importance in achieving good results.

Acknowledgement

The work of L. Borcea was partially supported by the AFSOR Grant FA9550-12-1-0117, by Air Force-SBIR FA8650-09-M-1523, the ONR Grant N00014-12-1-0256, and by the NSF Grants DMS-0907746, DMS-0934594. The work of T. Callaghan was partially supported by Air Force-SBIR FA8650-09-M-1523 and the NSF VIGRE grant DMS-0739420. The work of G. Papanicolaou was supported in part by AFOSR grant FA9550-11-1-0266.

Appendices

Appendix A Single target covariance matrix

We approximate here the function

and simplify notation as

We rewrite the integrand using a trigonometric identity, and completing the square