Low rank plus sparse decomposition of synthetic aperture radar data for target imaging and tracking

TL;DR

This paper presents a method using low rank plus sparse decomposition via RPCA to separate stationary and moving targets in SAR data, with a theoretical analysis guiding parameter choice and target velocity detection.

Contribution

The paper provides a theoretical analysis for optimal RPCA parameter selection in SAR data separation and establishes a lower bound for target velocity detection.

Findings

Effective separation of stationary and moving targets demonstrated

Rank of sparse component proportional to square root of target speed

Theoretical guidelines for parameter choice improve separation stability

Abstract

We analyze synthetic aperture radar (SAR) imaging of complex ground scenes that contain both stationary and moving targets. In the usual SAR acquisition scheme, we consider ways to preprocess the data so as to separate the contributions of the moving targets from those due to stationary background reflectors. Both components of the data, that is, reflections from stationary and moving targets, are considered as signal that is needed for target imaging and tracking, respectively. The approach we use is to decompose the data matrix into a low rank and a sparse part. This decomposition enables us to capture the reflections from moving targets into the sparse part and those from stationary targets into the low rank part of the data. The computational tool for this is robust principal component analysis (RPCA) applied to the SAR data matrix. We also introduce a lossless baseband transformation of the data, which simplifies the analysis and improves the performance of the RPCA algorithm. Our main contribution is a theoretical analysis that determines an optimal choice of parameters for the RPCA algorithm so as to have an effective and stable separation of SAR data coming from moving and stationary targets. This analysis gives also a lower bound for detectable target velocities. We show in particular that the rank of the sparse matrix is proportional to the square root of the target's speed in the direction that connects the SAR platform trajectory to the imaging region. The robustness of the approach is illustrated with numerical simulations in the X-band SAR regime.

Figures

Equations

$$\begin{array}{rl}{D}_r(s,t)& =\int \mathrm{d}{t}^{\prime }D\left(s,t-t^{\prime }+\tau (s,{\vec{\mathbf{\rho }}}_o)\right)\overline{f(-t^{\prime })}\\ & =\int \frac{\mathrm{d}\omega }{2\pi }\overline{\hat{f}(\omega )}\hat{D}(s,\omega )e^{-\mathrm{i}\omega [t+\tau (s,{\vec{\mathbf{\rho }}}_o)]},\end{array}$$

$$\begin{array}{rl}\hat{f}(\omega )& =\int dtf(t)e^{i\omega t}.\end{array}$$

$$\tau (s,{\vec{\mathbf{\rho }}}_o)=2\frac{|\vec{\mathbf{r}}(s)-{\vec{\mathbf{\rho }}}_o|}{c},$$

$$I^{\text{S}AR}(\vec{\mathbf{\rho }})=\sum _{j=-n/2}^{n/2}{D}_r(s_j,\tau (s_j,\vec{\mathbf{\rho }})-\tau (s_j,{\vec{\mathbf{\rho }}}_o)).$$

$$f_p(t)=\int d{t}^{\prime }f(t^{\prime })\overline{f(t^{\prime }-t)}.$$

$$D_r(s,t)=\sum _{i=1}^N{\sigma }_i{f}_p(t-\Delta {\tau }_i(s)).$$

$$\Delta {\tau }_i(s)=\frac{2(\Vert \vec{\mathbf{r}}(s)-{\vec{\mathbf{\rho }}}_i(s)\Vert -\Vert \vec{\mathbf{r}}(s)-{\vec{\mathbf{\rho }}}_o\Vert )}{c}.$$

$$f_p(t)=\cos ({\omega }_{o}t)f_B(t).$$

$$\begin{array}{rl}{\hat{f}}_p(\omega )& =\int dt{f}_p(t)e^{i\omega t}\\ & =\frac{1}{2}\left[{\hat{f}}_B(\omega +{\omega }_o)+{\hat{f}}_B(\omega -{\omega }_o)\right].\end{array}$$

$$D_r(s,t)=\sum _{i=1}^N{\sigma }_i\cos ({\omega }_o(t-\Delta {\tau }_i(s)))f_B(t-\Delta {\tau }_i(s)),$$

$${\hat{D}}_r(s,\omega )=\sum _{i=1}^N\frac{e^{-i\omega \Delta {\tau }_i(s)}}{2}\left[{\hat{f}}_B(\omega +{\omega }_o)+{\hat{f}}_B(\omega -{\omega }_o)\right].$$

$$\begin{array}{rl}& {\mathcal{D}}_{il}=D_r\left(s_{i-\frac{n}{2}-1},t_{l-1}\right),\\ & i=1,\dots ,n+1,l=1,\dots ,m+1,\end{array}$$

$$s_j=j\Delta s,j=-n/2,\dots ,n/2,$$

$$t_l=l\Delta t,l=1,\dots ,m.$$

$$\underset{\mathcal{L},\mathcal{S}\in {\mathbb{R}}^{n_1\times n_2}}{\min }||\mathcal{L}|{|}_{\ast }+\eta ||\mathcal{S}|{|}_1$$

$$\text{subject to}\mathcal{L}+\mathcal{S}=\mathcal{D}.$$

$$\eta =\frac{1}{\sqrt{\max \{n_1,n_2\}}}.$$

$${\hat{D}}_r(s,\omega )=\frac{e^{-i\omega \Delta \tau (s)}}{2}\left[{\hat{f}}_B(\omega +{\omega }_o)+{\hat{f}}_B(\omega -{\omega }_o)\right].$$

$${\hat{D}}_o(s,\omega )=\frac{e^{-i(\omega -{\omega }_o)\Delta \tau (s)}}{2}\left[{\hat{f}}_B(\omega )+{\hat{f}}_B(\omega -2{\omega }_o)\right].$$

$${\hat{D}}_{r\text{B}}(s,\omega )={\hat{h}}_{LP}(\omega ){\hat{D}}_o(s,\omega )=e^{-i(\omega -{\omega }_o)\Delta \tau (s)}{\hat{f}}_B(\omega ).$$

$${\hat{h}}_{LP}(\omega )=\{\begin{array}{ll}2& |\omega |\le 2B\\ 0& |\omega |>\alpha B\end{array},$$

$$D_{r\text{B}}(s,t)=e^{i{\omega }_o\Delta \tau (s)}{f}_B(t-\Delta \tau (s)).$$

$$D_r(s,t)=\Re \left\{e^{-i{\omega }_{o}t}{D}_{r\text{B}}(s,t)\right\}.$$

$$\begin{array}{rl}& {\left({\mathcal{D}}_{\text{B}}\right)}_{il}=D_{r\text{B}}(s_{i-\frac{n}{2}-1},t_{l-1}),\\ & i=1,\dots ,n+1,l=1,\dots ,m+1.\end{array}$$

$$\Vert \mathcal{L}{\Vert }_{\ast }+\eta \Vert \mathcal{S}{\Vert }_1+⟨Y,\mathcal{D}-\mathcal{L}-\mathcal{S}⟩+\frac{\mu }{2}\Vert \mathcal{D}-\mathcal{L}-\mathcal{S}{\Vert }_F^2,$$

$${\Theta }_{\eta }(a)=e^{i\arg a}\max (|a|-\eta ,0).$$

$$\begin{array}{r}\Vert A{\Vert }_{\ast }<\eta \Vert A{\Vert }_1,\text{for all low rank matrices }A,\end{array}$$

$$\begin{array}{r}\Vert A{\Vert }_{\ast }>\eta \Vert A{\Vert }_1,\text{for all sparse matrices }A.\end{array}$$

$$\begin{array}{rl}& \eta >\frac{\Vert A{\Vert }_{\ast }}{\Vert A{\Vert }_1},\text{ for all low rank matrices }A,\\ & \eta <\frac{\Vert A{\Vert }_{\ast }}{\Vert A{\Vert }_1},\text{ for all sparse matrices }A.\end{array}$$

$$\begin{array}{rl}& \text{ with }{\eta }_{\min }=\underset{A\text{ low rank}}{\sup }\frac{\Vert A{\Vert }_{\ast }}{\Vert A{\Vert }_1},{\eta }_{\max }=\underset{A\text{ sparse}}{\inf }\frac{\Vert A{\Vert }_{\ast }}{\Vert A{\Vert }_1}.\end{array}$$

$$\begin{array}{rl}& \underset{\eta }{\min }F(\eta )=\frac{{\eta }_{\min }}{\eta }+\frac{\eta }{{\eta }_{\max }}\\ & \Rightarrow {\eta }^{\ast }=\sqrt{{\eta }_{\min }{\eta }_{\max }},F({\eta }^{\ast })=2\sqrt{\frac{{\eta }_{\min }}{{\eta }_{\max }}}.\end{array}$$

Full text

1 Introduction

1.1 Synthetic aperture radar imaging

Synthetic aperture radar can provide high resolution imaging in several applications [6]. The SAR data are usually acquired by a moving platform that probes a remote area by sending broadband pulses, $f(t)$, and recording the echoes. Although there is only one transmit-receive element, high resolution images are obtained by coherently processing the data collected by the moving platform along a large synthetic aperture. Let $\vec{\bf r}(s)$ denote the platform location at time $s$, called the slow time. The data collected form a matrix $D(s,t)$, which records the echoes received at $\vec{\bf r}(s)$ when a pulse is emitted from the same location. We consider $n+1$ discrete slow times with sampling rate $\Delta s$, $s_{j}=j\Delta s$, with $j=-n/2,\cdots,n/2$. The fast time, $t\in(0,\Delta s)$ runs between successive pulse emissions. A schematic of a SAR imaging configuration is shown in Figure 1.

To maximize the power delivered to the region to be imaged, the probing signals used in SAR are long pulses of support $t_{c}\gg 1/B$ where $B$ denotes the bandwidth of the probing pulse. An example of such signals are linear frequency modulated chirps [6]. To concentrate the energy of the reflected echoes to a small time interval of size $1/B$ the received signal is convolved with the complex conjugate of the time-reversed emitted pulse. This step is called pulse compression and is followed by a range compression step which consists of migrating the pulse compressed data to a reference point $\vec{\boldsymbol{\rho}}_{o}$ in the imaging window. The range-compression step removes from the data the large phase $\omega\tau(s,\vec{\boldsymbol{\rho}}_{o})$. These two steps together are referred to as down-ramping,

[TABLE]

with $\widehat{f}(\omega)$ (res. $\widehat{D}(s,\omega)$) denoting the Fourier transform of $f(t)$ (resp. $D(s,t)$), defined as

[TABLE]

Here, $\tau(s,\vec{\boldsymbol{\rho}}_{o})$ is the round trip travel time between the platform location at slow time $s$, $\vec{\bf r}(s)$ and the reference point location $\vec{\boldsymbol{\rho}}_{o}$,

[TABLE]

with $c$ the speed of light.

The SAR image is formed by summing coherently the down-ramped data $D_{r}(s_{j},t)$ back-propagated to the imaging point $\vec{\boldsymbol{\rho}}$ using the travel times differences $\tau(s_{j},\vec{\boldsymbol{\rho}})-\tau(s_{j},\vec{\boldsymbol{\rho}}_{o})$,

[TABLE]

The usual SAR image function (1.4) assumes that only reflections from stationary targets are contained in the down-ramped data $D_{r}(s,t)$. Consequently, if moving targets are present in the region to be imaged then their reflections are not correctly back-propagated and this results in images with streaks or ghosts that depend on the reflectivity and the velocity of the moving targets. For a complex scene, with many stationary and moving targets the overall image is severely distorted and neither the stationary nor the moving targets can be imaged or tracked accurately. To address this problem several motion estimation and separation strategies have been developed. Before considering them we describe first the structure of the SAR data matrix.

1.2 The SAR data matrix

We consider a model for the SAR data matrix for a scene comprised of $N$ small point-like targets. We denote $\sigma_{i}$ the reflectivity of the $i$th target and $\vec{\boldsymbol{\rho}}_{i}(s)$ its location at slow time $s$. The SAR platform emits a pulse $f(t)$ and the down-ramped data, as in (1.1), are obtained by convolving the received echoes with $\overline{f(-t)}$. It is therefore as if the antenna emitted the pulse $f_{p}(t)$ defined as

[TABLE]

The down-ramped SAR data matrix can be modeled as a superposition of pulses $f_{p}$ shifted by $\Delta\tau_{i}(s)$ and multiplied by $\sigma_{i}$, the reflectivity of each target

[TABLE]

The difference travel time $\Delta\tau_{i}(s)$ is the round trip travel time from the antenna location $\vec{\bf r}(s)$ to the target $\vec{\boldsymbol{\rho}}_{i}(s)$ from which the round trip travel time from the antenna location $\vec{\bf r}(s)$ to the reference location $\vec{\boldsymbol{\rho}}_{o}$ is subtracted to account for the range compression step of down-ramping

[TABLE]

We used here the start-stop approximation which neglects the targets’ displacement during the round trip travel time. This is justified in radar because the electromagnetic waves travel at the speed of light that is many orders of magnitude larger than the speed of the targets and the platform.

It is usually assumed that $f_{p}(t)$ consists of a base-band waveform $f_{B}(t)$ modulated by a carrier frequency $\nu_{o}=\omega_{o}/(2\pi)$,

[TABLE]

Its Fourier transform is

[TABLE]

Here $\widehat{f}_{B}(\omega)$ is supported in the interval $[-\pi B,\pi B]$, where $B$ is the bandwidth while $\widehat{f}(\omega)$ is supported in $[-\omega_{o}-\pi B,-\omega_{o}+\pi B]\cup[\omega_{o}-\pi B,\omega_{o}+\pi B]$.

The down-ramped SAR data matrix takes the form

[TABLE]

and its Fourier transform is

[TABLE]

Our model for the SAR data matrix $\mathcal{D}\in{\mathbb{R}}^{(n+1)\times(m+1)}$ is defined by the discrete samples of (1.10),

[TABLE]

with slow times $s_{j}$ defined by

[TABLE]

and fast times $t_{l}$ defined as

[TABLE]

Here we assumed that the pulse repetition rate $\Delta s$ is an integer multiple of $\Delta t$ and set $m=\Delta s/\Delta t$.

1.3 Motion estimation and separation strategies

A simple way to separate moving targets from the stationary background is the change detection (CD) method. In video surveillance, this amounts to subtracting the data corresponding to two consecutive frames. Assuming that multi-frame SAR images, which show the same scene at different times, are available, CD can be used in SAR. We refer to [13] where CD is used for detecting moving targets in multi-frame SAR imagery and to [1] where CD was implemented for human motion tracking with through-the-wall radar. A similar method is the Coherent Change Detection (CCD) which also involves performing repeat pass radar data collection. However, the change detection part is now different. Change is detected in CCD by identifying zeros in the cross-correlation between pairs of images [20, 7, 5].

A specific form of CD that has been used extensively in SAR is the Displaced Phase Center Antenna (DPCA) method, in which the platform uses two antennas, synchronized so that the second antenna is following the same trajectory as the first with a certain time delay. Subtracting the two traces eliminates the stationary background echoes, leaving only the echoes from moving targets [18, 6]

Another method that has been used for moving target detection is autofocus [10, 14, 15, 16]. These algorithms exploit the fact that target motion introduces phase errors, which result in smearing of the target’s image. Autofocus algorithms can be used in the classical single channel SAR setup, and they are composed by two steps: first the conventional smeared images are formed and subsequently a phase error estimation and compensation of the image parts that correspond to moving targets is performed.

A lot of progress has been also achieved in space-time adaptive processing algorithms which aim to better detect the moving targets by suppressing the clutter of stationary targets. Examples are the multi-channel SAR [8], the space-time-frequency SAR [9], velocity SAR (VSAR) [11] and dual-speed SAR [21]. More recently compressed sensing and sparsity driven methods made their way into the moving target SAR imaging problem [19]. We refer the reader to [4] for a review on sparsity driven synthetic aperture imaging including moving target imaging.

1.4 Robust principal component analysis

We follow here the approach proposed in [2], where the robust principal component analysis was used for separating the data matrix $\mathcal{D}$ in (1.12) in two subsets: the echoes due to stationary targets that form the low rank part of the data matrix and those due to moving targets echoes which constitute the sparse part. The RPCA method consists of solving the following convex optimization problem,

[TABLE]

Here $||{\cal L}||_{*}$ denotes the nuclear norm, that is the sum of the singular values of ${\cal L}$, and $||{\cal S}||_{1}$ is the matrix $\ell_{1}$-norm of ${\cal S}$. Assuming the matrix $\mathcal{D}\in\mathbb{R}^{n_{1}\times n_{2}}$ is the sum of a low rank matrix ${\cal L}_{o}$ and a sparse matrix ${\cal S}_{o}$ then the solution of the optimization (1.15)-(1.16) recovers ${\cal L}_{o}$ and ${\cal S}_{o}$ exactly, under sufficient conditions. This optimization problem has been analyzed in [3], where the sufficient conditions for the decomposition to be exact are derived.

The parameter $\eta$ that balances the ratio between the nuclear norm of ${\cal L}$ and the 1-norm of ${\cal S}$ will play an important role in our analysis. The main idea being that we would like to adequately chose the value of $\eta$ given that ${\cal L}$ and ${\cal S}$ represent the echoes received from stationary and moving targets. The recommended value for $\eta$ proposed in [3] and used for SAR data in [2] is

[TABLE]

It was shown in [2] that, with this choice of $\eta$, the RPCA algorithm is sensitive to the window size of the data. Moreover, for complex scenes successful separation between moving and stationary targets’ echoes can be achieved only with adequate windowing of the SAR data.

1.5 Summary of main results

We determine in this paper an optimal range for the parameter $\eta,$ by taking into consideration the specific form that the ${\cal L}_{o}$ and ${\cal S}_{o}$ matrices take in the SAR problem. Specifically, by analytically computing the nuclear and $\ell_{1}$ norms for a single stationary and moving target, we obtain a range of values for $\eta$ for which robust separation results can be obtained with RPCA. Furthermore, our analysis shows that the range of acceptable values for $\eta$ increases with increasing velocity of the moving target. This is intuitively expected since for increasing target velocity the echoes from the moving target become more sparse or equivalently less low rank. Thus, they can be more easily separated from the stationary targets’ echoes. This implies in particular that targets that are moving with a speed above a certain threshold can be easily detected. More precisely, it is the projection of the velocity along the direction that connects the SAR platform to the imaging domain that matters. This is because for targets that are moving parallel to that direction, smaller variation in travel times is expected.

Our analysis suggests that the RPCA algorithm should be applied on the baseband data matrix $\mathcal{D}_{\text{\tiny B}}$ which is obtained from the down-ramped data, $\mathcal{D}$, using a lossless transformation that consists of shifting the signal in the frequency domain so that it is centered at $\omega=0$ instead of $\omega=\omega_{o}$. Going to baseband allows us to evaluate the performance of the RPCA algorithm for a wider range of moving target speeds. This transformation also improves the accuracy of our theoretical analysis and our numerical results indicate that it also improves the robustness of RPCA.

We derive a closed form expression for an optimal value of the RPCA parameter, $\eta^{*}$, for both the baseband and the original SAR data matrix. We show that using $\eta^{*}$ does yield results that are more robust, especially for the baseband case. Robustness is observed with respect to the window size (aperture) of the SAR data, the complexity of the scene and the moving target’s velocity.

2 Robust PCA for motion detection and SAR data separation

We present here the robust PCA algorithm adapted to the SAR problem. We first explain the baseband transformation. The performance of RPCA on the baseband data matrix is theoretically equivalent to the one on the original data matrix, but going to baseband greatly simplifies the analysis. We also see numerically that it yields more robust behavior, in a sense that will be explained in the analysis section 3.

2.1 Baseband data matrix

The proposed preprocessing takes the original data matrix and creates a baseband (BB) version of it by performing band pass filtering at every slow time $s$.

Assuming the model of (1.6) for a single scatterer, the Fourier transform of the down ramped data is

[TABLE]

Thus, the Fourier transform of $D_{o}(s,t)=e^{i\omega_{o}t}D_{r}(s,t)$ is

[TABLE]

We then filter out the signal centered around $2\omega_{o}$ to get a baseband signal centered around $\omega=0$,

[TABLE]

There is no intersection between the supports since we assume $\omega_{o}\gg B$ (narrow band regime). Thus, the filter can be e.g.,

[TABLE]

where $\alpha B<2\omega_{o}-B$ is small enough to eliminate the signal centered around $2\omega_{o}$. The factor of 2 ensures the maximum absolute value remains 1.

In the time domain we will get a filtered signal of the form

[TABLE]

No data is lost in the process, and the original data matrix can be reconstructed by taking the real part of the BB data multiplied by $e^{-i\omega_{o}t}$ ,

[TABLE]

We define $\mathcal{D}_{\text{\tiny B}}$, the baseband data matrix, by applying the Discrete Fourier Transform (DFT) to the data matrix $\mathcal{D}$ with respect to its columns, and applying a similar discrete filter to go to baseband. Assuming a small enough $\Delta t$ the DFT is a good approximation of the Fourier transform and we have

[TABLE]

2.2 Modification of RPCA for complex matrices

Following [17] we use the inexact augmented Lagrangian implementation of RPCA as described in Algorithm 1. This algorithm iteratively approximates the solution by alternatively minimizing the following augmented Lagrangian with respect to ${\cal L}$ and ${\cal S}$

[TABLE]

where $\langle A,B\rangle=\text{tr}(A^{T}B)$ is the regular matrix norm. Each of the minimization steps has an analytical solution given by either spectral- or element-wise thresholding.

The RPCA algorithm can be adapted for use on complex matrices by modifying the thresholding associated with $L^{1}$ minimization to preserve the argument of the number,

[TABLE]

The penalty term coefficient increases with every iteration, usually geometrically $\mu_{k+1}=\rho\mu_{k}$.

Based on [17], initial parameters were set to $\mu_{0}=1.2/\|\mathcal{D}\|_{2}$ and $\rho=1.4$.

2.3 Robustness of RPCA with optimal $\eta$

Before presenting our analysis, we wish to illustrate with a numerical example the robustness of RPCA when using the optimal parameter $\eta$. We consider a scene with five stationary targets and one moving target at a speed $15$ m/s. The moving target’s reflectivity is 5% of the stationary targets’, which makes the detection of the moving target quite challenging. We show the SAR data matrix and its ${\cal L}+{\cal S}$ decomposition in Figure 2. The performance of RPCA with the conventional choice for $\eta$, i.e., $\eta=\frac{1}{\sqrt{\max(n_{1},n_{2})}}$ is not very good when applied to the entire data matrix (see Figure 2(a)) but it is significantly improved when the data are decomposed into $5$ consecutive windows in slow time as seen in Figure 2(b). The plots in Figure 2(c) illustrate that RPCA achieves very good separation without any data windowing when the optimal parameter $\eta^{*}$ is used.

The corresponding Kirchhoff Migration (KM) images are shown in Figure 3. We observe that in the KM image of the original data matrix the moving target is not detectable (see Figure 3(a)). We also see that a lot of clutter is present in the KM image of the sparse matrix generated by the conventional choice of $\eta$ (see Figure 3(c)). This is the result of stationary target’s energy leakage into the sparse part. The KM image of the sparse matrix generated by the optimal choice of $\eta$ is shown in Figure 3(d). We observe a significant improvement in SNR compared to the results obtained with the conventional choice of $\eta$. As expected improvement in the ${\cal L}+{\cal S}$ separation implies increased SNR for both the moving and the stationary targets’ images.

3 Analysis of RPCA for motion estimation

We carry out in this section the analysis of the RPCA algorithm for motion estimation. This will allow us to determine the optimal range of values for the parameter $\eta$.

3.1 RPCA behavior as a function of $\eta$

The conventional value for the parameter $\eta$ is $\frac{1}{\sqrt{\max(n_{1},n_{2})}}$, where $n_{1},n_{2}$ are the dimensions of the data matrix to be analyzed (cf. [3]). This choice is motivated by assuming that the sparse part is random noise, sparsely supported in the data matrix. However, this is not the case for the SAR data.

By taking into account the specific form of ${\cal L}$ and ${\cal S}$ for the SAR problem, we derive an optimal range of values for $\eta$. Our analysis is based on the following observation: The ability of RPCA to separate between stationary and moving targets is due to the difference in their corresponding matrix norms. The two necessary conditions for the algorithm to achieve separation are:

[TABLE]

and

[TABLE]

These conditions suggest the following restrictions on $\eta$

[TABLE]

We can therefore define the following range of values for $\eta,\quad\eta_{\min}\leq\eta\leq\eta_{\max}$,

[TABLE]

Note that (3.24) is a necessary condition but not sufficient. Thus, it provides only an estimate of the range of allowed values for $\eta$. An optimal $\eta$ can be determined by minimizing the objective

[TABLE]

We next analyze each of the terms, i.e., $\eta_{\min}$ and $\eta_{\max}$ to arrive at a closed form expression for $\eta^{*}$. We then show that using this $\eta^{*}$ does yield results that are more robust, meaning that there is no need for fine tuning the windowing of the SAR data and ${\cal L}+{\cal S}$ separation is achieved for a wider range of moving target’s velocities.

We also note that our analysis implies that $\eta$ is independent of any global matrix normalization, i.e. $\frac{\|\lambda A\|_{*}}{\|\lambda A\|_{1}}=\frac{\|A\|_{*}}{\|A\|_{1}}$, and, assuming that there is little variance in the norm of the rows of the SAR data matrix, $\eta$ is independent of any row normalization as well.

3.2 Estimation of matrix norms in the SAR model

We choose as a representative of each class of matrices, i.e. sparse and low rank, a data matrix for a single target, either stationary or moving. We then compute the corresponding matrix norms in both cases and arrive at analytical expressions for $\eta_{\min}$ and $\eta_{\max}$ as in (3.24). We also show with numerical simulations that these matrices which correspond to a single stationary and moving target are indeed good representatives of the classes of low rank and sparse matrices for SAR. Using this norm evaluations we get analytical expressions for $\eta_{\max}$ and $\eta_{\min}$ and hence we can compute the optimal value of $\eta$ for the SAR motion separation problem. In the following sections we carry out evaluations of the matrix $\ell_{1}$ and nuclear norms for stationary and moving targets. The main results are presented here and the reader is referred to Appendix C for the detailed technical calculations.

3.2.1 Estimation of the $\ell_{1}$ norm for the original SAR data matrix

The $\ell_{1}$ norm is given by

[TABLE]

From here on, we assume for convenience a Gaussian baseband pulse

[TABLE]

and given our model form (1.10) we have

[TABLE]

The sum for every time trace can be well approximated as a Riemann integral,

[TABLE]

We notice that if $\Delta\tau(s_{i})\in[t_{0},t_{m}],\forall i,\ \Delta t\ll\{1/B,\omega_{o}\}$, then the absolute sum of elements is approximately the same for all time traces,

[TABLE]

Here $2S(a)=n\Delta s$, is the aperture size (in time units), along the trajectory. In Appendix B we provide an estimate of this integral using the Poisson summation formula. We get

[TABLE]

Since $B\ll\omega_{o}$, the approximation is highly accurate. We also notice that there is no difference in the value of the $\ell_{1}$ norm between stationary and moving targets. Thus, we infer that, in the SAR case, the main difference between stationary and moving targets is due to their corresponding nuclear norms.

3.2.2 Estimation of the $\ell_{1}$ norm for the baseband SAR data matrix

The $\ell_{1}$ norm is given by

[TABLE]

Given our model form (2.18) and assuming a Gaussian pulse as in (3.27) we have

[TABLE]

Similar to (3.28), (3.29), we can approximate this sum using a Riemann integral as follows,

[TABLE]

We notice again that there is no difference in the value of the $\ell_{1}$ norms between stationary and moving targets. Notice also that the $\ell_{1}$ norm is larger by a factor of $\pi/2$ compared to the original SAR data matrix (see (3.30)).

3.3 Spectral properties of the data matrix

In order to evaluate the spectrum of the SAR data matrix we look at the eigenvalues of $\mathcal{D}^{T}\mathcal{D}$, whose $(j,k)-$th element is the inner product of the $j$-th and $k$-th columns of $\mathcal{D}$. By asymptotic analysis we see that for a single target the form of

[TABLE]

is of a block Gaussian matrix. Taking into account the effect of the target’s speed qualitatively, we see that the magnitude of the diagonal elements of $\mathcal{M}$ will be inversely proportional to the slope of the target’s trace, while the total block size of the matrix $\mathcal{M}$ is proportional to the target’s speed. This is true because the faster the target is moving, the larger the column support the time traces will have in the data matrix. Taking $B\Delta t\rightarrow\infty$, we can approximate $\mathcal{M}$ by a diagonal matrix, whose size and elements are related by

[TABLE]

Hence in this case the spectrum of $\mathcal{M}$ is simply $\sqrt{\frac{\alpha}{N}}$, and the nuclear norm of $\mathcal{M}$ is given by

[TABLE]

i.e, the nuclear norm grows with the square root of $N$ and consequently with the square root of the target’s speed. Notice that the $\ell_{2}$ norm of the singular values does not change and is always equal to $\alpha$. This is what is expected, because this is the Frobenius norm of $\mathcal{D}$, which is insensitive to translation in the fast time and doesn’t discriminate between a moving and a stationary target. This preliminary analysis suggests that we should expect a growth rate of the nuclear norm with the target’s speed of the form $cv_{t}^{\beta}$ with $\beta\lesssim 0.5$.

More detailed calculations given in Appendix C show that, the spectrum has a Gaussian decay form as in

[TABLE]

where as mentioned already $N$, the column range of the target, is an increasing function of the velocity.

Because the sum of $\sigma_{k}^{2}$ is constrained by the Frobenius norm to be independent of $N$ and since $\sum\limits_{k=0}^{N}e^{-\frac{\pi^{2}k^{2}}{N^{2}B^{2}\Delta t^{2}}}\propto N$, we deduce that the proportionality factor is $\propto\frac{1}{N}$. Thus, the nuclear norm is proprtional to $\sqrt{N}$, as expected. As we shall see, it is this spectral property that makes RPCA a good tool for separating stationary from moving targets in SAR data. The nuclear norm is a sparsity-promoting norm in the spectrum. Thus, even though for our class of matrices the sum of the squares of singular values is invariant with respect to the target’s velocity, the nuclear norm depends on it.

As pointed out in Appendix C, the form of (3.37) is based on a linear approximation of $\Delta\tau(s)$. For targets that are moving at low velocities, or at directions close to parallel to the direction connecting the imaging region and the platform, the spectrum deviates from this Gaussian form. However, it turns out that the nuclear norm is still well approximated by considering that the spectrum has this Gaussian decay form with respect to $N$. This seems to stem from the fact that the singular values always obey the constraint, $\sum\limits_{k}\sigma_{k}^{2}=C$, almost independently of the target’s velocity or location. Thus, the nuclear norm of two spectra will be similar as long as the two spectra have similar support, which explains the robustness of our approximation for the nuclear norm evaluation.

3.4 Optimal Value for the RPCA parameter, $\eta^{*}$

Based on the evaluations of the matrix norms, we obtain the optimal $\eta$, $\eta^{*}=\sqrt{\eta_{\min}\eta_{\max}}$ with

[TABLE]

[TABLE]

where, to leading order,

[TABLE]

is the minimal column support of the moving target in the data matrix, which is in turn determined by the minimal velocity to be resolved. For the derivation of these expressions see Appendix C.

Note that, as expected, the column support $N(v_{t})$ grows with either a larger aperture or a smaller time resolution, as well as with the velocity component in the direction parallel to the vector connecting the platform to the imaging region. Thus, velocity resolution does not depend only on the speed but also on the direction in which the target is moving. Targets moving in a transverse direction would exhibit very small variation in their time delays independently of their speed.

We also see a growth of $\eta^{*}\propto(N(v_{t})B\Delta t)^{1/4}$. Notice that as $N(v_{t})B\Delta t\nearrow\infty$, we approach the diagonal limit of (3.35), which exhibits square root like growth of the nuclear norm with $N$. On the other hand, as $N(v_{t})B\Delta t\searrow 0$, the dynamic range becomes smaller, and it is harder to separate the targets, as the time traces have wider and wider support. In that case, the variance of the traces with an increasing slope of $\Delta\tau(s)$ would be less significant. We conclude that both the target’s trajectory/speed and the probing signal characteristics affect the performance of RPCA.

3.5 Baseband transformation and RPCA

By qualitatively estimating the nuclear norm of matrices associated with both stationary and moving targets, we have determined how the spectrum and, consequently, the nuclear norm depends on the target’s velocity. Using these evaluations, and the corresponding ones for the $\ell_{1}$ norms of the data matrices, we computed an optimal value for the parameter $\eta$ determined by the minimally resolved velocity.

Our asymptotic analysis results suggest that similar performance is expected for RPCA for both the original and the baseband data matrices. However, the asymptotic analysis, used to derive the optimal parameter, $\eta^{*}$, is highly simplified when using the baseband transformed matrix. Moreover, as shown in Figure 4, the agreement between the asymptotic evaluation of the nuclear norm and the corresponding numerically obtained value is better in the baseband case.

Additionally, our numerical simulations provide evidence that the RPCA algorithm in the baseband case has an improved performance in terms of the dynamic range that $\eta$ can take as illustrated by Figure 5. We observe in general a much smoother behavior in the baseband case. We also see very little variance of $\eta$ with changing the number of stationary targets compared with varying the velocity. This justifies our choice of considering a single stationary target as the representative of the class of stationary data matrices. We also see that the dynamic range of $\eta$ for the baseband matrices is larger, suggesting an overall better performance of the separation when applying RPCA to the baseband SAR data matrices.

4 Conclusions

In this paper, we considered the synthetic aperture radar (SAR) imaging problem and focused our attention in the separation of the echoes coming from moving targets from those due to the stationary objects in the scene. Depending on the complexity of the scene, the velocity and the scattering properties of the moving target, this problem can be quite challenging. SAR imaging is designed to account only for stationary targets and the presence of moving objects introduces artefacts which may significantly affect the the image’s SNR. Moreover, the moving target’s echoes should not be considered as noise but as a useful signal that if processed adequately will reveal the velocity and position of the target. We addressed this problem using the robust principal component analysis that separates a matrix $\mathcal{D}$ into its low rank, ${\cal L}$, and sparse ${\cal S}$ parts, i.e., $\mathcal{D}={\cal L}+{\cal S}$, by minimizing the objective functional $\|{\cal L}\|_{*}+\eta\|{\cal S}\|_{1}$. For the SAR problem, the low rank part consists of the stationary targets’ echoes while the sparse part contains the echoes of the moving targets. Our main contribution is a detailed analysis that allows us to provide a range of optimal values for the single parameter $\eta$, which weights the $\ell_{1}$ norm of ${\cal S}$ in the minimization functional. Our analysis suggests that the key quantity is the nuclear norm of the low rank part, $\|{\cal L}\|_{*}$, which is shown to be proportional to the square root of the velocity of the target. More precisely, it is the projection of the target velocity along the direction that connects the SAR platform to the imaging region that matters. This is intuitively expected since targets that are moving slowly or parallel to the SAR platform will appear as stationary. Our analysis also shows that better separation results are obtained with increasing bandwidth and fast time sampling rate. Another important remark is that it is beneficial to apply the RPCA after transforming the SAR data to baseband. This increases the robustness of the separation with respect to the synthetic aperture and the moving targets speed. Our numerical simulations carried out in the GOTCHA regime for X-band SAR surveillance are in very good agreement with our theoretical results.

Appendix A Parameters used in numerical simulations

Simulations were generated using a platform moving linearly at $V_{p}=300$m/s in the $y$ direction, located at $\vec{\bf r}(0)=[7.1\text{km},0,7.3\text{km}]$ at $s=0$. The emitted signal has a frequency $f_{0}=9.6$GHz, and a bandwith of 622MHz ($B=311$MHz). The total aperture size is 237 pulses, with $\Delta s$=0.015s.

For figures 1, and 2, the imaging area is centered around $\vec{\boldsymbol{\rho}}_{o}=[0,0,0]$. The stationary targets are located at $[4.67,-4.35,0]$, $[2.06,9.61,0]$, $[-3.02,10.64,0]$, $[1.27,-11.1,0]$, $[-4.4,-7.81,0]$. The moving target is located at $[-9.43,-3.07,0]$ at $s=0$. All target locations are given in meters.

For figures 4,6 and 7: The moving target is located at $5.8,7.29,0]$ at $s=0$.

For figure 5, the stationary targets location were drawn uniformly in a $10$m$\times 10$m imaging area. The moving target is moving at the $x$ direction.

Appendix B Estimates of the $L^{1}$ norm of the SAR data matrix

We wish to calculate $\sum\limits_{i,j}|\mathcal{D}_{ij}|$, which is

[TABLE]

Assuming $\tau(s_{i})$ is behaving regularly ($\Delta\tau(s_{i})\in[t_{\min},t_{\max}]\forall i$), the sum over different time traces (rows) is the same and independent of the local translation $\Delta\tau(s_{j})$,

[TABLE]

The integral can be recast as a series, $\sum\limits_{n=-\infty}^{\infty}g(n)$,

[TABLE]

Then, using Poisson’s summation formula, we have

[TABLE]

and

[TABLE]

Thus,

[TABLE]

Since in our setting $\omega_{o}\gg B$ we can approximate the integral well as $2\sqrt{\frac{2}{\pi B^{2}}}$. Notice the $\frac{2}{\pi}$ factor with respect to the regular Gaussian integral $\sqrt{\frac{2\pi}{B^{2}}}$, which is the $\ell_{1}$ norm for the baseband data matrix. Therefore we deduce that the $\ell_{1}$ norm is smaller by a factor of $2/\pi$ compared to the baseband data case.

Appendix C Estimates of the spectrum of the SAR data matrix

C.1 Baseband case

Given the data matrix $D_{r{\text{\tiny B}}}(s,t)$, we want to evaluate its spectrum, based on the number of targets (stationary and moving) and the size of the aperture. Following the model presented in Section  2.1, we evaluate the matrix $D_{r{\text{\tiny B}}}^{H}D_{r{\text{\tiny B}}}$, where the general form of $D_{r{\text{\tiny B}}}$ for a single target is

[TABLE]

with

[TABLE]

When evaluating the matrix $D_{r{\text{\tiny B}}}^{H}D_{r{\text{\tiny B}}}$, the complex phase of (2.18) cancels out. We are left with

[TABLE]

with $\phi(s)$ defined as

[TABLE]

We can aproximate the last sum in (C.39) as an integral

[TABLE]

The values of $\Delta\tau$ range as the same value of $t$ (to be shown), thus $B^{2}\max|\Delta\tau|^{2}\gg 1$, and we can evaluate the integral using Laplace’s method

[TABLE]

The equation for $s^{*}$ is

[TABLE]

So,

[TABLE]

The first option is the only possibility if $\Delta\tau$ is close to linear. It would also dominate, as it achieves by definition the minimum value $g(x)=-(t_{j}+t_{k})x+x^{2}$ and $\phi(s)=g(\Delta\tau(s))$. Assuming $\exists s^{*}\in[s^{-},s^{+}]$ we have $\Delta\tau(s^{*})=\frac{t_{j}+t_{k}}{2}$, and

[TABLE]

Hence,

[TABLE]

Note that for $\Delta\tau(s)$ quasi-linear, $\Delta\tau^{\prime}(s)$ is close to a constant, and the matrix has a Töplitz form. If no such stationary point exists, the element would be strongly suppressed. So we expect, for a single stationary target, $D_{r{\text{\tiny B}}}^{H}D_{r{\text{\tiny B}}}$ to have the form of a matrix, with a localized block around the diagonal. Notice $\Delta\tau(s^{*})=\left(\frac{t_{j}+t_{k}}{2}\right)$, so the matrix $D_{r{\text{\tiny B}}}^{H}D_{r{\text{\tiny B}}}(j,k)$ has the approximate general form

[TABLE]

i.e, a Töplitz times Hankel form on a subset of indices of size $N=\frac{\Delta\tau_{\max}-\Delta\tau_{\min}}{\Delta t}$. $\Delta\tau_{\min},\Delta\tau_{\max}$ are the extreme values of $\Delta\tau(s),s\in[s^{-},s^{+}]$.

We can approximate

[TABLE]

If $\Delta\tau(s)$ is approximately linear, then $\phi^{\prime\prime}(s)$ is independent of $s$, $C_{j+k}\equiv C$, and the matrix has a Töplitz form.

The dominant $s$ dependent term is $\frac{2}{c}\left[-\frac{(\vec{\bf r}(0)-\vec{\boldsymbol{\rho}}_{o})\cdot v_{t}s}{\|\vec{\bf r}(s)-\vec{\boldsymbol{\rho}}_{o}\|}\right]$, thus the analysis holds as long the target’s motion has a significant portion perpendicular to the platforms movement.

C.1.1 Evaluation of the matrix’s spectrum

Using results from [12], we can approximate the spectrum for large $N$ using a circulant matrix. The first row of the circulant matrix would be

[TABLE]

The spectrum of $\mathcal{C}$ is given by the Discrete Fourier transform of $\mathcal{C}_{1,j}$, which we again approximate for large $N$ as an integral

[TABLE]

When evaluating the integral we need to take $\omega_{k}\rightarrow\omega_{k}\mod[-\pi,\pi]$, thus the right end of the spectrum actually corresponds to negative values of $\omega$. This would give us the following spectrum (in decreasing value)

[TABLE]

Using the fact that $\sum\limits_{k}\sigma_{k}^{2}=\|D_{r{\text{\tiny B}}}\|_{F}^{2}$, and we have similar calculation to (3.33) $\|D_{r{\text{\tiny B}}}\|_{F}^{2}=\frac{2S(a)\sqrt{\pi}}{\Delta sB\Delta t}$

[TABLE]

Formally

[TABLE]

where $\vartheta_{3}$ is the Jacobi theta function of the third kind. $\vartheta_{3}(0,e^{-\frac{1}{\sigma^{2}}})$ can be well approximated as

[TABLE]

Thus, for large $N$ and $B\Delta t\ll\pi$, we can approximate the sum using the Jacobi theta function as

[TABLE]

i.e., $\alpha\approx\frac{2S(a)\sqrt{\pi}}{\Delta sB\Delta t\left(\frac{NB\Delta t}{2\sqrt{\pi}}+\frac{1}{2}\right)}$. We can compute $N(v_{t})=\frac{\Delta\tau_{\max}-\Delta\tau_{\min}}{\Delta t}$ using our data model.

We can approximate $N(v_{t})$, for targets that have a non negligible component of their velocity in the direction of $\vec{\bf r}-\vec{\boldsymbol{\rho}}_{o}$ using (C.41) as

[TABLE]

If $(\vec{\bf r}(0)-\vec{\boldsymbol{\rho}}_{o}(0))\cdot v_{t}$ gets too small compared to higher order terms, then the linear behavior is no longer accurate, but using (C.41), we get an accurate quadratic approximation.

Thus, the expression for the spectrum would be

[TABLE]

and therefore,

[TABLE]

Since $N\propto|v_{t}|$ we can expect a $\propto\sqrt{v_{t}}$ growth in the nuclear norm with increasing velocity.

For the approximation of (C.42) to be valid we need

[TABLE]

This means that for smaller $N$ we use explicit summation to normalize the spectrum. Figure 4 shows numerical simulations that suggest that this approximation of the nuclear norm holds well even when the matrix deviates from the Töplitz form. Thus our estimate of the spectrum is robust for a broad range of target directions.

C.2 Spectrum of the original data matrix

We can repeat the analysis of the previous section for the original data matrix. Givent that the data model for a single target is

[TABLE]

Evaluating the matrix, we have the expression

[TABLE]

Notice the first term behaves exactly like the term in the baseband case. We approximate the second term again as an integral and obtain

[TABLE]

Stationary phase analysis gives the dominant stationary point $s^{*}_{\pm}=\frac{t_{j}+t_{k}}{2}\mp i\frac{\omega_{o}}{B^{2}}.$

Plugging this back into the expression we get:

[TABLE]

Since $\omega_{o}\gg B$ this term is heavily attenuated, so that we can neglect it, resulting to

[TABLE]

for $t_{j},t_{k}\in[\Delta\tau_{\min},\Delta\tau_{\max}]$, and 0 otherwise. We next approximate $D^{H}D$ as a circulant matrix, and use the DFT to calculate its spectrum, denoted as $\mu_{k}^{2}$. The effect of the cosine factor would be to create two copies of the spectrum of the baseband matrix, $\sigma_{k}^{2}$, each multiplied by a factor of $\frac{1}{4}$. Thus, the spectrum is

[TABLE]

Note that this is consistent with

[TABLE]

Which is the expected result, from the element-wise Frobenius norm computation. This would also mean

[TABLE]

As shown in Figures 4 and 6, simulations suggest that the actual evaluation of the nuclear norm is more accurate for the baseband case over a wider range of parameters.

Acknowledgments

The work of M. Leibovich and G. Papanicolaou was partially supported by AFOSR FA9550-18-1-0519. The work of C. Tsogka was partially supported by AFOSR FA9550-17-1-0238 and AFOSR FA9550-18-1-0519.