Motion Estimation and Imaging of Complex Sceneswith Synthetic Aperture Radar

# Motion Estimation and Imaging of Complex Scenes with Synthetic Aperture Radar

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. (papanicolaou@stanford.edu)
###### Abstract

We study synthetic aperture radar (SAR) imaging and motion estimation of complex scenes consisting of stationary and moving targets. We use the classic SAR setup with a single antenna emitting signals and receiving the echoes from the scene. The known motion estimation methods for such setups work only in simple cases, with one or a few targets in the same motion. We propose to extend the applicability of these methods to complex scenes, by complementing them with a data pre-processing step intended to separate the echoes from the stationary targets and the moving ones. We present two approaches. The first is an iteration designed to subtract the echoes from the stationary targets one by one. It estimates the location of each stationary target from a preliminary image, and then uses it to define a filter that removes its echo from the data. The second approach is based on the robust principle component analysis (PCA) method. The key observation is that with appropriate pre-processing and windowing, the discrete samples of the stationary target echoes form a low rank matrix, whereas the samples of a few moving target echoes form a high rank sparse matrix. The robust PCA method is designed to separate the low rank from the sparse part, and thus can be used for the SAR data separation. We present a brief analysis of the two methods and explain how they can be combined to improve the data separation for extended and complex imaging scenes. We also assess the performance of the methods with extensive numerical simulations.

## 1 Introduction

Synthetic aperture radar (SAR) is an important technology capable of computing high resolution images of a scene on the ground, using measurements made from an antenna mounted on a platform flying over it [8, 18, 9]. The antenna periodically emits a signal and records the echoes, the data that are processed to form an image. The problem setup is illustrated in Figure 1. The data are indexed by the slow-time and the fast-time . The slow-time parameterizes the antenna position at the moment it emits the signal. The fast-time parameterizes the time between consecutive illuminations at and , so that

The recordings are approximately and up to a multiplicative factor a linear superposition of the probing signals time-delayed by the round-trip travel time between the antenna position and the locations of the scatterers (targets) in . Assuming that the waves propagate at the constant speed of light , the round-trip travel time is given by

 τ(s,→ρ)=2|→r(s)−→ρ|/c. (1.1)

To form an image, the data are first convolved with the time reversed emitted signal , a process known as match-filtering or pulse compression. Then, they are back-propagated to search points using the travel times . The image is obtained by superposing the pulse compressed and back-propagated data over a slow time interval defining the synthetic aperture , the distance traveled by the antenna.

SAR systems that cover long apertures and emit broad-band high frequency signals give very good resolution of images of stationary scenes. The resolution can be of the order of ten centimeters at ranges of ten kilometers away. However, the images are degraded when there is motion in the scene. Depending on how fast the targets move, they may appear blurred and displaced, or they may not be visible at all. The motion needs to be estimated, and then compensated in the image formation, in order to bring these targets in focus. This is a complicated task for complex imaging scenes, consisting of many stationary targets and moving ones.

There are basically two approaches to motion estimation. The first determines the motion from the phase modulations of the return echoes [3, 11, 10, 1, 23, 26, 15, 17, 20, 22, 12], assuming either a single target in the scene, or many targets moving the same way. It does not work well in complex scenes. The estimation is sensitive to the presence of strong stationary targets and to targets that have different motion. The second approach uses multiple receiver and/or transmitter antennas [16, 25, 24]. It forms a collection of images with the echoes measured by each receiver-transmitter pair. Then, it extracts the target velocities from the phase variation of the images with respect to the receiver/transmitter offsets.

We are interested in the first approach that uses the classic SAR setup with a single antenna. To extend its applicability to complex scenes, we complement it with a data pre-processing step intended to separate the echoes from the stationary targets and the moving ones. A successful separation allows motion estimation to be carried with the echoes from the moving targets alone, using the algorithms in [3, 11, 10, 1, 23, 26, 15, 17, 20, 22, 12].

We present two data separation methods. The first extends the ideas in [4, 5], developed in a different context of imaging in scattering layered media. The method seeks to subtract the stationary target echoes one by one, using a combination of travel time transformations and a linear annihilation filter. Each travel time transformation is relative to one stationary target at a time, whose location can be estimated from a preliminary image. After the transformation, the echo from this target becomes essentially independent of the slow time, and so it can be annihilated with a linear filter such as a difference operator in . The travel time transformation is undone after the filtering, and the method moves on to the next stationary target.

The second method is based on the robust principle component analysis (PCA) method [7]. We have shown in [2] with analysis and numerical simulations that after appropriate pre-processing and windowing, the samples of the echoes from stationary targets form a low rank matrix. Moreover, the samples of the echoes from a few moving targets form a high rank but sparse matrix. Thus, the echoes can be separated with the robust PCA method, which is designed to decompose a matrix into its low rank and sparse parts [7]. We review briefly the results in [2]. We also explain that in practice it may be useful to use a combination of the two approaches in order to achieve better data separation for extended and complex imaging scenes.

The paper is organized as follows. We begin in section 2 with a brief description of basic SAR data processing and image formation. The data separation methods are described in section 3. We analyze them in section 4 for simple imaging scenes. Then we show with numerical simulations in section 5 that the methods can handle complex scenes. We end with a summary in section 6.

## 2 Preliminaries

We begin with two data pre-processing steps: pulse compression and range compression. They are commonly used in SAR, but in addition, they are an important first step in our data separation methods. We also describe the generic SAR image formation. The goal of the paper is to explain how to complement the image formation with data separation and filtering in order to estimate motion and focus images of complex scenes.

### 2.1 Pulse and range compression

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

 f(t)=cos(ωot)fB(t). (2.1)

The Fourier transform of the signal is

 ˆf(ω)=∫∞−∞dtf(t)eiωt=12[ˆfB(ω+ωo)+ˆfB(ω−ωo)], (2.2)

where the support of is the interval , and is the bandwidth. Ideally, the waveform should be a pulse of short duration, so that travel times can be easily estimated. It should also carry enough power so that the received echoes are stronger than the ambient noise. However, antennas have instantaneous power limitations and they can transmit large net power only by spreading it over long signals, like chirps. Consequently, the measured echoes are long and it is impossible to read travel times directly from them. This difficulty is overcome by compressing the echoes, as if they were due to short emitted pulses.

Pulse compression is achieved by convolving (match filtering) the data with the complex conjugate of the time reversed emitted signal

 Dp(s,t)=∫dt′D(s,t′)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯f(t′−t). (2.3)

Mathematically, this is equivalent to the antenna emitting the signal

 fp(t)=∫dt′f(t′)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯f(t′−t), (2.4)

and receiving the echoes . The signal is a pulse of duration of the order , as shown for example in [18].

Range compression is another important data pre-processing step. It consists of shifting the fast time in the argument of by the travel time to a reference point ,

 Dr(s,t′)=Dp(s,t′+τ(s,→ρo)). (2.5)

Here is the shifted fast-time, satisfying

Range compression is beneficial for numerical computations and it is essential for our data separation methods because it removes partially the dependence of with respect to the slow time. We illustrate this in Figure 2, where we compare the pulse compressed echo from a stationary target before and after the range compression.

The target is at location m and the reference point is at m. The SAR platform moves at speed m/s along a linear aperture m at range km from . The imaging scene lies on a flat surface at zero altitude. We show in the left plot in Figure 2 the amplitude of the pulse compressed echo as a function of and . Note that its peak lies on the curve defined by

 (ct2)2=|→r(s)−→ρ|2.

This curve is a hyperbola because the aperture is linear and the platform moves at constant speed. The right plot in Figure 2 shows the amplitude of the range and pulse compressed echo. Its peak lies on the curve

 ct′2=|→r(s)−→ρ|−|→r(s)−→ρo|,

which looks almost as a vertical line segment, because the dependence on the slow time has been approximately removed by the range compression.

We work henceforth with the range and pulse compressed data and drop the prime from the shifted fast time to simplify the notation. Borrowing terminology from geophysics, we call the data traces or simply the traces.

### 2.2 Image formation

A SAR image is formed by superposing the traces over the aperture and back propagating them to the imaging points using travel times

 I(→ρI) =n/2∑j=−n/2Dr(sj,τ(sj,→ρI)−τ(sj,→ρo)). (2.6)

The slow time samples discretize the interval which defines the aperture along the trajectory of the SAR platform. The sampling is uniform and equals the time interval between tho consecutive signal emissions from the antenna. To simplify the notation, we assume that the origin of the slow time corresponds to the center of the aperture. Thus, is an even integer and

 2S(a)=nΔs.

Equation (2.6) is a simplification of the imaging function used in practice, which includes a weight factor that accounts for geometrical spreading. The weight makes a difference only when the aperture is large. We neglect it in this paper because we focus attention on motion estimation, which is necessarily done over small successive sub-apertures. Targets may have complicated trajectories over long data acquisition times, but we can approximate their motion by a uniform translation over short durations defining small sub-apertures . The idea is that the speed of translation can be estimated sub-aperture by sub-aperture, and then it can be compensated in the image formation. By compensation we mean that if we seek to image a moving target and we estimate its speed to be , we can form an image as

 I→u(→ρI)=n/2∑j=−n/2Dr(sj,τ(sj,→ρI+sj→u)−τ(sj,→ρo)). (2.7)

We repeat the process for successive sub-apertures, and assemble the final image by superposing the results. This superposition requires proper weighting to compensate geometrical spreading effects over long SAR platform trajectories.

Motion estimation and its compensation in imaging works for simple scenes with one or a few targets that move at approximately the same speed. To our knowledge, there is no approach that estimates motion in complex scenes with a single antenna. Moreover, even if a target speed were estimated, the compensation in (2.7) would bring that target in focus but it would blur the others. To address the difficulty of imaging with motion estimation in complex scenes, we propose a data pre-processing step for separating the traces from the stationary targets and the moving ones. If we could achieve such separation, we could carry the motion estimation on the data traces from the moving targets alone. Moreover, we could image the stationary and moving targets separately, and then combine the results to obtain an image of the complex scene.

## 3 Data separation

We present two complementary approaches for separating traces from stationary and moving targets. They are analyzed in section 4, in simple setups. Here we describe how they work and illustrate the ideas with numerical simulations, for a complex scene named hereafter “Scene 1”. It consists of twenty stationary targets and two moving targets with speeds of m/s and m/s, respectively. We refer to section 5 for details on the setup of the numerical simulations. The traces are displayed in the left plot of Figure 6.

### 3.1 Travel time transformations and data separation

The basic idea of the first data separation approach is that if we have an estimate of the location of a stationary target, we can shift the time by the travel time to remove approximately the dependence on the slow time of the trace from that target. Then, we can subtract the trace from the data by exploiting its independence of . To determine , we need a preliminary image of the scene. The assumption is that the stationary targets dominate the scene, and therefore they are visible in the preliminary image even in the presence of moving targets. See Figure 3 for an illustration.

Let be the map taking the range compressed data to the data in the transformed coordinates

 [T→ρe+Dr](s,t)=Dr(s,t+Δτ(s,→ρe)), (3.1)

where

 Δτ(s,→ρe)=τ(s,→ρe)−τ(s,→ρo). (3.2)

Recall that is the reference point in the range compression. To illustrate the effect of the map (3.1), we show on the left in Figure 4 the trace from a target at location . The plot on the right shows the trace after the travel time transformation (3.1), with the ideal estimate . The transformation makes the trace independent of , and thus it can be annihilated by a difference operator in ,

 Ds[T→ρe+Dr](s,t)=1Δs{[T→ρe+Dr](s+Δs,t)−[T→ρe+Dr](s,t)}≈∂∂s[T→ρe+Dr](s,t). (3.3)

This is an extension of the ideas proposed in [4, 5] for annihilating echoes from a layered scattering medium. Other annihilation operators than (3.3) may be used, as well. In practice, where noise plays a role, they should be complemented with a smoothing process. We refer to [4, 5] for a detailed analysis of annihilator filters. For the simulations in this paper the annihilator works well, and there is no noise added to the data.

After taking the derivative, we map the traces back to the coordinates using the travel time transformation . The data filter is a linear operator denoted by , given by composition of , and ,

 [Q→ρeDr](s,t)=[T→ρe−DsT→ρe+Dr](s,t)=[DsDr(s,t′+Δτ(s,→ρe))]t′=t−Δτ(s,→ρe). (3.4)

We demonstrate in Figure 5 the travel-time transformations and annihilation steps for a scene with one stationary and one moving target. Note that the trace from the stationary target located at is completely removed by the mapping , because .

If there are stationary targets in the scene, we can subtract their traces from the data one by one, by iterating the procedure above. Let be the estimated locations of the targets, with . The data filter is the composition of the linear operators (3.4),

 [QDr](s,t)=[Q→ρeN∘Q→ρeN−1∘…∘Q→ρe1Dr](s,t). (3.5)

The estimates are not exact in practice, but the filter (3.4) is robust with respect to the estimation errors. We show this in section 4, where we present a brief analysis of the filters. In fact, our simulations show that it is typical that one application of removes at once the traces from all the stationary targets located near . Thus, the needed number of iterations in (3.5) may be much smaller than the number of stationary targets in practice.

Because antennas measure discrete samples of the data , we need to use interpolation when making the travel-time transformations. We do so by implementing the travel-time transformations as phase modulations in the frequency domain

 [T→ρe+Dr](s,t)=12π∫dω ˆDr(s,ω)e−iω(t+Δτ(s,→ρ)), (3.6)

where

 ˆDr(s,ω)=∫dt eiωtDr(s,t).

We implement (3.6) in the discrete setting using the fast Fourier transform.

### 3.2 Data separation with robust PCA

We introduced recently in [2] a method for separating the traces from stationary and moving targets using robust principal component analysis (PCA). We review here the data separation method, and give some details of its analysis in section 4.

The robust PCA method was introduced and analyzed in [7]. It decomposes a matrix into a low rank matrix and a sparse matrix . In our context, the entries in are the samples of the traces

 Mjℓ=Dr(sj,tℓ), (3.7)

where

 sj=jΔs,j=−n/2,…,n/2,2S(a)=nΔs, (3.8)

and

 tℓ=ℓΔs,ℓ=−m/2,…,m/2,Δs=mΔt. (3.9)

The sampling is at uniform intervals and , determined in practice by the SAR antenna system. To simplify the notation, we take the origin of at the center of the aperture, and the origin of in the middle of the interval between two consecutive signal emissions.

To explain why we can use robust PCA for the data separation, consider the matrix of sampled traces for Scene 1. We display it in the left plot of Figure 6. The traces that appear almost vertical correspond to the stationary targets. This is because these targets are not too far from the reference point , and the range compression removes most of their dependence. We expect that these traces form a low rank matrix . We display them in the middle plot of Figure 6. The range compression does not remove the slow time dependence of the traces from the moving targets, so they appear sloped in Figure 6. They form a high rank but sparse matrix , shown in the right plot of Figure 6.

The robust PCA method is designed to separate from in It does so by solving an optimization problem called principle component pursuit [7]

 minL,S∈Rn+1×m+1∥L∥∗+η∥S∥1,subject toL+S=M. (3.10)

Here is the nuclear norm of , the sum of the singular values of , is the matrix 1-norm of and

 η=1√max{n+1,m+1}. (3.11)

The optimization (3.10) can be applied to any matrix . But it is shown in [7] that if is a matrix with a low rank part and a high rank sparse part , then the algorithm recovers exactly and . We refer to [7] for sufficient (not necessary) conditions satisfied by , so that it can be succesfully decomposed by principle component pursuit.

Robust PCA is well suited for our purpose, but it cannot be applied as a black box to the matrix of traces. We need to complement it with properly calibrated data windowing. To illustrate why, we present two examples.

Example 1: The first example considers a scene with one moving target and seven stationary ones. The moving target is at location m at time , corresponding to the center of the aperture, and velocity m/s. The stationary scatterers are located at m, m, m, and m. in the imaging region. The reference point is m.

When we apply the robust PCA method to the matrix of traces shown in the top left plot in Figure 7, we obtain the “low rank”and “sparse” parts shown in the top middle and top right plots. We note from the right plot that the “sparse” component contains the trace from the moving target but also remnants of the traces from the stationary targets. The data separation is not successful because the matrix is sparse to begin with.

However, if we apply robust PCA to the matrix of traces in a smaller fast-time window, we get a successful separation. This is shown in the bottom row of Figure 7. The windowed has a much larger number of non-zero entries relative to its size, and robust PCA separates the traces from the moving target and the stationary ones.

Example 2: The second example considers a scene with thirty stationary scatterers and the same moving target as in example 1. The matrix of traces is displayed in Figure 8. The results with robust PCA applied to the entire are in the top row of Figure 9. Again, we see that the “sparse” part contains remnants of the traces from the stationary targets. This is not because is sparse, as was the case in example 1, but because the traces from the stationary targets do not form a matrix of sufficiently low rank. When we window , we work with a few nearby traces at a time. These traces are similar to each other so they form a low rank matrix that can be separated with robust PCA. The final result is given by the concatenation of the matrices in each window. It is shown in the bottom row of Figure 9.

These two examples show that windowing is a key step to successful data separation with robust PCA. To choose appropriate window sizes, it is necessary to understand how the rank and sparsity of the data matrix depend on the location and density of the stationary targets, and on the velocity of moving targets. This is addressed in the analysis of [2], which we review in section 4.4.

### 3.3 Discussion

We use numerical simulations for Scene 1 to compare the two methods described in sections 3.1 and 3.2. The results are in Figures 10 and 11.

The first method uses the filter defined in (3.5) to annihilate the traces from the stationary targets. We see from Figure 10 that this is accomplished, but the filter removes some parts of the traces from the moving targets as well. All our numerical simulations show that separates nicely the traces from the moving targets if the scene does not have too many stationary targets at almost the same range. Otherwise, it also removes parts of the traces from the moving targets.

The data separation with the robust PCA approach is better, as seen in Figure 11. This method is more sophisticated than the other, because it involves careful windowing of the traces, but it works better for complex scenes that are not too large. Robust PCA relies heavily on range compression removing most of the slow-time dependence of the traces from the stationary targets. But for extended scenes there is no single reference point that accomplishes this task, and data separation with robust PCA deteriorates. The filter (3.5) still works for extended scenes, because the travel time transformations in (3.4) are relative to the location of the targets that we wish to annihilate, and not some arbitrary reference point.

In principle, we can use a combination of the two approaches to improve the data separation with robust PCA for extended scenes. The idea is to begin with a preliminary image and identify roughly the stationary targets to be separated from the scene. Then we can window the image around groups of such targets and work with one subscene at a time. It is important to observe here that if the aperture is not too large, as is the case in motion estimation, the image is approximately the Fourier transform of the data with a phase correction. See for example [3] for a detailed explanation of this fact. Thus, we can obtain approximately the traces from the smaller subscene by Fourier transforming the windowed image. To separate these traces, we can apply the travel time transformation (3.1) with a reference point in the subscene, followed by the robust PCA approach as described above.

### 3.4 Velocity estimation and separation of traces from multiple moving targets

Let us assume in this section that the traces from the stationary targets have been removed, and denote by the same the traces from the moving targets. Since there may be multiple moving targets, we would like to separate these traces further, based on differences between the target velocities. We can do so by extending the annihilation filters described in section 3.1, as we explain next.

We need a key observation, shown with analysis in section 4. The slope of the trace from a moving target is determined approximately by its speed along the unit vector

 →mo=→r(0)−→ρo|→r(0)−→ρo|, (3.12)

pointing from the SAR platform location at the middle of the aperture to the reference point in the imaging scene. Moreover, the curvature of the trace is determined by the speed projected on the plane orthogonal to .

Because we assume a flat imaging surface, the target speed is a vector in the horizontal plane

 →u=(u,0),u∈R2. (3.13)

We parametrize it by the range velocity and the cross-range velocity . The range velocity is

 u=→u⋅→mo=u⋅mo, (3.14)

where is the projection of on the imaging plane. The cross-range velocity is defined by

 u⊥=→t⋅Po→u=u⋅t−u→mo⋅→t(0), (3.15)

where is the orthogonal projection

 Po=I−→mo→mTo, (3.16)

and is the unit vector tangent to the SAR platform trajectory, at the center of the aperture. Note that in many setups and are almost orthogonal, so is approximately the target speed along the vector , the projection of on the imaging plane.

Next, we define the transformation mapping , the generalization of (3.1),

 [T→ue,→ρe+Dr](s,t)=Dr(s,t+Δτ(s,→ρe+s→ue)), (3.17)

with given in (3.2). Here is the estimated location of the moving target at the origin of the slow time, and is the estimated velocity. If the line segment

 →ρe+s→ue,s∈[−S(a),S(a)],

is close to the target trajectory, then the transformation (3.17) removes approximately the dependence of its trace on the slow time. Thus, we can subtract it from the other traces using the same difference operator . The annihilation filter is the generalization of (3.4),

 [Q→ue,→ρeDr](s,t) = [T→ue,→ρe−DsT→ue,→ρe+Dr](s,t) (3.18) = [DsDr(s,t′+Δτ(s,→ρe+s→ue))]t′=t−Δτ(s,→ρe+s→ue),

where the mapping undoes the travel time transformation. Additionally we can apply robust PCA to . This should result in the low rank matrix containing the trace corresponding to the moving target with velocity and the sparse matrix containing the rest of the moving target traces.

The implementation of (3.17) and (3.18) requires the estimates and . There are two scenarios. The first assumes that we know from tracking the target in previous sub-apertures, and seeks only . The second seeks both and . In either case, the estimation begins with that of the range velocity, which controls the slope of the trace. Then, we can estimate if we need to, and the cross-range speed .

We estimate as the maximizer of the objective function

 (3.19)

The idea is that when is close to the range speed of the target, the transformation makes the trace an approximate line segment, and when we sum it over , we get a peak around some fast time , for . If we do not know , we can set it equal to in (3.19). Then, we can form an image using (2.7) with , and estimate from there. The motion compensation in (2.7) brings the moving target approximately in focus, and this is why we expect to be able to estimate its location from the image.

Once we have estimated and , we seek the cross-range speed as the minimizer of the objective function

 g⊥(u⊥;ue,→ρe)=m/2∑ℓ=−m/2n/2∑j=−n/2∣∣D2s[T→u,→ρe+Dr](sj,tℓ)∣∣, (3.20)

which quantifies the curvature of the trace. Here approximates the second derivative in , and is defined by

 u⋅mo=ue,u⋅t(0)=u⊥+ue→mo⋅→t(0). (3.21)

Because has a weaker effect on the traces than , it is more difficult to estimate it. Even if there is no noise, the estimation is sensitive to remnants of the traces from the stationary targets that have not been perfectly annihilated, the traces from the other moving targets, and the discrete sampling of the data. If we took the maximum over of the sum

 F(u⊥,ℓ)=n/2∑j=−n/2∣∣[D2sT→u,→ρe+Dr](sj,tℓ)∣∣

as we did in the objective function (3.19), we would seek a saddle point of to estimate . That is to say, we would minimize with respect to and maximize with respect to . Our numerical simulations show that such an estimation is not robust. To stabilize it, we sum over in (3.20) and reduce the problem to that of minimizing in .

As an illustration, we present the results for Scene 1. Figure 12 displays using the true locations and velocities for the two moving targets. Note how after the transformation the traces become vertical line segments. If we were to search for the range speed without annihilating the traces of the stationary targets, we would get the objective function plotted on the left in Figure 13. The prominent peak of is at speed , because the stationary targets dominate the scene. The smaller peaks are due to the two moving targets. When we compute the objective function with the traces from the moving targets alone, we get two large peaks, at the range speeds of the targets. The speed in (3.17) is , with estimated from the first and second peak of , respectively.

Recall that the cross-range speed has a weaker effect on the traces than the range speed, it determines only their curvature. We illustrate this in the top right plot of Figure 14, where we zoom around the trace of one target, rotated with transformation (3.17) at the estimated range speed . The bottom row in Figure 14 shows the traces transformed with (3.17), after the cross-range speed correction. The trace is now straightened, as seen in the bottom right plot. The objective function (3.20) is shown in the middle and right plots of Figure 15. We plot it for the two range speeds, the peaks of shown on the left. The function has minima at the cross-range speeds of the two targets.

## 4 Analysis

In this section we give a brief analysis of our data separation approaches. The analysis is for simple scenes, but the methods work in complex scenes, as shown already in the previous section. More numerical results are in section 5. We begin with the data model and the scaling regime. Then, we analyze the annihilation filters defined in section 3.1. The analysis extends to the filters described in section 3.4, once we explain how the target velocity affects the traces. We end the section with a brief review of the analysis in [2] of the data separation with robust PCA.

### 4.1 Scaling Regime

We illustrate our scaling regime using the GOTCHA Volumetric SAR data set. The relevant scales are the central frequency , the bandwidth , the typical range from the SAR platform to the targets, the aperture , the magnitude of the velocity of the targets, the speed of the SAR platform and , the diameter of the imaging set . We let , and assume that

 B≪νo,a≪L,RI≪L. (4.1)

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

 |u|

The GOTCHA parameters satisfy these assumptions:

• The central frequency is GHz and the bandwidth is MHz.

• The SAR platform trajectory is circular, at height km, with radius km and thus km. We consider sub-apertures of one circular degree, which corresponds to m, and image in domains of radius m.

• The platform speed is m/s and the targets move with velocities m/s.

• The pulse repetition rate is per degree, which means that a pulse is sent every m, and s.

For a stationary target, we obtain from basic resolution theory [3, 18] that the range resolution is cm, and the cross range resolution is m, with one degree aperture and central wavelength cm.

### 4.2 Data model

We use a simple data model that approximates the scatterers by point targets in the horizontal plane and neglects multiple scattering between them. We also suppose that the target trajectories can be approximated locally, for the slow time interval defining a single sub-aperture, by the straight line segment

 →ρ(s)=→ρ(0)+s→u,

with constant speed , and . The traces are modeled by

 Dr(s,t)≈N∑q=1σq(ωo)(4π|→r(s)−→ρq(s)|)2fp(t−(τ(s,→ρq(s))−τ(s,→ρo))), (4.3)

where is the number of targets. They are located at and with reflectivity , for . The model neglects the displacement of the targets during the round trip travel time. This is justified in radar applications because the waves travel at the speed of light which is several orders of magnitude larger than the speed of the targets. We refer to [3] for a derivation of model (4.3) in the scaling regime described above. We simplify it further by approximating the amplitude factors

 14π|→r(s)−→ρq(s)|≈14πL,

and assuming that the targets have the same reflectivity

 σq(ωo)=σ(ωo),q=1,…,N.

The model of the traces becomes

 Dr(s,t)≈σ(ωo)(4πL)2M(s,t), (4.4)

with

 M(s,t)=N∑q=1fp(t−Δτ(s,→ρq(s))), (4.5)

and given in (3.2). Thus, the traces are approximately, up to a multiplicative constant, equal to the matrix , for and . We neglect the multiplicative constant hereafter, and call the matrix of traces.

We could take any pulse in (4.5), but to simplify the calculations in the remainder of the section, we assume it Gaussian, modulated by a cosine at the central frequency,

 fp(t)=cos(ωot)e−B2t2/2. (4.6)

### 4.3 The annihilation filter

Recall that the filter defined by (3.4) is a linear operator. Because we model the traces (4.5) as a superposition of pulses delayed by the travel times to the targets, we can study the effect of on one trace at a time.

Consider the trace from a stationary or moving target at location , for some integer satisfying , and denote it by

 Mj(s,t)=fp(t−Δτ(s,→ρj(s))). (4.7)

The transformation defined in (3.1) maps the trace to

 [T→ρe+Mj](s,t)=fp(t+Δτ(s,→ρe)−Δτ(s,→ρjs))), (4.8)

and then we filter it with the difference operator ,

 Ds[T→ρe+Mj](s,t)≈A(s,→ρj(s),→ρe)f′p(t+Δτ(s,→ρe)−Δτ(s,→ρj(s))). (4.9)

Here we approximate the finite difference by the derivative in , denote by the derivative of the pulse, and let be the “annihilation factor”

 A(s,→ρj(s),→ρe)=dds[τ(s,→ρj(s))−τ(s,→ρe)]. (4.10)

Using the scaling assumptions (4.1-4.2), we obtain after a straightforward calculation that

 dds[τ(s,→ρj(s))−τ(s,→ρe)]=2c⎡⎣−→u⋅→me+(V→t−→u)⋅PeΔ→ρjL−s(2V→t−→u)⋅Pe→uL⎤⎦+ O(au⊥cL)+O(aVRIcL2), (4.11)

where

 Δ→ρj=→ρj(0)−→ρe,→me=→r(0)−→ρe|→r(0)−→ρe|,Pe=I−→me→mTe.

The annihilation factor of a stationary target is given by

 A(s,→ρj,→ρe)=2Vc→t⋅PeΔ→ρjL+O(aVRIcL2). (4.12)

It is approximately linear in the cross-range component of the error of the estimated location of the target. For a moving target the factor is

 A(s,→ρj(s),→ρe)≈2Vc⎡⎣→t⋅PeΔ→ρjL−uV−2VsLu⊥V⎤⎦+O(au⊥cL)+O(aVRIcL2). (4.13)

It is much larger than (4.12) for stationary targets that are in the vicinity of , and satisfy the condition

 (4.14)

Consequently, the filter defined by (3.4) emphasizes the traces from the moving targets and suppresses the traces from the stationary targets near , as stated in section 3.1.

###### Remark 1

In the GOTCHA regime, and for one degree apertures as described in section 4.1, the cross-range resolution of a preliminary image is approximately m. Thus, we can assume that for the target that we wish to annihilate m, and condition (4.14) is satisfied for targets moving at range speed that exceeds m/s.

Note also from (4.11), with replaced by , that the slope of the trace is approximately given by the range speed , as assumed in section 3.4, as long as

 u≫V→t⋅PoΔ→ρjL∼VRIL.

If the target is a slower mover, than its trace looks similar to that of a stationary target and it cannot be separated. The curvature of the trace is determined by the cross-range velocity of the target, as seen from the third term in the right hand side of equation (4.11).

### 4.4 Analysis of data separation with robust PCA

In this section we review briefly from [2] the analysis of the rank of the matrix of traces for simple scenes. The goal is to show that indeed, traces from stationary targets should be expected to give a low rank contribution whereas those from the moving targets should give a high rank contribution. This is the key assumption in the data separation approach described in section 3.2.

#### 4.4.1 One target

We begin with the case of one target at location . It moves at constant speed over the slow time interval , so we can write that

 →ρ(s)=→ρ+s→u,|s|≤S(a), (4.15)

where

 →ρ=→ρ(0).

The matrix denoted by is defined by the discrete samples of

 M(s,t)=cos[ωo(t−Δτ(s,→ρ(s)))]exp[−B22(t−Δτ(s,→ρ(s)))2], (4.16)

where is defined by (3.2), and we used the pulse (4.6). The rank of is the same as that of the symmetric, square matrix , given by

 Y=MMT. (4.17)

Its entries are the discrete samples of the function

 Y(s,s′)=m/2∑q=−m/2M(s,tq)M(s′,tq)≈1Δt∫∞−∞dt M(s,t)M(s′,t). (4.18)

We assume here that is small enough to approximate the Riemann sum over by the integral over . Since the traces vanish for , we can extend the integral to the whole real line.

We showed in [2] that under our scaling assumptions, the function (4.18) takes the form

 Y(s,s′)≈√π2BΔtcos[ωoα(s−s′)]exp[−(Bα)2(s−s′)24], (4.19)

with dimensionless parameter

 α=2→u⋅→moc−2V→t⋅Po(→ρ−→ρo)cL+2→u⋅Po(→ρ−→ρo)cL. (4.20)

That is to say, the matrix (4.17) is approximately Toeplitz, with diagonals given by the entries in the sequence ,

 Yjℓ≈yj−ℓ,

where

 yj=√π2BΔte−(BΔs|α|j)24cos(ωoΔs|α|j). (4.21)

There are many slow time samples in an aperture , so the matrix is also large, and we can approximate its rank using the asymptotic Szegö theory [19, 6].

Note that multiplication of a large Toeplitz matrix with a vector can be interpreted approximately as a convolution, which is why the spectrum can be estimated from the symbol , the Fourier transform of (4.21),

 ˆy(θ)=∞∑j=−∞yjeijθ,θ∈(−π,π). (4.22)

Szegő’s first limit theorem [6] gives that

 limn→∞N(n;β1,β2)n+1=12π∫π−π1[β1,β2](ˆy(θ))dθ, (4.23)

where is the indicator function of the interval and is the number of eigenvalues of that lie in this interval. Using this result, we estimated the rank in [2] as

 rank[Y]:=N(n;ϵ∥ˆy∥∞,∞),∥ˆy∥∞=supθ∈(−π,π)|ˆy(θ)|. (4.24)

Here is a small threshold parameter, and is of the order of the largest singular value of . The Szegő theory [19, 6] gives that this singular value is approximated by the maximum of the symbol. The right hand side in (4.23) is calculated explicitly in [2], and the rank estimate is

 rank[Y]n+1 ≈ 12π∫π−π1[ϵ∥ˆy∥∞,∞)(ˆy(θ))dθ (4.25) = min(2|α|BΔs√log1/ϵπ,1).

It grows linearly with , meaning that the larger the offset in the cross-range direction is, and the faster the target moves, the higher the rank.

We show in Figure 16 an illustration of this result. The left plot is for a stationary target and the right plot is for a moving target. The green line shows the numerically computed rank of the matrix , while the blue line shows its asymptotic estimate (4.25). Note that both curves show the same growth. In the case of the stationary target there is a mismatch of the magnitude of the computed and estimated rank. This is because is not large enough in the simulation. A good approximation is obtained when for , because then we can approximate the series (4.22) 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, if we increase the aperture, thus increasing , the result improves. Figure 17 illustrates this fact by showing how the rank normalized by the size of the matrix has the predicted asymptote as increases.

In the right plot of Figure 16 we compare the computed rank (in green) and the asymptotic estimate (4.25) (in blue) for a moving target located at at . The rank is plotted as a function of the range component of the velocity. Note that, as expected, the rank increases with the velocity, and therefore with . For the moving target, the asymptotic estimate is in closer agreement to the computed one. This is because in this case the entries in the sequence decay faster with . Finally, we can observe that the rank is much larger for the moving target than the stationary target, even for small velocities.

#### 4.4.2 Two targets

When there are two targets in the scene, we obtain from model (4.5) and the choice (4.6) of the pulse that

 M(s,t)=2∑j=1cos[ωo(t−Δτ(s,→ρj))]exp[−B22(t−Δτ(s,→ρj))2]. (4.26)

We restrict the discussion to the case of two stationary targets. The extension to moving targets is straightforward. It amounts to adjusting the parameters and defined below in (4.28), by adding two linear terms in the target velocity, as in equation (4.20).

We obtain after a calculation given in detail in [2] that the entries of the matrix (4.17) are the discrete slow time samples of the function

 Y(s,s′) ≈√π2BΔt{2∑j=1cos[ωoαj(s−s′)]exp[−(Bαj)2(s−s′)24] +cos[ωo(α1s−α2s′+β)]exp[−B2(α1s−α2s′+β)24] (4.27)