1 Introduction

Robust Principle Component Analysis (RPCA) has been one of the most popular ways in dealing with huge amount of data in statistics. It also has important applications in many other areas, such as video surveillance, face recognition, latent semantic indexing, ranking and collaborative filtering. We refer to [2, 13, 28] for the review of RPCA. Given matrix D∈ℜm×n, RPCA is to find a low rank matrix L 0 and a sparse matrix S 0, such that

$$ D = L_0 + S_0. $$
(1)

In [2, 28], the authors study its convex relaxation problem

$$ \everymath{\displaystyle }\begin{array}{l@{\quad }l} \min_{L,S\in \mathcal{R}^{m\times n}} & \|L\|_*+ \lambda\|S\|_1 \\ \noalign {\vspace {4pt}} \mathrm{s.t.}& D = L+S, \end{array} $$
(2)

where ∥L∥∗ is the nuclear norm of L, and ∥S∥1=∑ i,j |S ij |. They give the exact recovery result, i.e., under certain assumptions, the optimal solution of (2) is exactly (L 0,S 0) with high probability. [28] carries on to investigate the noise case where D can be decomposed as

$$ D = L_0 + S_0 + N_0, \quad \|N_0\|\leqslant \delta, $$
(3)

with noise matrix N 0 and its corresponding convex relaxation problem

$$ \everymath{\displaystyle }\begin{array}{l@{\quad }l} \min_{L,S\in \mathcal{R}^{m\times n}} & \|L\|_*+ \lambda\|S\|_1 \\ \noalign {\vspace {4pt}} \mathrm{s.t.}& \|D - L -S\|\leqslant \delta. \end{array} $$
(4)

Similar recovery results are given there.

With the theoretical support, the main concern for RPCA is how to solve the resulting convex relaxation problems (2) and (4). Several numerical algorithms are proposed, including the Accelerated Proximal Gradient method (APG) [11, 21], the Approximate Feasible Direction method based on the dual problem (AFD) [11], and the Alternating Direction Method (ADM) [12, 16, 26]. For ADM, several solvers are developed, including IALM (Inexact Augmented Lagrangian Multiplier Method) in [12], LRSD (Low-Rank and Sparse matrix Decomposition) in [26] and LMaFit (Low-rank MAtrix Fitting) in [16]. They are widely used to solve (2) or (4) in application problems, such as video surveillance, face recognition, latent semantic indexing, ranking and collaborative filtering, robust photometric stereo, image alignment, and so on [13]. A common feature of such applications is that the data matrix D is a real matrix. However, as we will mention later, data in the moving target indication problem are complex, which is one of the differences between the GMTI problem and other application problems. Another difference is that in many real application problems, the assumptions for exact recovery are satisfied, whereas in the GMTI problem, some assumptions fail. It is these differences that make the GMTI problem special, and motivate us to investigate the problem.

Moving target indication, known as MTI, is the main aim in radar system. It covers two important classes of problems: ground moving target indication (GMTI) and air moving target indication (AMTI). Our interest here is the GMTI problem. As one of the radar systems, wide-area surveillance radar system is to detect moving targets on the ground from a wide area. Several radar systems already have such feature, including PAMIR (Phased Array Multifunctional Imaging Radar) and Joint STARS (Joint Surveillance Target Attack Radar System), see [1, 3–6, 8] for more details. Essentially, the wide-area surveillance radar system works as follows. Through antenna scanning, moving targets can be detected at different times and from different aspect angles, which makes the detection of moving targets possible, and one can further track the detected targets. The GMTI problem is therefore to extract some useful information, including the location of the moving target as well as its velocity, from the echoes received by antennas. To that end, the echoes, which consist of the information of not only moving targets, but also ground clutter and noises, must be carefully analyzed. In other words, to identify moving targets, ground clutter must be suppressed efficiently. Among the existing methods for ground clutter suppression, STAP (Space-Time Adaptive Processing) is one of the most efficient, which is basically a two-dimensional filter [6, 20]. However, it suffers from a high computational burden, and the performance is not good under the nonhomogeneous circumstances, mainly due to the mis-estimation of the interference covariance matrix. In [24], the authors first introduced robust PCA to model the GMTI problem. Noting that data of echo information are complex matrices, the authors use APG [11, 21] to solve the convex relaxation problem (4) for the real part and the imaginary part respectively, and unify the output matrices of APG into complex matrices as the final results. Simulation results demonstrate the advantages over STAP method. However, as mentioned in [24], the assumptions for the exact recovery of (4) are not satisfied due to the special structure of the sparse matrix (information matrix for moving targets), leaving the exact recovery of (4) a question. Due to the same reason, the RPCA model can not fully represent the GMTI problem.

Motivated by the above observations, in this paper, we will give some further investigation of the GMTI problem. We find that the sparse matrix in the GMTI problem actually follows a row-sparsity pattern (also referred as joint sparsity or group sparsity in compressive sensing). Based on such observation, we propose a so-called ‘structured RPCA’ model, which will better fit the GMTI problem, compared to the generic RPCA model. We modify IALM [12], one of the ADM solvers, to solve the corresponding convex relaxation problem. Further, to include more information of the structure of the GMTI problem, we also propose the row-modulus RPCA model to extract the moving target information. To solve the convex relaxation of the model, we give a closed-form solution for the resulting subproblem. Simulation results confirm the improvement of the proposed models.

The main contribution of our paper is twofold. Firstly, we propose the structured RPCA model to describe the problem, and discuss numerical methods to solve the convex relaxation problem. Secondly, we explore the row-modulus feature (will be defined later) to establish the so-called row-modulus RPCA model, which gives better extraction results. We discuss the numerical method to solve it, and a closed-form solution is derived for the subproblem.

The rest of the paper is organized as follows. In Sect. 2, we give a detailed description of the GMTI problem, and analyze the mathematical properties. In Sect. 3, we propose the structured RPCA model. Two methods are discussed to solve the convex relaxation of the model. In Sect. 4, by exploring the row-modulus structure, we propose the row-modulus RPCA model and discuss the numerical algorithm. In particular, a closed-form solution is derived for the resulting subproblem. Simulation results are shown in Sect. 5, with a comparison to the RPCA model in [24] to verify the efficiency of the two proposed models. It can be seen that the two models have great advantages over the one in [24]. In particular, the row-modulus RPCA model is overall better than the others. Final conclusions and some future questions are included in Sect. 5.

Some notations are given as follows.

For \(a\in \mathcal{C}\), a can be written as \(a = \mathit{Re}(a) + \mathfrak{j} \cdot \mathit{Im}(a)\), where \(\mathfrak{j}\) is the imaginary unit, Re(a) and Im(a) are the real part and the imaginary part of a respectively. In addition, a can also be written as \(a = |a|\texttt{exp}(\mathfrak{j}\theta)\), where |a| is the modulus (amplitude) and θ is the phase. For \(A\in \mathcal{C}^{m\times n}\), |A| is the matrix by taking entrywise modulus of A. For \(x\in \mathcal{C}^{n}\), x H is the conjugate transpose of x, and x T is the transpose of x. The l 2 norm is defined by \(\|x\|_{2}:=\sqrt{x^{H}x}\). For \(A\in \mathcal{C}^{m\times n}\), ∥A∥1,2 is the sum of the l 2 norm of the rows of A, i.e., \(\|A\|_{1,2}=\sum_{i = 1}^{m}\|A_{i\cdot}\|_{2} \), with A i⋅ is the i-th row of A. \(\|A\|_{F}:=\sqrt{\mathop {\mathit {tr}}(AA^{H})}\), and  tr(A) is the trace of A. ∥X∥∗ is the nuclear norm, which is the sum of the singular values of \(X\in\mathcal{C}^{m\times n}\). The vectorization operator \(\operatorname {vec}(\cdot): \mathcal{C}^{p\times q}\to \mathcal{C}^{m}\), with m=pq is defined as the column vector formed by stacking each column of \(X\in\mathcal{C}^{p\times q}\) one after another. For \(x\in \mathcal{C}^{n}\), the support of x is defined as Ω(x):={i:|x i |≠0}. Without notification, we use ∥⋅∥ as the Frobenius norm for matrices, and l 2 norm for vectors.

2 Ground Moving Targets Indication Problem

In this section, we will detail the GMTI problem in wide-area surveillance radar system, and give the mathematical reformulation of the problem.

The multi-channel wide-area surveillance radar system works as follows. As shown in Fig. 1, the airborne radar system equipped with n channels flies at the velocity of V along direction Y. The echo data are preprocessed per scan and per antenna look direction, respectively. At time t, data collected by n channels are first pre-processed through a series of stages, then transformed to the range-Doppler domain for further processing. We refer to [23, 24] for more details of data pre-processing.

Fig. 1
figure 1

Mechanism of GMTI

Full size image

Our main concern is to detect ground moving targets based on data in range-Doppler domain. For convenience, let \(\overline{Z}^{k}\in\mathcal{C}^{p\times q}\) be the data from the k-th channel after pre-processing with \(\overline{Z}^{k}_{ij}\) the echo data of the i-th Doppler cell and the j-th range cell. \(\overline{Z}^{k}_{ij}\) is composed of the stationary target \(\overline{C}^{k}_{ij}\) and the moving target \(\overline{T}^{k}_{ij}\), as well as the independent noise \(\overline{N}^{k}_{ij}\), i.e.,

$$ \overline{Z}^k_{ij} = \overline{C}^k_{ij} + \overline{T}^k_{ij} + \overline{N}^k_{ij}, \quad k = 1,\cdots, n, $$
(5)

where

$$ \overline{C}^k_{ij} = A_{ij}^c \hbox{exp}\bigl(\mathfrak{j}\phi_{ij}^k\bigr) $$

is the ground clutter component, indicating the information of the stationary target with \(A_{ij}^{c}\) the amplitude and \(\phi_{ij}^{k}\) the phase, and

$$ \overline{T}^k_{ij} = A_{ij}^t \hbox{exp}\bigl(\mathfrak{j} \overline{\psi}_{ij}^k\bigr) $$

is the moving target component, indicating the information of the moving target with \(A_{ij}^{t}\) the amplitude and \(\overline{\psi}_{ij}^{k}\) the phase. \(\overline{N}^{k}_{ij}\) is the noise for the corresponding cell. Here we would like to highlight that the moving target matrix \(\overline{T}^{k}\) is sparse, and the nonzero entries in \(\overline{T}^{k}\) are randomly distributed.

For each (i,j) cell, there is a constant phase difference due to the space between neighboring channels, that is,

$$ \phi_{ij}^k-\phi_{ij}^{k-1} \ \hbox{is constant}, \quad \hbox{and}\quad \overline{\psi}_{ij}^k-\overline{\psi}_{ij}^{k-1} \ \hbox{is constant}, \quad k = 2,\cdots, n. $$

In particular, for the uniformly spaced channel system, where the distance between each two neighboring channels are the same, there are

$$ \phi_{ij}^k-\phi_{ij}^{k-1} =c_{ij}, \quad \hbox{and}\quad \overline{\psi}_{ij}^k-\overline{\psi}_{ij}^{k-1}=\bar{d}_{ij}, \quad k = 2,\cdots,n, $$

for some nonzero constants c ij , and \(\bar{d}_{ij}\), i=1,⋯,p, j=1,⋯,q. In the rest of the paper, without loss of generality, we will always discuss the uniformly spaced channel system.

For the phase difference in \(\overline{C}^{k}_{ij}\), one can use the adaptive 2-D calibration method [7, 10] for channel removing phases as well as channel equalization. However, the phase difference in \(\overline{T}_{ij}^{k}\) can not be removed, due to the velocity of the moving target. Therefore, after removing the phase difference in \(\overline{C}^{k}_{ij}\), we get the following data:

$$ C^k_{ij} = C^{k-1}_{ij} , \quad k = 2,\cdots, n, $$
(6)
$$ Z^k_{ij} = C^k_{ij} +T^k_{ij} + N^k_{ij}, \quad k = 1,\cdots, n, $$
(7)

where

$$ T^k_{ij} = A_{ij}^t \hbox{exp}\bigl(\mathfrak{j} \psi_{ij}^k\bigr), $$

and the phase differences are

$$ \psi_{ij}^k- \psi_{ij}^{k-1}=d_{ij}\neq 0, \quad k = 2,\cdots, n. $$

Property 2.1

Note that for k=1,⋯,n, C k ’s are identical to each other, and T k ’s are sparse matrices, sharing the same support.

If we vectorize Z k, k=1,⋯,n, and put them as columns of a new matrix \(G\in \mathcal{C}^{m\times n}\) with m=pq, there is

$$ G = G_C + G_T+ G_N, $$
(8)

where

$$\begin{array}{c} G_C = \bigl[\operatorname {vec}\bigl(C^1\bigr),\cdots, \operatorname {vec}\bigl(C^n\bigr)\bigr], \qquad G_T =\bigl[\operatorname {vec}\bigl(T^1\bigr),\cdots, \operatorname {vec}\bigl(T^n\bigr)\bigr], \\ \noalign {\vspace {8pt}} G_N = \bigl[\operatorname {vec}\bigl(N^1\bigr),\cdots, \operatorname {vec}\bigl(N^n\bigr)\bigr]. \end{array} $$

We summarize the data procession in Fig. 2.

Fig. 2
figure 2

Data processing

Full size image

Property 2.2

Based on Property 2.1, it is not difficult to get that

  1. (i)

    In the ideal case, G C is a rank one matrix, with identical columns.

  2. (ii)

    G T is a sparse matrix, with few rows nonzeros, and most rows zeros.

  3. (iii)

    Besides (ii), the amplitude (modulus) matrix |G T |, is also a rank one matrix, since for each nonzero row i, the entries |G T | ij , j=1,⋯,n, are the same.

  4. (iv)

    For each nonzero row i in the phase matrix Ψ of G T , the differences of each neighboring pair are the same, i.e.,

    $$ \varPsi_{ij} -\varPsi_{ij-1} =d_i, \quad j = 2,\cdots, n. $$
    (9)

Now we are ready to give the mathematical reformulation of the GMTI problem.

GMTI Problem

Given matrix \(G\in \mathcal{C}^{m\times n}\), find a low rank matrix \(L_{0}\in\mathcal{C} ^{m\times n}\) and a sparse matrix \(S_{0}\in\mathcal{C}^{m\times n}\) such that

  1. (a)

    G=L 0+S 0+N 0, ∥N 0∥⩽δ,

  2. (b)

    L 0 and S 0 have properties in Property 2.2.

The rest of the paper is devoted to the investigation of the GMTI problem from the perspective of matrix decomposition. We will mainly focus on the mathematical models and the numerical algorithms. Before proposing our models, we first give a brief review of the RPCA model in [24].

2.1 RPCA Model

With Property 2.2 (i), the RPCA model in [24] is to solve the convex relaxation problem (4) with the real part and the imaginary part as the input matrices, respectively, then form the output matrices into complex matrices as the final results. The outline of the method is given as follows.

Algorithm 2.1

RPCA Model (rpca)

Step 1:

Solve (4) with input D:=Re(G) to get the optimal solution (\(L_{R}^{*},S^{*}_{R}\)).

Step 2:

Solve (4) with input D:=Im(G) to get the optimal solution (\(L_{I}^{*},S^{*}_{I}\)).

Step 3:

Let \(G_{C}^{*}= L_{R}^{*} + \mathfrak{j}L_{I}^{*}\), \(G_{T}^{*}= S_{R}^{*} + \mathfrak{j}S_{I}^{*}\). Output \((G_{C}^{*}, G_{T}^{*})\).

To give the recovery results for (4), we need the following preparation. Let L 0∈ℜm×n with rank r have the singular value decomposition \(L_{0} = U\varSigma V^{H}=\sum_{i = 1}^{r}\sigma_{i}u_{i}v_{i}^{H}\), where \(\varSigma = \operatorname {Diag}(\sigma_{1}, \cdots, \sigma_{r})\) is the diagonal matrix with nonzero singular values σ 1⩾⋯⩾σ r , and U∈ℜm×r, V∈ℜn×r are formed by the corresponding left and right singular vectors. Assume the following conditions hold.

  1. (A1)

    Incoherence condition with parameter μ:

    $$ \max_i \bigl\|U^He_i \bigr\| \leqslant \sqrt{\frac{\mu r}{m}}, \qquad \max_i \bigl\|V^He_i\bigr\| \leqslant \sqrt{\frac{\mu r}{n}}, \qquad \bigl\|UV^H\bigr\|_{\infty} \leqslant \sqrt{ \frac{\mu r}{mn}}. $$
    (10)
  2. (A2)

    The sparsity pattern of S 0 is selected uniformly at random.

The exact recovery result for the convex relaxation problem (2) is given as follows [2].

Theorem 2.1

Suppose \(L_{0} \in \mathcal{R}^{n\times n}\) obeys (10), and that the support set of S 0 is uniformly distributed among all sets of cardinality m. Then there is a numerical constant c such that with probability at least 1−cn −10(over the choice of support of S 0), principle component pursuit (2) with \(\lambda = \frac{1}{\sqrt{n}}\) is exact, i.e., \(\hat{L} = L_{0}\), \(\hat{S} = S_{0}\), where \((\hat{L},\hat{S})\) is the optimal solution of (RPCA), provided that

$$ \operatorname {rank}(L_0)\leqslant \rho_r n \mu^{-1} (\log n )^{-2} \quad\hbox{and}\quad m\leqslant {\rho_s} n^2. $$
(11)

Above, ρ r and ρ s are positive numerical constants. In the general rectangular case where L 0 is m 0×n 0, (2) with \(\lambda = \frac{1}{\sqrt{n_{1}}}\) succeeds with probability at least \(1-cn_{1}^{-10}\) provided that rank(L 0)⩽ρ r n 2 μ −1(logn 1)−2 and m⩽ρ s n 1 n 2, where n 1=max(m 0,n 0) and n 2=min(m 0,n 0).

Remark 2.1

If Re(G C ), Im(G C ) satisfy assumption (A1), and Re(G T ), Im(G T ) satisfy assumption (A2), Algorithm 2.1 can return G C and G T with high probability. However, in the GMTI problem, as mentioned in Property 2.2(ii), the sparse support in Re(G T ) and Im(G T ) are obviously not uniformly distributed, violating assumption (A2). It leaves the exact recovery a question. It is this observation that motivates the following sections, which are also the main contributions of this paper.

3 Structured RPCA Model for the GMTI Problem

In this section, we propose the structured RPCA model to better fit the problem. We also modify IALM [12], one of the ADM solvers, to solve the convex relaxation of the problem. Finally, another alternative model aiming at dealing with complex matrices is also discussed.

3.1 The Structured RPCA Model

In this subsection, we will propose the structured RPCA model, and discuss the numerical algorithm to solve the convex relaxation problem.

Here and in the following, we say that \(X\in\mathcal{C}^{m\times n}\) is sparse in the sense of rows by meaning that the vector x:=[∥X 1⋅∥2,⋯,∥X m⋅∥2]T∈ℜm is sparse. We mean the row-support of the row-sparse matrix X by meaning the support of x, i.e., Ω(X):={i: ∥X i⋅∥2≠0}.

Based on Property 2.2(ii), we know that the sparsity of G T is actually in the sense of rows. Such structure is often referred as joint-sparsity in compressive sensing. Interested readers are referred to [14, 17–19] for further theories and algorithms. Therefore, we can take each row as a block, and consider the row-sparsity of G T . This results the following model, referred as the structured RPCA model.

Structured RPCA Model

Given D∈ℜm×n, find a low rank matrix L 0 and a row-sparse matrix S 0 such that

$$ D = L_0 + S_0 + N_0, \quad \|N_0\|\leqslant \delta. $$
(12)

Here ‘structure’ means that S 0 has the row-sparsity structure. Obviously, the structured RPCA model can better fit the GMTI problem, compared with the generic RPCA model. It is expected to obtain better simulation results as long as the resulting convex relaxation problem is efficiently solved. Similar to (4), the corresponding convex relaxation optimization problem takes the following form:

$$ \everymath{\displaystyle }\begin{array}{l@{\quad }l} \min_{L,S\in \mathcal{R}^{m\times n}} & \|L\|_*+ \lambda\|S\|_{1,2} \\ \noalign {\vspace {4pt}} \mathrm{s.t.}& \|D - L -S\|\leqslant \delta. \end{array} $$
(13)

The method based on the structured RPCA model is outlined as follows.

Algorithm 3.1

Structured RPCA model (srpca)

Step 1:

Solve (13) with input D:=Re(G) to get the optimal solution (\(L_{R}^{*},S^{*}_{R}\)).

Step 2:

Solve (13) with input D:=Im(G) to get the optimal solution (\(L_{I}^{*},S^{*}_{I}\)).

Step 3:

Let \(G_{C}^{*}= L_{R}^{*} + \mathfrak{j}L_{I}^{*}\), \(G_{T}^{*}= S_{R}^{*} + \mathfrak{j}S_{I}^{*}\), output \((G_{C}^{*}, G_{T}^{*})\).

Remark 3.1

Note that the row-support for the row-sparse matrix G T is randomly distributed. However, it is pointed out in [9] that, exact recovery is not possible if nonzero entries are across entire row or column. Despite of that, the promising simulation results still confirm the improvement of the structured RPCA model over the RPCA model.

Next, we will discuss how to solve the convex relaxation optimization problem (13). As mentioned in Introduction, existing methods for (2) and (4) include APG [11, 21], AFD [11], and ADM [12, 16, 26]. In [24], the authors choose APG to solve (4). However, in [22], the author conducts systematic comparison and concludes that ADM (with solver IALM [12]) performs overall better than APG and AFD. In fact, the idea of ADM has been widely used in many algorithms, such as [16, 27], and so on. Therefore, we choose the IALM solver in [12] to solve (13).

To apply IALM, we have to first derive the extended version of shrinkage operator, which is achieved based on the Moreau decomposition theorem [15, Theorem 31.5]. A proof can also be found in [25, Lemma 3.4]. Given z∈ℜp, and ϵ>0, consider the Moreau-Yosida regularization:

$$ \min_{x\in \Re^p } \epsilon \| x \|_2 + \frac{1}{2}\|x- z\|_2^2. $$
(14)

The solution \(\mathcal{T}_{\epsilon}(z)\) takes the following form:

$$ \mathcal{T}_{\epsilon}(z)= \left \{ \begin{array}{l@{\quad }l} 0,& \|z\|_2\leqslant \epsilon,\\ \noalign {\vspace {2pt}} z-\epsilon\frac{z}{\|z\|_2}, &\hbox{otherwise}. \end{array} \right . $$
(15)

In particular, if p=1, \(\mathcal{T}_{\epsilon}(z)\) reduces to the following the shrinkage operator:

$$ \mathcal{S}_{\epsilon}(z) = \left \{ \begin{array}{l@{\quad }l} z-\epsilon, & z >\epsilon,\\ \noalign {\vspace {2pt}} z+\epsilon, & z<-\epsilon,\\ \noalign {\vspace {2pt}} 0,& \hbox{otherwise}. \end{array} \right . $$
(16)

For vectors and matrices, \(\mathcal{S}_{\epsilon}(\cdot)\) is performed in an elementwise way.

Remark 3.2

For the complex case, i.e., \(z\in \mathcal{C}^{p}\) and \(x\in\mathcal{C}^{p}\) in (14), \(\mathcal{T}_{\epsilon}(z)\) is the same as in (15).

Now we are ready to give the idea of the IALM solver for (13). Recall the augmented Lagrangian function of (13) is given by:

$$ \mathcal{L}_{\mu} (L,S;Z)= \|L\|_*+\lambda\|S \|_{1,2}-\langle Z, D-L-S\rangle + \frac{\mu}{2}\|D-L-S\|_F^2. $$
(17)

For more details of IALM, please refer to [12] and [22].

Algorithm 3.2

IALM for (13)

Input: Observation matrix D∈ℜm×n, λ>0.

  1. 1

    Z 0=D/J(D) with J(D):=max(∥D∥2,λ −1∥D∥∞); S 0=0; μ 0>0; ρ>1; k=0.

  2. 2

    While not converged do

  3. 3

    Solve \(L^{k+1} = \arg\min_{L\in \Re^{m\times n}} \mathcal{L}_{\mu_{k}} (L,S^{k};Z^{k})\) by

    $$ L^{k+1} = U\mathcal{S}_{\mu_k^{-1}}( \varSigma)V^T, \quad \hbox{where} \ D-S^k+\mu_k^{-1}Z^k = U\varSigma V^T. $$
    (18)
  4. 4

    Solve \(S^{k+1} = \arg\min_{S\in \Re^{m\times n}} \mathcal{L}_{\mu_{k}} (L^{k+1},S;Z^{k} )\) by (19).

  5. 5

    Z k+1=Z k+μ k (D−L k+1−S k+1); μ k+1=ρμ k . k:=k+1.

  6. 6

    end while

  7. 7

    Output: (L k,S k).

Note that the only difference between IALM in our paper and that in [12] is Step 4. In our IALM solver, S k+1 is calculated row by row using the extended shrinkage operator. Specifically,

For \(S_{i\cdot}^{k+1}\), i=1,⋯,m, there is

(19)

3.2 Another Alternative for Structured RPCA Model

We also propose another alternative to realize the structured RPCA model, referred as the structured RPCA-c model, which deals with complex matrices directly.

Structured RPCA-c Model

Given matrix \(D\in\mathcal{C}^{m\times n}\), find a low rank matrix L 0 and a row-sparse matrix S 0 such that

$$ D = L_0 + S_0 + N_0, \quad \|N_0\|\leqslant \delta. $$
(20)

The corresponding convex relaxation problem is

$$ \everymath{\displaystyle }\begin{array}{l@{\quad }l} \min_{L,S\in \mathcal{C}^{m\times n}} & \|L\|_*+ \lambda\|S\|_{1,2} \\ \noalign {\vspace {3pt}} \mathrm{s.t.}& \|D - L -S\|\leqslant \delta. \end{array} $$
(21)

For \(X, Y \in \mathcal{C}^{m\times n}\), define 〈X,Y〉:=Re(tr(X H Y)). We can directly apply Algorithm 3.2 to solve (21) with some slight modification in Step 4 caused by the complex numbers. We name the resulting method as srpca-c.

4 Row-Modulus RPCA Model

In this section, we will propose a new model, referred as the row-modulus RPCA model, by taking into account of (i)–(iii) in Property 2.2. To make use of (iii) in Property 2.2, the entries in the same row must have the same modulus, therefore, we can describe the model as follows.

Row-Modulus RPCA Model

Given \(D\in\mathcal{C}^{m\times n}\), find a low rank matrix L 0 and a row-sparse matrix S 0 such that

$$ \begin{array}{l} D= L_0 + S_0 + N_0, \quad \|N_0\|\leqslant \delta, \\ \noalign {\vspace {5pt}} \bigl|(S_0)_{ij}\bigr| = \bigl|(S_0)_{ik}\bigr|, \quad j,\ k \in\{1,\cdots, n\},\ \forall i, \end{array} $$
(22)

where (22) means that for the entries in the same row, they have the same modulus. The corresponding convex relaxation problem is

$$ \everymath{\displaystyle }\begin{array}{l@{\quad }l} \min_{L,S\in \mathcal{C}^{m\times n}} & \|L\|_*+ \lambda\|S\|_{1,2}\\ \noalign {\vspace {3pt}} \mathrm{s.t.}& \|D - L -S\|\leqslant \delta,\\ \noalign {\vspace {3pt}} &|S_{ij}| = |S_{ik}|, \quad j, \ k \in\{1,\cdots, n\}, \ \forall i. \end{array} $$
(23)

By introducing vector s∈ℜm, we have the following equivalent form:

$$ \everymath{\displaystyle }\begin{array}{l@{\quad }l} \min_{L,S\in \mathcal{C}^{m\times n},s\in\Re^m} & \|L\|_*+ \lambda\|S\|_{1,2} \\ \noalign {\vspace {3pt}} \mathrm{s.t.}& \|D - L -S\|\leqslant \delta,\\ \noalign {\vspace {3pt}} &|S_{ij}| = s_i, \quad j\in\{1,\cdots, n\}, \ i = 1,\cdots, m,\\ \noalign {\vspace {3pt}} &s_i\geqslant 0, \quad i = 1,\cdots, m. \end{array} $$
(24)

Define \(\mathcal{D}:=\{S\in \mathcal{C}^{m\times n}: \ |S_{ij}| = s_{i}, \ s_{i}\geqslant 0, \ j =1,\cdots, n, \ i = 1,\cdots, m \}\), and \(\mathcal{D}_{i}:=\{x\in \mathcal{C}^{ n}: \ |x_{j}| = s_{i},\ s_{i}\geqslant 0,\ j =1,\cdots, n\}\), i=1,⋯, m.

Next, we will discuss the numerical algorithm for solving (24). We still follow the frame of IALM. The augmented Lagrangian function is the same as in (17). During each iteration, L k+1 is calculated as in (18). The only difference is in S k+1, where one needs to solve the following problem:

That is,

(25)

Next, we will focus on the derivation of \(S^{k+1}_{i\cdot}\). Given \(y\in\mathcal{C}^{n}\) and ϵ>0, consider the following problem:

$$ \everymath{\displaystyle }\begin{array}{l@{\quad }l} \min_{x\in \mathcal{C}^n, a\in\Re } & \epsilon \|x\|_2 + \frac{1}{2}\|x- y\|_2^2\\ \noalign {\vspace {4pt}} \mathrm{s.t.}& |x_j| = a, \quad j = 1,\cdots, n,\\ & a\geqslant 0. \end{array} $$
(26)

We will show that (26) has a closed-form solution.

Note that for any feasible point x of (26), \(\|x\|_{2} = \sqrt{n} a \). Then (26) is equivalent to the following problem

$$ \everymath{\displaystyle }\begin{array}{l@{\quad }l} \min_{ x\in \mathcal{C}^n, a\in \Re} & f(a, x):= \frac{1}{2} na^2 -a \bigl(\langle x, y\rangle-\epsilon \sqrt{n} \bigr) \\ \noalign {\vspace {5pt}} \mathrm{s.t.}& |x_j| = 1, \quad j = 1,\cdots, n,\\ \noalign {\vspace {5pt}} & a\geqslant 0. \end{array} $$
(27)

f(a,x) can be further reformulated as

$$ f(a, x):= \frac{1}{2} n \biggl(a -\frac{1}{n}\bigl(\langle x, y\rangle-\epsilon \sqrt{n}\bigr)\biggr)^2 -\frac{1}{2n}\bigl(\langle x, y\rangle-\epsilon\sqrt{n}\bigr)^2. $$
(28)

The following proposition gives the optimal solution of (27).

Proposition 4.1

The optimal solution (a ∗,x ∗) of (27) is given by

$$ \left \{ \begin{array}{l@{\quad }l} a^*=0,\ x^* \ \hbox{is any vector with unit entries}, & \hbox{if} \ \|y\|_1- \epsilon \sqrt{n}\leqslant 0, \\ \noalign {\vspace {4pt}} a^* = \frac{1}{n}(\|y\|_1 - \epsilon \sqrt{n}),\ x^*= \hat{x} , &\hbox{otherwise}, \end{array} \right . $$
(29)

where \(\hat{x}\) is the entrywise division \(\frac{y}{|y|}\) with the convention \(\frac{0}{0}=0\).

Proof

By \(\langle \hat{x} ,y\rangle = \|y\|_{1}\) and \(-\langle \hat{x} ,y\rangle = -\|y\|_{1}\), we have 〈x,y〉∈[−∥y∥1,∥y∥1], i.e., \(\langle x ,y\rangle - \epsilon \sqrt{n} \in[-\|y\|_{1} - \epsilon \sqrt{n},\|y\|_{1}- \epsilon \sqrt{n}]\). In order to discuss the optimal solution, we consider the following two cases.

Case 1. \(\|y\|_{1}- \epsilon \sqrt{n}\leqslant 0\). It is trivial that the optimal solution is \(a^{*} =\allowbreak \max\{0, \frac{1}{n}(\langle x ,y\rangle- \epsilon \sqrt{n})\}=0 \), x ∗ is any vector satisfying

$$ |x_j| = 1, \ j = 1,\cdots, n. $$
(30)

The optimal value is f(a ∗,x ∗)=0.

Case 2. \(\|y\|_{1}- \epsilon \sqrt{n} > 0\). If x is chosen such that \(\langle x ,y\rangle - \epsilon \sqrt{n} \in [-\|y\|_{1} - \epsilon \sqrt{n},0 ]\), then a=0 and f(a,x)=0. Otherwise, if x is chosen such that \(\langle x ,y\rangle - \epsilon \sqrt{n} \in[0,\|y\|_{1} - \epsilon \sqrt{n}]\), then \(a = \frac{1}{n}(\langle x ,y\rangle - \epsilon \sqrt{n}) \), resulting the first part of f zero. Note that \(\langle \hat{x} ,y\rangle - \epsilon \sqrt{n}>0\). In such situation, to minimize the second part of f, x ∗ has to be \(\hat{x} \), and it follows that \(a^{*} =\frac{1}{n}(\langle \hat{x} ,y\rangle - \epsilon \sqrt{n}) = \frac{1}{n}(\|y\|_{1} - \epsilon \sqrt{n}) \). Therefore, (a ∗,x ∗) with \(a^{*} = \frac{1}{n}(\|y\|_{1} - \epsilon \sqrt{n}) \) and \(x^{*}= \hat{x} \) is the optimal solution of (27).

The proof is finished. □

Remark 4.1

With Proposition 4.1, the optimal solution of (26) can be given by (a ∗,a ∗ x ∗). Consequently, \(S^{k+1}_{i\cdot}\) in (25) can be efficiently calculated.

Now we are ready to give the method for (24), which is named as rm-rpca.

Algorithm 4.1

Row-Modulus RPCA model (rm-rpca)

Input: Observation matrix \(D\in\mathcal{C}^{m\times n}\), λ>0.

  1. 1

    Z 0=D/J(D) with J(D):=max(∥D∥2,λ −1∥D∥∞); S 0=0; μ 0>0; ρ>1; k=0.

  2. 2

    While not converged do

  3. 3

    Solve \(L^{k+1} = \arg\min_{L\in \mathcal{C}^{m\times n}} \mathcal{L}_{\mu_{k}} (L,S^{k};Z^{k})\) by

    $$ L^{k+1} = U\mathcal{S}_{\mu_k^{-1}}(\varSigma)V^H, \quad \hbox{where} \ D-S^k+\mu_k^{-1}Z^k = U\varSigma V^H. $$
    (31)
  4. 4

    Solve \(S^{k+1} = \arg\min_{S\in \mathcal{D}}\mathcal{L}_{\mu_{k}} (L^{k+1},S;Z^{k} )\) by (25) and Remark 4.1.

  5. 5

    Z k+1=Z k+μ k (D−L k+1−S k+1); μ k+1=ρμ k . k:=k+1.

  6. 6

    end while

  7. 7

    Output: (L k,S k).

5 Simulation Results

In this section, we conduct the simulation to verify the performance of three numerical algorithms, including srpca, srpca-c, rm-srpca for three proposed models: the structured RPCA model, the structured RPCA-c model and the row-modulus RPCA model. We also compare with rpca, the numerical algorithm for the RPCA model in [24]. We choose IALM in [12], one of the ADM solvers, to solve the convex relaxation problem (4) in rpca. The simulation is performed in Matlab 2010a on a computer with Pentium(R) Core(TM) i3-2350 CPU 2.30 GHz and RAM 3.41GB.

We analyze data from a five-channel wide-area surveillance radar system, which are the same as used in [24]. We simulate the radar echoes from one antenna look direction (side-looking) and extract the data from 1000th to 1299th range cell to do the experiment. For the linear frequency modulation (LFM) signal, the overall received echoes from the entire scene are assumed to be the superposition of the returns from all point targets including stationary targets and moving targets. 128 stationary points are disposed in each range cell with the distribution of Gaussian and Clutter to Noise Ratio (CNR) is about 33 dB. In these 300 range cells, 6 moving targets are disposed at the 1045th, 1075th, 1105th, 1195th, 1225th, 1255th range cells with the radical velocities 3.8 m/s, 3 m/s, 2.3 m/s, −2.3 m/s, −3 m/s, −3.8 m/s respectively, and the Signal to Noise Ratio (SNR) is set to 20 dB. In other words, n=5, \(T^{k}\in\mathcal{C}^{128\times 300}\), and \(G\in\mathcal{C}^{38400\times 5}\). Figure 3 shows the modulus of the observation data by the third channel in the range-Doppler domain before clutter suppression, where one can not identify the six moving targets. Figure 4 shows the modulus of the true echoes of the moving targets for the third channel.

Fig. 3
figure 3

Observation data matrix in the range-Doppler domain of the third channel before clutter suppression

Full size image
Fig. 4
figure 4

True moving target data matrix in the range-Doppler domain of the third channel

Full size image

Parameter Setting

We first show some preliminary results on parameters μ and λ since they are closely related to the performance of the algorithms. For μ in each algorithm, we mainly consider the following two update rules:

  1. R1.

    μ k+1=ρ 1 μ k (default update rule in IALM [12]);

  2. R2.

    \(\mu_{k+1} =\left \{\begin{array}{l@{\quad }l} \rho_{2}\mu_{k}, &\hbox{ if }\|D-L^{k+1}-S^{k+1}\|<10^{2}\cdot\texttt{tol}, \\ \noalign {\vspace {3pt}} \rho_{3}\mu_{k} &\hbox{otherwise}, \end{array} \right .\)

where ρ 1=1.5, ρ 2=1.69, ρ 3=1.3, tol=10−7 is the default tolerance. In terms of λ, we pick up five values for each model. For rpca, λ 1,⋯,λ 5 are chosen as \(\frac{1}{\sqrt{100}}\), \(\frac{1}{\sqrt{200}}, \cdots,\frac{1}{\sqrt{500}}\), respectively, while for the rest, λ 1,⋯,λ 5 are chosen as \(\frac{1}{\sqrt{40}}\), \(\frac{1}{\sqrt{60}}, \cdots, \frac{1}{\sqrt{120}}\). The residuals are reported in Table 1, where \(\mathit{res} := \|G_{T}-G_{T}^{*}\|\), the residual between the true moving target matrix and the estimated one. We find that the update rule R2 leads to smaller residuals for all the algorithms. In terms of λ, the results demonstrate that smaller λ leads to smaller residuals. Therefore, we choose λ 5 in each algorithm.

Table 1 res for different parameters
Full size table

Simulation Results

Based on the above analysis, we choose R2 to update μ k for all the algorithms, \(\lambda =\frac{1}{\sqrt{500}}\) for rpca and \(\lambda =\frac{1}{\sqrt{120}}\) for others. Recall that the main aim of the GMTI problem is to extract information of the moving targets, including the location, the energy and the velocity. The above three can be represented by the support of the row-sparse matrix G T , the modulus of each nonzero row in G T and the phase difference between entries in each nonzero row in G T . Therefore, in the following, we will qualify the results of different algorithms from the above three aspects.

(a) Row-Support

As mentioned before, the row-support of the estimated sparse matrix \(G_{T}^{*}\) indicates the location of the moving targets. Typical simulation results of the four algorithms: rpca, srpca, srpca-c and rm-rpca are shown in Fig. 5. It demonstrates that the six moving targets can be successfully identified by each algorithm.

Fig. 5
figure 5

Moving targets extraction results: the third channel

Full size image

(b) Row-Modulus and Phase Difference

To compare row-modulus and phase difference, we do more comparisons on different residuals and report them in Table 2. Let Ω be the row-support of G T . \(\mathit{res}_{a}= \| |P_{\varOmega}(G_{T})|-|P_{\varOmega}(G_{T}^{*})|\|\), indicates the residual of modulus between the true moving targets and the estimated ones. \(\mathit{res}_{p}= \| P_{\varOmega}(\varPsi_{d}- \varPsi_{d}^{*})\|\), is the residual between the true phase difference matrix Ψ d and the estimated phase difference matrix \(\varPsi_{d}^{*}\), where Ψ d ∈ℜm×(n−1) is defined as (Ψ d ) ij =Ψ ij+1−Ψ ij , i=1,⋯,m, j=1,⋯,n−1, and similarly, \((\varPsi_{d}^{*})_{ij} = \varPsi^{*}_{ij+1}-\varPsi_{ij}^{*}\), i=1,⋯,m, j=1,⋯,n−1. One may notice that, for rpca, srpca and srpca-c, the modulus for the entries in the same row may be different. To make the comparison reasonable, one natural way is to take the averaged modulus across each nonzero row as the estimated modulus of the moving target. Based on such analysis, we compare the estimated averaged modulus with the true modulus. Denote g as the modulus vector of the true moving targets, and g ∗ the averaged estimated modulus of the moving targets. res aa =∥g−g ∗∥. Similarly, we can calculate the averaged phase difference, and compare with the true phase difference, i.e., res pa =∥ψ−ψ ∗∥, where ψ is the true phase difference vector, and ψ ∗ is the estimated averaged phase difference vector.

Table 2 Detailed residuals
Full size table

Based on the results in Table 2, we can find that: (i) compared to rpca, the decrease in res p and res pa implies some improvement in estimating phase difference for srpca and srpca-c, but no significant improvement in row-modulus; (ii) rm-rpca demonstrates great advantages in estimating both row-modulus and phase difference since it returns the smallest residuals among the four algorithms. Here we would like to highlight that the phase difference is crucial as it is closely related to the velocity of the moving targets, which is very important in radar system. In that sense, the two proposed models, especially the row-modulus RPCA model, can extract more accurate information of the moving targets.

6 Conclusions

In this paper, we discussed the mathematical model for the GMTI problem in multi-channel wide-area surveillance radar system. Based on the RPCA model in [24], we go a further step to propose two models: the structured RPCA model and the row-modulus RPCA model, which can better fit the GMTI problem by taking use of the row-sparsity information of the sparse matrix. Three algorithms: srpca, srpca-c and rm-rpca are discussed to solve the resulting convex relaxation optimization problems. Simulation results demonstrate the advantages of the proposed models. In particular, the row-modulus RPCA is the best compared to the structured RPCA model and the RPCA model. However, several questions still remain unsolved. One is that although we included (i)–(iii) in Property 2.2 in the row-modulus RPCA model, it is still far from satisfaction. Property 2.2(iv) is not fully explored to model the problem. Therefore, how to take full use of such information to extract moving target information is still a question. Another one is how to establish the equivalent result between the original problem and the relaxed convex optimization problem. It will provide the theoretical support to our investigation on numerical algorithms. Both of them are very interesting topics and worth further investigation.