Distributed Compressed Sensing Based Ground Moving Target Indication for Dual-Channel SAR System

Abstract

The dual-channel synthetic aperture radar (SAR) system is widely applied in the field of ground moving-target indication (GMTI). With the increase of the imaging resolution, the resulting substantial raw data samples increase the transmission and storage burden. We tackle the problem by adopting the joint sparsity model 1 (JSM-1) in distributed compressed sensing (DCS) to exploit the correlation between the two channels of the dual-channel SAR system. We propose a novel algorithm, namely the hierarchical variational Bayesian based distributed compressed sensing (HVB-DCS) algorithm for the JSM-1 model, which decouples the common component from the innovation components by applying variational Bayesian approximation. Using the proposed HVB-DCS algorithm in the dual-channel SAR based GMTI (SAR-GMTI) system, we can jointly reconstruct the dual-channel signals, and simultaneously detect the moving targets and stationary clutter, which enables sampling at a further lower rate in azimuth as well as improves the reconstruction accuracy. The simulation and experimental results show that the proposed HVB-DCS algorithm is capable of detecting multiple moving targets while suppressing the clutter at a much lower data rate in azimuth compared with the compressed sensing (CS) and range-Doppler (RD) algorithms.

1. Introduction

Synthetic aperture radar (SAR) is a kind of high resolution imaging radar, which is capable of working at long distances, all-weather, and day and night. In the SAR system, the scene is observed by an antenna at different positions in the azimuth direction. The coherent information recorded at the different positions is used to synthesize a very long antenna to improve the azimuth resolution. The dual-channel SAR system is widely applied in the field of ground moving-target indication (GMTI). In practice, the Doppler shift of the stationary clutter spreads the Doppler bandwidth of the clutter returns due to the motion of the platform. Accordingly, the Doppler bandwidth of the clutter will mask that of moving targets [1]. Dual-channel SAR can effectively suppress the clutter, which is propitious to slowly moving targets detection [2]. Some traditional approaches, such as the displaced phase center antenna (DPCA) approach, and along-track interferometry (ATI) detection [3], have achieved good results. However, relying on the Nyquist sampling theorem, dual-channel SAR is obliged to transmit and store substantial raw data samples, which limits its application in practice.

Compressed sensing (CS) [4,5,6,7], as a novel sparse reconstruction technique, is able to recover the signals that have been sampled below the traditional Nyquist sampling rate. It provides a good solution for the transmission and storage of substantial raw data samples. CS aims at minimizing the number of measurements and attempts to recover the original signal from these necessary measurements by solving an ${\ell }_0$-norm optimization problem, or equivalently a convex ${\ell }_1$-norm optimization problem under the restricted isometry property (RIP) condition [8].

A lot of valuable literature on CS-based SAR imaging and the GMTI system has been published in recent years. A CS-based SAR imaging algorithm is proposed in [9], which demonstrates that CS has the potential to eliminate the need for the pulse compression matched filter at the receiver, and reduce the required receiver analog-to-digital conversion bandwidth. In [10], CS is applied to randomly select lines in azimuth from the signals after range compression, which allows the implementation of wide-swath modes without reducing the azimuth resolution. Sun et al. [11] propose a CS-based method for joint sparse recovery of all channel or sub-aperture images. In [12], a fast CS-based SAR imaging algorithm is proposed to save the computational cost both in time and memory. In [13], a phase error correction method is proposed for the CS-based radar imaging system based on approximated observation, which improves the defocus caused by the phase error. In [14], CS is utilized in tomographic SAR system, extended to the SAR elevation direction for 3-D imaging. In the field of GMTI, in [15] CS is used to estimate the velocities and positions of moving targets. In [16,17], a CS-based space-time adaptive processing (STAP) algorithm is proposed to detect the moving targets. The above CS-based algorithms can efficiently reduce data rate and improve image quality, as well as improve anti-jamming performance. However, for the dual-channel SAR based GMTI (SAR-GMTI) system, with the increase of the imaging resolution, the resulting substantial raw data samples aggravate transmission and storage load. Thus it is required to develop an efficient algorithm that can further reduce the data rate.

In this work, we tackle the problem by exploiting the correlation between two channels of the dual-channel SAR system. The antennas of dual-channel SAR are configured in the along-track direction. Thus the SAR images acquired by each channel are highly correlated during a flight along the track. Using distributed compressed sensing (DCS) based algorithms, we can jointly reconstruct the dual-channel signals, which enables sampling at a further lower rate in azimuth as well as improves the reconstruction accuracy.

DCS [18] exploits both intra-signal and inter-signal correlations. Three different joint sparsity models have been proposed, i.e., joint sparsity model 1 (JSM-1), joint sparsity model 2 (JSM-2) and joint sparsity model 3 (JSM-3). In JSM-1, a sparse signal consists of a common component and an innovation component. Such signals may exist when large-scale phenomena affects all signal sources and local phenomena affects specific individual source. In JSM-2, all signals are constructed from the same sparse set of basis vectors, but with different coefficient values. JSM-3 extends JSM-1 so that the common component needs no longer to be sparse in any basis. In this work, the echo signal of dual-channel SAR is the combination of the signals from both stationary clutter and moving targets. After time delay and phase compensation, the signal from stationary clutter is common to all channels, while the signals from moving targets are specific and different across different channels. Thus the JSM-1 model is adopted, considering that the signal from stationary clutter can be treated as the common component, while those from moving targets can be treated as the innovation components.

Many reconstruction algorithms have been proposed for the JSM-1 model, either in a central manner or a distributed manner. The central recovery methods assume the presence of a fusion center where all information is gathered. The joint recovery strategy (JRS) [18,19] is to stack the common component and all innovation components together as one single sparse signal, and reconstruct the single sparse signal via ${\ell }_0$-norm minimization or ${\ell }_1$-norm minimization. There are well developed algorithms such as subspace pursuit (SP), basic pursuit (BP), iterative hard thresholding (IHT), etc., to solve the minimization problem. However, JRS is computationally expensive since the dimension of the sparse signal is multiplied by stacking. As another branch of the methods, the distributed recovery performs the reconstruction in the network, with no fusion center. The Texas Hold ’Em algorithm [20] separates the estimation of the common component and innovation components by the measurement decomposition process. Valsesia et al. [21] exploit the side information from one of the signals to reconstruct the innovation components. Recently Sundman et al. [22] propose the parallel pursuit with side information (SIPP) algorithm, which extends SP by using the side information of common support set at each iteration. SIPP can reconstruct the signals in a distributed manner, thus achieving much less computing time than JRS. This implies that the information of support set can be utilized to accelerate the reconstruction process. In [23], an alternating direction method of multipliers (ADMM) based algorithm is proposed for the JSM-1 model in a decentralized framework, which exploits the local cooperation between nodes.

In this work, we consider dealing with the JSM-1 model in a different way, i.e., in a Bayesian framework. Bayesian methods essentially approximate the posterior distributions of sparse signals according to the prior knowledge and data, which is effective in dealing with uncertain models and large noises. We consider using the variational Bayesian algorithm for the JSM-1 model. Variational Bayesian (VB) inference [24] obtains the optimal solution by iterating over a set of mutually dependent equations, providing a good approximation to the exact posterior with a readily monitored convergence.

Variational Bayesian has been widely applied in sparse signal reconstruction [25,26]. In this work, we propose a novel algorithm, namely the hierarchical variational Bayesian based distributed compressed sensing (HVB-DCS) algorithm for the JSM-1 model, which decouples the common component from the innovation components by applying variational approximation. The proposed algorithm is capable of approximating the complex signal in a hierarchical framework. Moreover, we apply the proposed HVB-DCS algorithm for signal reconstruction in the dual-channel SAR based GMTI system. Based on the received signal from both the moving targets and the stationary clutter, we decouple the common component (stationary clutter) from the innovation components (moving targets). We derive separable PDF functions for common component and innovation components, respectively. We can then obtain the reconstructed signal (including the common and innovation components) by applying the maximum a posteriori (MAP) estimation. The performance of the algorithm is verified by both simulation on point targets, and experiment on real SAR data.

The main contribution of this paper is threefold. First, we model the received signals of the dual-channel SAR-GMTI system in the JSM-1 model of DCS, by exploiting the correlation between the two channels. Secondly, we propose a novel HVB-DCS algorithm for the JSM-1 model, which decouples the common component from the innovation components by applying variational approximation. The proposed HVB-DCS algorithm approximates the complex signal in a hierarchical framework. Thirdly, we apply the proposed HVB-DCS algorithm for signal reconstruction in the dual-channel SAR-GMTI system. In the framework of the JSM-1 model, the proposed algorithm can simultaneously generate the stationary clutter (common component) and the moving targets (innovation components), which omits the DPCA process adopted in the traditional GMTI system. The simulation and experimental results show that the proposed HVB-DCS algorithm is capable of detecting multiple moving targets while suppressing the clutter at a much lower data rate compared with the CS and range-Doppler (RD) algorithms.

The organization of this paper is as follows. First, we briefly introduce the theory of DCS in Section 2. We then present the signal model of dual-channel SAR system in Section 3. The DCS-based dual-channel SAR-GMTI system is introduced in Section 4, where the echo signals of dual-channel SAR-GMTI system are represented in the JSM-1 model. The HVB-DCS algorithm is proposed in Section 5, which decouples the common component from the innovation components by applying variational approximation, in a hierarchical framework. The simulation and experimental results are shown in Section 6, and the conclusions are summarized in Section 7.

2. Brief Introduction to Distributed Compressed Sensing

CS [4] is an emerging field for signal and image processing. CS aims to recover a sparse signal with a reduced number of linear measurements, which is significantly less than that required by the traditional bandwidth constraint based on the Shannon-Nyquist sampling theorem. CS theory has been widely applied in various fields, such as magnetic resonance imaging [27], target tracking [28], biomedical monitoring [29], computer-generated holography [30], machine learning [31], signal acquisition in wireless sensor networks [32], image/video compression [33], etc.

As an important branch of CS, DCS [18] considers a scenario where correlated sparse signals are captured from multiple sources. DCS aims to exploit the inter-signal and intra-signal correlations of those sparse signals, so that the ensemble of signals can be reconstructed from fewer measurements than the standard CS approach requires. Baron et al. in [18] propose three various joint sparsity models (JSMs) for the correlated signals, namely JSM-1, JSM-2 and JSM-3. Among them, JSM-1 is capable of modeling a wide range of scenarios where large-scale phenomena affect all signal sources and local phenomena affect specific individual source. In this work, we focus on the JSM-1 model.

2.1. Joint Sparsity Model-1

In the JSM-1 framework, an ensemble of vectors (signals) ${\left\{{\mathbf{x}}_j\right\}}_{j=1}^J$ are jointly sparse if each vector is sparse and comprised of a common component and an innovation component, as

$${\mathbf{x}}_j={\mathbf{z}}_c+{\mathbf{z}}_j, j\in \{1,2,\cdots ,J\},$$

(1)

where ${\mathbf{z}}_c\in {\Re }^N$ denotes the common component of sparse vector ${\mathbf{x}}_j$, and ${\mathbf{z}}_j\in {\Re }^N$ denotes the innovation component of ${\mathbf{x}}_j$. Both ${\mathbf{z}}_c$ and ${\mathbf{z}}_j$ are supposed to be sparse.

The DCS problem consists in reconstructing jointly sparse vectors from their measurements, which are obtained by

$${\mathbf{y}}_j={\mathbf{A}}_j{\mathbf{x}}_j,$$

(2)

where ${\mathbf{A}}_j\in {\Re }^{M_j\times N}$ is the measurement matrix, ${\mathbf{y}}_j\in {\Re }^{M_j}$ is the measurement vector, and $M_j$ is the number of measurements.

2.2. The Joint Recovery Strategy for the JSM-1 Model

JRS is the first proposed algorithm for the JSM-1 model. The main idea is to stack the common component ${\mathbf{z}}_c$ and all innovation components ${\left\{{\mathbf{z}}_j\right\}}_{j=1}^J$ together to form one sparse vector $\tilde{\mathbf{z}}$ with multiple dimension $(J+1)N$, and the similar process for measurement vectors ${\left\{{\mathbf{y}}_j\right\}}_{j=1}^J$ and measurement matrices ${\left\{{\mathbf{A}}_j\right\}}_{j=1}^J$, as

$$\begin{array}{c}\tilde{\mathbf{z}}:=\left[\begin{array}{c}{\mathbf{z}}_c\\ {\mathbf{z}}_1\\ {\mathbf{z}}_2\\ ⋮\\ {\mathbf{z}}_J\end{array}\right],\tilde{\mathbf{y}}:=\left[\begin{array}{c}{\mathbf{y}}_1\\ {\mathbf{y}}_2\\ ⋮\\ {\mathbf{y}}_J\end{array}\right],\tilde{\mathbf{A}}:=\left[\begin{array}{ccccc}{\mathbf{A}}_1& {\mathbf{A}}_1& \mathbf{0}& \cdots & \mathbf{0}\\ {\mathbf{A}}_2& \mathbf{0}& {\mathbf{A}}_2& \cdots & \mathbf{0}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\mathbf{A}}_J& \mathbf{0}& \mathbf{0}& \cdots & {\mathbf{A}}_J\end{array}\right].\end{array}$$

(3)

The J measurement Equations (2) can be formulated in one compact form as

$$\tilde{\mathbf{y}}=\tilde{\mathbf{A}}\tilde{\mathbf{z}},$$

(4)

resulting in the formulation of a standard CS problem. Then the sparse vector $\tilde{\mathbf{z}}$ is reconstructed based on Equation (4) via well developed algorithms of CS, such as orthogonal matching pursuit (OMP) and BP. Finally, the jointly sparse vectors ${\left\{{\mathbf{x}}_j\right\}}_{j=1}^J$ are reconstructed by adding the reconstructed common component to each reconstructed innovation component.

However, the computational costs of JRS for the JSM-1 model increase dramatically when the dimension of sparse vector increases. So it is necessary to develop another more efficient recovery strategy for the JSM-1 model.

3. Signal Model of the Dual-Channel SAR System

This section introduces the signal model of a dual-channel SAR system, which is the same as that in [34]. The readers can refer to [34] for more details about the signal model. The geometry relationship between the flying platform and moving target is shown in Figure 1, where v and H denote velocity and height of the platform respectively. Channel one periodically transmits and receives pulses, while channel two receives pulses. The distance between the two channels is d. $t_m$ is the slow time in azimuth. PRF represents the pulse repetition frequency. The target is denoted by P. $R_B$ denotes the nearest slant range between the platform and target, whereas $R_1\left(t_m\right)$ is the instantaneous slant range between channel one and the moving target, and the similar definition for $R_2\left(t_m\right)$.

In order to facilitate the analysis, the slant range history geometry of the moving target is shown in Figure 2. $v_a$ and $v_r$ denote the along-track and cross-track velocities of the moving target on the slant range plane, respectively. The instantaneous slant ranges $R_1\left(t_m\right)$ and $R_2\left(t_m\right)$ are respectively defined as

$$\begin{array}{c}{\displaystyle R_1(t_m;R_B)=\sqrt{{(v{t}_m-v_a{t}_m)}^2+{(R_B-v_r{t}_m)}^2}\approx R_B-v_r{t}_m+\frac{{(v-v_a)}^2}{2R_B}{t}_m^2,}\hfill \end{array}$$

(5)

and

$$\begin{array}{c}{\displaystyle R_2(t_m;R_B)=\sqrt{{(v{t}_m-v_a{t}_m-d)}^2+{(R_B-v_r{t}_m)}^2}\approx R_B-v_r{t}_m+\frac{{\left[(v-v_a)t_m-d\right]}^2}{2R_B},}\hfill \end{array}$$

(6)

where

$$\begin{array}{c}{\displaystyle t_m\in \{t_1,t_2\cdots ,t_i,\cdots ,t_M\}-\frac{t_1+t_M}{2},}\hfill \end{array}$$

and M denotes the number of pulses transmitted by antennas.

For channel one, the received signal of moving target after demodulation and range compression can be expressed as

$$\begin{array}{c}{\displaystyle s_1(t,t_m)=\sigma G\mathrm{sinc}\left(\Delta B\left(t-\frac{2R_1\left(t_m;R_B\right)}{c}\right)\right){\omega }_a\left(t_m\right)exp\left(-j\frac{4\pi }{\lambda }{R}_1\left(t_m;R_B\right)\right).}\hfill \end{array}$$

(7)

The received signal for channel two is similarly defined as

$$\begin{array}{c}{s}_2(t,t_m)=\sigma G\mathrm{sinc}\left(\Delta B\left(t-\frac{R_1\left(t_m;R_B\right)+R_2\left(t_m;R_B\right)}{c}\right)\right){\omega }_a\left(t_m\right)exp(-j\frac{2\pi }{\lambda }\left(R_1\left(t_m;R_B\right)+R_2\left(t_m;R_B\right)\right)).\hfill \end{array}$$

(8)

In the above expressions, t is the fast time in range, $\sigma$ denotes the complex reflectivity of the target, and G is the range compression gain. $\Delta B$ is the bandwidth of transmitted signal. $\mathrm{sinc}(·)$ is the envelope after range compression. ${\omega }_a(·)$ is the azimuth windowing function. c and $\lambda$ represent the speed of light and the wavelength of transmitted signal, respectively.

4. Distributed Compressed Sensing Based Dual-Channel SAR-GMTI System

4.1. Sparse Representation for Individual Channel

In this section, the random sampling mode in the azimuth direction [34] is adopted for the dual-channel SAR system, which is shown in Figure 3. In this mode, the pulses in azimuth are randomly transmitted and received, which results in time gaps within a coherent processing interval (CPI) where no echoes are recorded. Moreover, we assume that the sensors remain still between transmission and reception of a pulse.

For channel one, we can obtain Equation (9) by substituting Equation (5) into Equation (7).

$$\begin{array}{c}{\displaystyle s_1(t,t_m)=\sigma G\mathrm{sinc}\left(\Delta \mathrm{B}\left(t-\frac{2R_B}{c}\right)\right){\omega }_a\left(t_m\right)exp\left(-j\frac{4\pi }{\lambda }{R}_B\right)exp\left(j2\pi f_{dc}{t}_m+j\pi {\gamma }_m\left(R_B\right)t_m^2\right),}\hfill \end{array}$$

(9)

where $f_{dc}=2v_r/\lambda$ is the Doppler centroid frequency and ${\gamma }_m\left(R_B\right)=-2{(v-v_a)}^2/\left(\lambda R_B\right)$ is the corresponding Doppler chirp rate. Thus Equation (9) can be expanded as

$$\begin{array}{c}{\displaystyle s_1(t,t_m)=\sigma G\mathrm{sinc}\left(\Delta \mathrm{B}\left(t-\frac{2R_B}{c}\right)\right){\omega }_a\left(t_m\right)exp\left(-j\frac{4\pi }{\lambda }{R}_B\right)exp\left(j2\pi \frac{v_r^2{R}_B}{\lambda {(v-v_a)}^2}\right)}\hfill \\ {\displaystyle \times exp\left(-j2\pi \frac{{\left(v-v_a\right)}^2}{\lambda R_B}{\left(t_m-\frac{v_r{R}_B}{{\left(v-v_a\right)}^2}\right)}^2\right).}\hfill \end{array}$$

(10)

Assuming that the cross-track velocity of the moving point target is relatively slow, the range migration through resolution cells does not take place.

Remark: If the cross-track velocity of the moving target is high enough, the Doppler centroid frequency $f_{dc}$ may exceed the limit of PRF, which induces the Doppler ambiguity problem and an additional range walk [35]. For the condition of high cross-track velocity, an additional range walk correction step is required after range compression. As the cross-track velocity is unknown, it is required to search the Doppler ambiguity number and estimate the cross-track velocity. In recent years, many algorithms have been proposed to estimate the Doppler ambiguity number and the cross-track velocity, e.g., the road slope-aided algorithm [36] and the range walk correcting-based algorithm [37]. In the future work, we will combine the proposed HVB-DCS algorithm and the above-mentioned estimation algorithms to tackle the high cross-track velocity problem.

The received signal at a given range bin data, after demodulation and range compression, can be rewritten as

$$\begin{array}{c}\hfill \begin{array}{c}\hfill s_1\left(t_m\right)={\rho }^{\left(1\right)}{\omega }_a\left(t_m\right)exp\left(-j2\pi \frac{{\left(v-v_a\right)}^2}{\lambda R_B}{\left(t_m-\frac{v_r{R}_B}{{\left(v-v_a\right)}^2}\right)}^2\right),\end{array}\end{array}$$

(11)

where

$$\begin{array}{c}\hfill \begin{array}{c}{\rho }^{\left(1\right)}=\sigma G\mathrm{sinc}\left(\Delta \mathrm{B}\left(t-\frac{2R_B}{c}\right)\right)exp\left(-j\frac{4\pi }{\lambda }{R}_B\right)exp\left(j2\pi \frac{v_r^2{R}_B}{\lambda {\left(v-v_a\right)}^2}\right).\hfill \end{array}\end{array}$$

(12)

We assume that there are no more than $N_0$ scattering centers which can be distinguished in the synthetic aperture time $T_0$. The coordinates of these scatterers in the azimuth direction are $x_n(n=1,2,\cdots ,N_0)$. The received signal at a given range bin data, after demodulation and range compression, can then be represented as a combination of the echo signals from $N_0$ scattering centres as,

$$\begin{array}{c}\hfill \begin{array}{c}{s}_1\left(t_m\right)=\sum _{n=1}^{N_0}{\rho }_n^{\left(1\right)}{\omega }_a(t_m-x_n/v)exp\left(-j2\pi \frac{{\left(v-v_a\right)}^2}{\lambda R_B}{\left(t_m-x_n/v-\frac{v_r{R}_B}{{\left(v-v_a\right)}^2}\right)}^2\right),\hfill \end{array}\end{array}$$

(13)

where

$${\rho }_n^{\left(1\right)}={\sigma }_{n}G\mathrm{sinc}\left(\Delta \mathrm{B}\left(t-\frac{2R_B}{c}\right)\right)exp\left(-j\frac{4\pi }{\lambda }{R}_B\right).$$

(14)

The along-track velocity of the target defocuses the image, thus in order to facilitate the analysis, we only consider the cross-track velocity and let $v_a=0$.

Next, we represent the signal received by channel one in a standard CS framework. We first build the measurement matrix according to Equation (13), as

$$\begin{array}{c}{\mathsf{\Phi }}_1={\left[{\mathbf{s}}_0(-\frac{N}{2}+1),{\mathbf{s}}_0(-\frac{N}{2}),\cdots ,{\mathbf{s}}_0\left(0\right),\cdots ,{\mathbf{s}}_0(\frac{N}{2}-1),{\mathbf{s}}_0\left(\frac{N}{2}\right)\right]}_{M\times N},\hfill \end{array}$$

(15)

where,

$${\mathbf{s}}_0\left(i\right)={\left[s_0(t_1-i\Delta \tau ),s_0(t_2-i\Delta \tau ),\cdots ,s_0(t_M-i\Delta \tau )\right]}_M^T,$$

(16)

$$s_0\left(t_m-i\Delta \tau \right)=\left\{\begin{array}{ccc}\hfill exp\left(j\pi {\gamma }_m\left(R_B\right){\left(t_m-i\Delta \tau \right)}^2\right)& ,\hfill & \hfill \left|t_m-i\Delta \tau \right|\le T/2,\\ \hfill 0& ,\hfill & \hfill \left|t_m-i\Delta \tau \right|>T/2,\end{array}\right.$$

(17)

$$\begin{array}{c}{\displaystyle i\in \left\{-\frac{N}{2}+1,\cdots ,0,\cdots ,\frac{N}{2}\right\},}\hfill \end{array}$$

${\gamma }_m\left(R_B\right)=-2{(v-v_a)}^2/\left(\lambda R_B\right)$, T is the full synthetic aperture time, $\Delta \tau =1/\mathrm{PRF}$, and $N\ge (T_0+T)/\Delta \tau$. N is the number of samples in the azimuth direction. M is far less than N. The measurement matrix captures the contribution to the received signal of a point target.

We then define the measurement vector ${\mathbf{s}}_1\in {\mathcal{C}}^M$ as the received signal for a given range bin after range compression by channel one, ${\mathbf{s}}_1={\left[s_1\left(t_1\right),s_1\left(t_2\right),\cdots ,s_1\left(t_M\right)\right]}^T$, and the sparse vector ${\mathbf{\rho }}_1\in {\mathcal{C}}^N$ as the complex image for a given range bin, ${\mathbf{\rho }}_1={\left[{\rho }_1^{\left(1\right)},{\rho }_2^{\left(1\right)},\cdots ,{\rho }_N^{\left(1\right)}\right]}^T$. Thus we can obtain the standard equation in CS for channel one as

$$\begin{array}{c}{\displaystyle {\mathbf{s}}_1={\mathsf{\Phi }}_1{\mathbf{\rho }}_1.}\hfill \end{array}$$

(18)

Similarly, for channel two, Equation (6) is substituted into Equation (8), resulting in

$$\begin{array}{cc}{s}_2\left(t,t_m\right)& {\displaystyle =\sigma G\mathrm{sinc}\left(\Delta \mathrm{B}\left(t-\frac{2R_{\mathrm{B}}}{c}\right)\right){\omega }_a\left(t_m\right)exp\left(-j\frac{4\pi }{\lambda }{R}_B\right)exp\left(j2\pi \frac{v_r^2{R}_B}{\lambda {\left(v-v_a\right)}^2}\right)exp\left(-j\pi \frac{d^2}{2\lambda R_B}\right)}\hfill \\ & {\displaystyle \times exp\left(-j2\pi \frac{{\left(v-v_a\right)}^2}{\lambda R_B}{\left(t_m-\frac{v_r{R}_B}{{\left(v-v_a\right)}^2}-\frac{d}{2\left(v-v_a\right)}\right)}^2\right)exp\left(j\frac{2\pi }{\lambda }\frac{v_{r}d}{\left(v-v_a\right)}\right).}\hfill \end{array}$$

(19)

Comparing Equation (19) (for channel two) with Equation (10) (for channel one), we can find that there is a phase difference $exp\left(-j\pi d^2/\left(2\lambda R_B\right)\right)$ and a time delay $d/\left(2\left(v-v_a\right)\right)$ between the two channels. Thus the measurement matrix for channel two can be represented as

$$\begin{array}{c}{\mathsf{\Phi }}_2={\left[{{\mathbf{s}}_0}^{\prime }\left(-\frac{N}{2}+1\right),{{\mathbf{s}}_0}^{\prime }\left(-\frac{N}{2}\right),\cdots ,{{\mathbf{s}}_0}^{\prime }\left(0\right),\cdots ,{{\mathbf{s}}_0}^{\prime }\left(\frac{N}{2}-1\right),{{\mathbf{s}}_0}^{\prime }\left(\frac{N}{2}\right)\right]}_{M\times N},\hfill \end{array}$$

(20)

where,

$${\mathbf{s}}_0^{\prime }\left(i\right)={\left[s_0^{\prime }\left(t_1-i\Delta \tau \right),s_0^{\prime }\left(t_2-i\Delta \tau \right),\cdots ,s_0^{\prime }\left(t_M-i\Delta \tau \right)\right]}_M^T,$$

(21)

$${s_0}^{\prime }(t_m-i\Delta \tau )=\left\{\begin{array}{ccc}\hfill exp\left(j\pi {\gamma }_m\left(R_B\right){\left(t_m-i\Delta \tau -\frac{d}{2\left(v-v_a\right)}\right)}^2\right)& ,\hfill & \hfill \left|t_m-i\Delta \tau -\frac{d}{2(v-v_a)}\right|\le T/2,\\ \hfill 0& ,\hfill & \hfill \left|t_m-i\Delta \tau -\frac{d}{2(v-v_a)}\right|>T/2,\end{array}\right.$$

(22)

$$\begin{array}{c}{\displaystyle i\in \left\{-\frac{N}{2}+1,\cdots ,0,\cdots ,\frac{N}{2}\right\}.}\hfill \end{array}$$

We define the measurement vector ${\mathbf{s}}_2\in {\mathcal{C}}^M$ as the received signal for a given range bin after range compression by channel two, ${\mathbf{s}}_2={\left[s_2\left(t_1\right),s_2\left(t_2\right),\cdots ,s_2\left(t_M\right)\right]}^T$ and the sparse vector ${\mathbf{\rho }}_2\in {\mathcal{C}}^N$ as ${\mathbf{\rho }}_2={\left[{\rho }_1^{\left(2\right)},{\rho }_2^{\left(2\right)},\cdots ,{\rho }_N^{\left(2\right)}\right]}^T$. Thus, we obtain the standard equation in CS for channel two as

$$\begin{array}{c}{\displaystyle {\mathbf{s}}_2={\mathsf{\Phi }}_2{\mathbf{\rho }}_2.}\hfill \end{array}$$

(23)

As we can see, the measurement matrices are built using the echo signal models of moving targets as in Equations (13) and (19). In practice, it is difficult to construct the measurement matrices since the along-track velocity $v_a$ and cross-track velocity $v_r$ of the moving target are unknown. Thus, in order to facilitate the analysis and construction of the measurement matrices, we set the values of $v_a$ and $v_r$ to zero. The effects of the nonzero cross-track velocity and nonzero along-track velocity on SAR imaging are analyzed in detail through simulations performed in Section 6.1.

4.2. JSM-1 Model

For the $i^{th}$ scattering center, its corresponding reflection coefficients in channel one and channel two can be respectively represented as

$$\begin{array}{c}{\displaystyle {\rho }_i^{\left(1\right)}={\sigma }_{i}G\mathrm{sinc}\left(\Delta B\left(t-\frac{2R_B}{c}\right)\right)exp\left(-j\frac{4\pi }{\lambda }{R}_B\right)exp\left(j\frac{2\pi v_r^2{R}_B}{\lambda v^2}\right),}\hfill \end{array}$$

(24)

and

$$\begin{array}{c}{\displaystyle {\rho }_i^{\left(2\right)}={\sigma }_{i}G\mathrm{sinc}\left(\Delta B\left(t-\frac{2R_B}{c}\right)\right)exp\left(-j\frac{4\pi }{\lambda }{R}_B\right)exp\left(j\frac{2\pi v_r^2{R}_B}{\lambda v^2}\right)}\hfill \\ {\displaystyle \times exp\left(-j\pi \frac{d^2}{2R_B\lambda }\right)exp\left(j\frac{2\pi }{\lambda }\left(v_r\frac{d}{v}\right)\right).}\hfill \end{array}$$

(25)

Equations (24) and (25) demonstrate the close connection between the two channels, i.e., except the items of a fixed phase difference $exp\left(-j\pi d^2/\left(2\lambda R_B\right)\right)$, and a varying phase difference $exp\left(j\frac{2\pi }{\lambda }\left(v_r\frac{d}{v}\right)\right)$ changing with cross-track velocity in Equation (25), the other items on the reflection coefficients of the two channels are the same.

First, we consider compensating for the fixed phase difference on the $i^{th}$ reflection coefficient of channel two, ${\rho }_i^{\left(2\right)}$.

$$\begin{array}{cc}\hfill {\rho }_i^{{\left(2\right)}^{\prime }}& {\displaystyle ={\rho }_i^{\left(2\right)}exp\left(j\pi d^2/2\lambda R_B\right)={\sigma }_{i}G\mathrm{sinc}\left(\Delta B\left(t-\frac{2R_B}{c}\right)\right)}\hfill \\ & {\displaystyle ·exp\left(-j\frac{4\pi }{\lambda }{R}_B\right)exp\left(j\frac{2\pi v_r^2{R}_B}{\lambda v^2}\right)exp\left(j\frac{2\pi }{\lambda }\left(v_r\frac{d}{v}\right)\right),i=1,\cdots ,N.}\hfill \end{array}$$

(26)

The compensating procedure in Equation (26) is repeated for N reflection coefficients of channel two, resulting in a new sparse vector, i.e., ${\mathbf{\rho }}_2^{\prime }={\left[{\rho }_1^{{\left(2\right)}^{\prime }},{\rho }_2^{{\left(2\right)}^{\prime }},\cdots ,{\rho }_N^{{\left(2\right)}^{\prime }}\right]}^T$, namely the compensated reflection coefficient vector. Furthermore, we can obtain the relationships between the reflection coefficient vectors of channel two, respectively after and before the compensating procedure, together with the reflection coefficient vector of channel one, as

$$\begin{array}{c}\left[\begin{array}{c}{\mathbf{\rho }}_1\\ {\displaystyle {\mathbf{\rho }}_2^{\prime }}\end{array}\right]=\left[\begin{array}{cc}\mathbf{I}& \mathbf{0}\\ \mathbf{0}& {\displaystyle exp\left(j\pi \frac{d^2}{2\lambda R_B}\right)\mathbf{I}}\end{array}\right]\left[\begin{array}{c}{\mathbf{\rho }}_1\\ {\displaystyle {\mathbf{\rho }}_2}\end{array}\right],\hfill \end{array}$$

(27)

where $\mathbf{I}$ denotes an identity matrix of appropriate size.

Secondly, we build the JSM-1 model based on the original reflection coefficient vector of channel one, ${\mathbf{\rho }}_1$, and the compensated reflection coefficient vector of channel two, ${\mathbf{\rho }}_2^{\prime }$. The nonzero elements of ${\mathbf{\rho }}_1$ consist of the reflection coefficients from both the stationary clutter and moving targets, which are denoted as ${\rho }_s^{\left(1\right)}$ and ${\rho }_m^{\left(1\right)}$, respectively. Similarly, the nonzero elements of ${\mathbf{\rho }}_2^{\prime }$ consist of the reflection coefficients from both the stationary clutter and moving targets, which are denoted as ${\rho }_s^{{\left(2\right)}^{\prime }}$ and ${\rho }_m^{{\left(2\right)}^{\prime }}$, respectively. On the one hand, the reflection coefficients of stationary clutter are common to both ${\mathbf{\rho }}_1$ and ${\mathbf{\rho }}_2^{\prime }$, since $v_r=0$ and we have ${\rho }_s^{\left(1\right)}={\rho }_s^{{\left(2\right)}^{\prime }}$. On the other hand, the reflection coefficients of moving targets are different for channel one and channel two, since $v_r\ne 0$, thus ${\rho }_m^{\left(1\right)}\ne {\rho }_m^{{\left(2\right)}^{\prime }}$. Let ${\mathbf{z}}_c$ denote the vector consisting of the reflection coefficients of stationary clutter, ${\mathbf{z}}_1$ and ${\mathbf{z}}_2$ denote the vectors consisting of the reflection coefficients of moving targets observed in channel one and channel two, respectively. Then we have

$$\begin{array}{c}\left[\begin{array}{c}{\mathbf{\rho }}_1\\ {\displaystyle {\mathbf{\rho }}_2^{\prime }}\end{array}\right]=\left[\begin{array}{ccc}\mathbf{I}& \mathbf{I}& \mathbf{0}\\ \mathbf{I}& \mathbf{0}& \mathbf{I}\end{array}\right]\left[\begin{array}{c}{\mathbf{z}}_c\hfill \\ {\mathbf{z}}_1\hfill \\ {\mathbf{z}}_2\hfill \end{array}\right].\hfill \end{array}$$

(28)

We define $\mathbf{y}={\left[{\mathbf{s}}_1^T{\mathbf{s}}_2^T\right]}^T$, and can obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\mathbf{y}\hfill & =\left[\begin{array}{cc}{\mathsf{\Phi }}_1& \mathbf{0}\\ \mathbf{0}& {\mathsf{\Phi }}_2\end{array}\right]\left[\begin{array}{c}{\mathbf{\rho }}_1\hfill \\ {\mathbf{\rho }}_2\hfill \end{array}\right]=\left[\begin{array}{cc}{\mathsf{\Phi }}_1& \mathbf{0}\\ \mathbf{0}& {\mathsf{\Phi }}_2\end{array}\right]\left[\begin{array}{cc}\mathbf{I}& \mathbf{0}\\ \mathbf{0}& {\displaystyle exp\left(-j\pi \frac{d^2}{2\lambda R_B}\right)\mathbf{I}}\end{array}\right]\left[\begin{array}{c}{\mathbf{\rho }}_1\\ {\displaystyle {\mathbf{\rho }}_2^{\prime }}\end{array}\right]\hfill \\ & =\left[\begin{array}{cc}{\mathsf{\Phi }}_1& \mathbf{0}\\ \mathbf{0}& {\mathsf{\Phi }}_2\end{array}\right]\left[\begin{array}{cc}\mathbf{I}& \mathbf{0}\\ \mathbf{0}& {\displaystyle exp\left(-j\pi \frac{d^2}{2\lambda R_B}\right)\mathbf{I}}\end{array}\right]\left[\begin{array}{ccc}\mathbf{I}& \mathbf{I}& \mathbf{0}\\ \mathbf{I}& \mathbf{0}& \mathbf{I}\end{array}\right]\left[\begin{array}{c}{\mathbf{z}}_c\hfill \\ {\mathbf{z}}_1\hfill \\ {\mathbf{z}}_2\hfill \end{array}\right]\hfill \\ & =\left[\begin{array}{ccc}{\mathsf{\Phi }}_1& {\mathsf{\Phi }}_1& \mathbf{0}\\ {\mathsf{\Phi }}_{2}exp\left(-j\pi \frac{d^2}{2\lambda R_B}\right)& \mathbf{0}& {\mathsf{\Phi }}_{2}exp\left(-j\pi \frac{d^2}{2\lambda R_B}\right)\end{array}\right]\left[\begin{array}{c}{\mathbf{z}}_c\hfill \\ {\mathbf{z}}_1\hfill \\ {\mathbf{z}}_2\hfill \end{array}\right]\hfill \end{array}.\end{array}$$

(29)

In Equation (29), the first equality is based on Equations (18) and (23), and the second equality holds considering the inverse of Equation (27).

Equation (29) is a standard JSM-1 model. By jointly reconstructing the common and innovation components of the above, the moving targets ${\mathbf{z}}_1$ and ${\mathbf{z}}_2$ can be separated. However, the computational costs of some existing algorithms for solving the above model increase dramatically when the dimension of sparse vector increases. So we develop a more efficient recovery strategy for this JSM-1 model.

5. The Hierarchical Variational Bayesian Based DCS Algorithm

5.1. Proposed Algorithm

In Bayesian modeling, all unknowns are treated as stochastic quantities with assigned probability distributions. Considering the effects of noise, the standard JSM-1 model can be presented as

$${\mathbf{y}}_k={\mathbf{A}}_k({\mathbf{z}}_c+{\mathbf{z}}_k)+{\mathbf{v}}_k,$$

(30)

where ${\mathbf{y}}_k\in {\Re }^M$, ${\mathbf{A}}_k\in {\Re }^{M\times N}$ and ${\mathbf{v}}_k\in {\Re }^M$ denote the measurement vector, the measurement matrix and noise vector, of channel k, respectively; ${\mathbf{z}}_c,{\mathbf{z}}_k\in {\Re }^N$ denote the common component and the innovation component respectively.

By applying a suitable prior distribution to ${\mathbf{z}}_c$ and ${\mathbf{z}}_k$, the sparsity can be guaranteed. However, it is important to allow the flexibility to model local characteristics of the signal, as the simple stationary sparse prior distribution is unable to meet the demand. For this reason, we propose an HVB-DCS algorithm for the JSM-1 model.

We adopt zero-mean Gaussian prior distributions for the common component and innovation components, respectively, which are given as

$$\begin{array}{c}{\displaystyle p\left({\mathbf{z}}_c|{\mathbf{\alpha }}_c\right)=\prod _{n=1}^N\mathcal{N}\left(z_c^n|0,{\left({\alpha }_c^n\right)}^{-1}\right),}\hfill \end{array}$$

(31)

and

$$\begin{array}{c}{\displaystyle p\left({\mathbf{z}}_k|{\mathbf{\alpha }}_k\right)=\prod _{n=1}^N\mathcal{N}\left(z_k^n|0,{\left({\alpha }_k^n\right)}^{-1}\right),}\hfill \end{array}$$

(32)

where $z_c^n$, $z_k^n$, ${\alpha }_c^n$, and ${\alpha }_k^n$ are the $n^{th}$ element of ${\mathbf{z}}_c$, ${\mathbf{z}}_k$, ${\mathbf{\alpha }}_c$, and ${\mathbf{\alpha }}_k$, respectively. The precision parameters ${\mathbf{\alpha }}_c$ and ${\mathbf{\alpha }}_k$ are constrained by treating them as random variables, with Gamma prior distributions as

$$\begin{array}{c}{\displaystyle p\left({\mathbf{\alpha }}_c;a,b\right)=\prod _{n=1}^{N}Gamma\left({\alpha }_c^n|a,b\right),}\hfill \end{array}$$

(33)

and

$$\begin{array}{c}{\displaystyle p\left({\mathbf{\alpha }}_k;c_k,d_k\right)=\prod _{n=1}^{N}Gamma\left({\alpha }_k^n|c_k,d_k\right).}\hfill \end{array}$$

(34)

Furthermore, we assume that the $m^{th}$ element of noise vector, $v_k^m(m=1,\cdots ,M),$ obeys an independent and identically distributed (i.i.d) zero-mean Gaussian distribution with inverse variance $\beta$, i.e., $v_k^m\sim \mathcal{N}\left(v_k^m\right|0,{\beta }^{-1})$. Thus we can obtain the posterior distribution of noise vector ${\mathbf{v}}_k$ as the product of that of the individual element as,

$$p\left({\mathbf{v}}_k\right|\beta )=\prod _{m=1}^M\mathcal{N}\left(v_k^m\right|0,{\beta }^{-1}).$$

(35)

We further assume a Gamma distribution as prior for the noise inverse variance $\beta$

$$\begin{array}{c}{\displaystyle p(\beta ;e,f)=Gamma\left(\beta \right|e,f),}\hfill \end{array}$$

(36)

where e and f represent the shape and scale parameters of the Gamma distribution, respectively.

We define $\mathbf{Y}=\{{\mathbf{y}}_1,{\mathbf{y}}_2,\dots ,{\mathbf{y}}_k\}$, $\mathbf{Z}=\{{\mathbf{z}}_c,{\mathbf{z}}_1,\dots ,{\mathbf{z}}_k,{\mathbf{\alpha }}_c,\dots ,{\mathbf{\alpha }}_k,\beta \}$, and $\mathbf{\theta }=\{a,b,c_1\dots c_k,d_1\dots d_k,e,f\}$ as the sets of measurement vectors, the hidden variables, and the hyperparameters of the imposed prior, respectively. Figure 4 describes the relationships among the measurements (indicated as a doubly circled node), the hidden variables (indicated as single circled nodes) and the hyperparameters (indicated as square nodes). The directed edges of the graphical model represent the dependencies among the variables. For instance, the measurements $\mathbf{Y}$ relies on the hidden variables ${\mathbf{z}}_c$, ${\mathbf{z}}_1$, ${\mathbf{z}}_2$, ⋯, ${\mathbf{z}}_k$ and $\beta$, while the hidden variable ${\mathbf{z}}_c$ stochastically depends on the hidden variable ${\mathbf{\alpha }}_c$, and ${\mathbf{\alpha }}_c$ further relies on the model parameters $a,b$.

Next, we are to pursuit the posterior distribution of the hidden variable, $p\left(\mathbf{Z}\right|\mathbf{Y};\mathbf{\theta })$. In [24], VB introduces a variational distribution $q\left(\mathbf{Z}\right)$ to approximate the true posterior distribution $p\left(\mathbf{Z}\right|\mathbf{Y};\mathbf{\theta })$. First, the log marginal likelihood $logp(\mathbf{Y};\mathbf{\theta })$ can be represented as

$$\begin{array}{c}logp\left(\mathbf{Y};\mathbf{\theta }\right)=F\left(q\left(\mathbf{Z}\right),\mathbf{\theta }\right)+KL\left(q\left(\mathbf{Z}\right)\left|\right|p\left(\mathbf{Z}|\mathbf{Y};\mathbf{\theta }\right)\right),\hfill \end{array}$$

(37)

where $F\left(q\left(\mathbf{Z}\right),\mathbf{\theta }\right)$ is the free energy,

$$\begin{array}{c}F\left(q\left(\mathbf{Z}\right),\mathbf{\theta }\right)=\int {\displaystyle q\left(\mathbf{Z}\right)log\left(\frac{p\left(\mathbf{Z},\mathbf{Y};\mathbf{\theta }\right)}{q\left(\mathbf{Z}\right)}\right)d\mathbf{Z}},\hfill \end{array}$$

(38)

and $\mathrm{KL}\left(\mathrm{q}\left(\mathbf{Z}\right)\left|\right|\mathrm{p}\left(\mathbf{Z}|\mathbf{Y};\mathbf{\theta }\right)\right)$ is the Kullback-Leibler (KL) divergence between the true posterior $p(\mathbf{Z}|\mathbf{Y};\mathbf{\theta })$ and the variational distribution $q\left(\mathbf{Z}\right)$,

$$\begin{array}{c}\mathrm{KL}\left(\mathrm{q}\left(\mathbf{Z}\right)\left|\right|\mathrm{p}(\mathbf{Z}|\mathbf{Y};\mathbf{\theta })\right)=\int {\displaystyle \mathrm{q}\left(\mathbf{Z}\right)\log \left(\frac{\mathrm{q}\left(\mathbf{Z}\right)}{\mathrm{p}\left(\mathbf{Z}\right|\mathbf{Y};\mathbf{\theta })}\right)\mathrm{d}\mathbf{Z}}.\hfill \end{array}$$

(39)

The goal is to approximate the true posterior distribution by minimizing the KL. Due to the fact that $\mathrm{KL}\left(\mathrm{q}\left(\mathbf{Z}\right)\left|\right|\mathrm{p}(\mathbf{Z}|\mathbf{Y};\mathbf{\theta })\right)\ge 0$, that objective can be achieved by maximizing $F\left(q\left(\mathbf{Z}\right),\mathbf{\theta }\right)$. For the JSM-1 model, we should not only find separable functions that approximate the posterior distribution of ${\mathbf{z}}_c$ and ${\mathbf{z}}_k$, but also make the integral $F\left(q\left(\mathbf{Z}\right),\mathbf{\theta }\right)$ tractable.

In order to meet the requirements, we can assume that $q\left(\mathbf{Z}\right)$ has a factorized form Equation (40). This factorized form stems from theoretical physics where it is called mean field theory [38].

$${\displaystyle \begin{array}{c}q\left(\mathbf{Z}\right)=q\left({\mathbf{z}}_c\right)q\left({\mathbf{z}}_1\right)\cdots q\left({\mathbf{z}}_k\right)q\left({\mathbf{\alpha }}_c\right)q\left({\mathbf{\alpha }}_1\right)\cdots q\left({\mathbf{\alpha }}_k\right)q\left(\beta \right).\hfill \end{array}}$$

(40)

Optimizing the free energy in Equation (38) is realized by taking functional derivatives with respect to each of $q(·)$ distributions while fixing the other distributions and setting $\partial F\left(q\right)/\partial q(·)=0$. Furthermore, the computation of $\partial F\left(q\right)/\partial q(·)=0$ (Assuming that $q\left(\mathbf{Z}\right)={\displaystyle \prod _{i=1}^K}q\left({\mathbf{z}}_i\right)$) can be expressed as

$$q\left({\mathbf{z}}_j\right)\propto exp\left({〈lnp(\mathbf{Y},\mathbf{Z};\mathbf{\theta })〉}_{i\ne j}\right),$$

(41)

where ${〈·〉}_{i\ne j}$ denotes an expectation with respect to all factors except $q\left({\mathbf{z}}_j\right)$.

Next, the approximate posterior distribution of each part in Equation (40) is calculated according to Equation (41). We first calculate the variational distribution for common component ${\mathbf{z}}_c$, via Equation (42).

$$\begin{array}{c}q\left({\mathbf{z}}_c\right)\propto exp\left(\langle lnp\left(\mathbf{Y},\mathbf{Z};\mathbf{\theta }\right){\rangle }_{q\left(\mathbf{Z}\right)/q\left({\mathbf{z}}_c\right)}\right)\propto exp\left(\sum _{k=1}^K\langle lnp\left({\mathbf{y}}_k|{\mathbf{z}}_c,{\mathbf{z}}_k,\beta \right){\rangle }_{q\left({\mathbf{z}}_k\right)q\left(\beta \right)}+\langle lnp\left({\mathbf{z}}_c|{\mathbf{\alpha }}_c\right){\rangle }_{q\left({\mathbf{\alpha }}_c\right)}\right)\hfill \\ \propto \mathcal{N}\left({\mathbf{z}}_c;{\mathbf{\mu }}_c,{\mathsf{\Sigma }}_c\right).\hfill \end{array}$$

(42)

Thus ${\mathbf{z}}_c$ is confirmed to be Gaussian distributed, with covariance matrix ${\mathsf{\Sigma }}_c$ and mean vector ${\mathbf{\mu }}_c$, where

$$\begin{array}{c}{\mathsf{\Sigma }}_c={\left({\displaystyle \sum _{k=1}^K}\langle \beta \rangle {\mathbf{A}}_k^T{\mathbf{A}}_k+{\mathsf{\Lambda }}_{\langle {\mathbf{\alpha }}_c\rangle }\right)}^{-1},\hfill \end{array}$$

(43)

$$\begin{array}{c}{\mathbf{\mu }}_c=\langle \beta \rangle {\mathsf{\Sigma }}_c{\displaystyle \sum _{k=1}^K}{\mathbf{A}}_k^T({\mathbf{y}}_k-{\mathbf{A}}_k{\mathbf{\mu }}_k).\hfill \end{array}$$

(44)

In the above equations, ${\mathsf{\Lambda }}_{\langle {\alpha }_c\rangle }\in {\Re }^{N\times N}$ is a diagonal matrix with hyperparameters ${\alpha }_c^n$ ($n=1,\dots ,N$), and K is the number of channels.

Secondly, we calculate the variational distributions for innovation components ${\mathbf{z}}_k,k=1,\cdots ,K$ via Equation (45).

$$\begin{array}{cc}q\left({\mathbf{z}}_k\right)& \propto exp\left({〈lnp\left(\mathbf{Y},\mathbf{Z};\mathbf{\theta }\right)〉}_{q\left(\mathbf{Z}\right)/q\left({\mathbf{z}}_k\right)}\right)\propto exp\left({〈lnp\left({\mathbf{y}}_k|{\mathbf{z}}_c,{\mathbf{z}}_k,\beta \right)〉}_{q\left({\mathbf{z}}_c\right)q\left(\beta \right)}+{〈lnp\left({\mathbf{z}}_k|{\mathbf{\alpha }}_k\right)〉}_{q\left({\mathbf{\alpha }}_k\right)}\right)\hfill \\ & {\displaystyle \propto exp\left({\displaystyle -\frac{\langle \beta \rangle }{2}{〈{∥{\mathbf{y}}_k-{\mathbf{A}}_k\left({\mathbf{z}}_c+{\mathbf{z}}_k\right)∥}_2^2〉}_{q\left({\mathbf{z}}_c\right)}-\frac{1}{2}{\mathbf{z}}_k^T{\mathsf{\Lambda }}_{\langle {\mathbf{\alpha }}_k\rangle }{\mathbf{z}}_k}\right)\propto \mathcal{N}\left({\mathbf{z}}_k;{\mathbf{\mu }}_k,{\mathsf{\Sigma }}_k\right).}\hfill \end{array}$$

(45)

Thus, ${\mathbf{z}}_k$ is also Gaussian distributed, with covariance matrix ${\mathsf{\Sigma }}_k$ and mean vector ${\mathbf{\mu }}_k$, where

$$\begin{array}{c}{\mathsf{\Sigma }}_k={\left(\langle \beta \rangle {\mathbf{A}}_k^T{\mathbf{A}}_k+{\mathsf{\Lambda }}_{\langle {\mathbf{\alpha }}_k\rangle }\right)}^{-1},\hfill \end{array}$$

(46)

$$\begin{array}{c}{\mathbf{\mu }}_k=\langle \beta \rangle {\mathsf{\Sigma }}_k{\mathbf{A}}_k^T\left({\mathbf{y}}_k-{\mathbf{A}}_k{\mathbf{\mu }}_c\right).\hfill \end{array}$$

(47)

In the above equations, ${\mathsf{\Lambda }}_{\langle {\mathbf{\alpha }}_k\rangle }\in {\Re }^{N\times N}$ is a diagonal matrix with hyperparameters ${\alpha }_k^n$ ($k=1,\dots ,K$; $n=1,\dots ,N$).

Thirdly, we calculate the variational distribution for the prior of common component.

$$\begin{array}{cc}q\left({\alpha }_c^n\right)& \propto exp\left({〈lnp\left(z_c^n|{\alpha }_c^n\right)+lnp\left({\alpha }_c^n;a,b\right)〉}_{q\left(z_c\right)q\left(\beta \right)}\right)\hfill \\ & \propto exp\left({\displaystyle \frac{1}{2}ln{\alpha }_c^n-\frac{1}{2}〈{\left(z_c^n\right)}^2〉{\alpha }_c^n+\left(a-1\right)ln{\alpha }_c^n-b{\alpha }_c^n}\right)\hfill \\ & \propto exp\left({\displaystyle \left(\frac{1}{2}+a-1\right)ln{\alpha }_c^n-\left(b+\frac{1}{2}〈{\left(z_c^n\right)}^2〉\right){\alpha }_c^n}\right)\propto Gamma\left({\alpha }_c^n;\tilde{a},{\tilde{b}}_n\right).\hfill \end{array}$$

(48)

Thus ${\alpha }_c^n$ is distributed as $Gamma({\alpha }_c^n;\tilde{a},{\tilde{b}}_n)$, where

$$\begin{array}{c}{\displaystyle \tilde{a}=a+\frac{1}{2},{\tilde{b}}_n=b+\frac{1}{2}〈{\left(z_c^n\right)}^2〉.}\hfill \end{array}$$

(49)

Similarly, we can obtain the variational distributions for the prior of innovation components, as

$$q\left({\alpha }_k^n\right)\propto Gamma({\alpha }_k^n;{\tilde{c}}_k,{\tilde{d}}_k^n),$$

(50)

where

$$\begin{array}{c}{\displaystyle {\tilde{c}}_k=c_k+\frac{1}{2},{\tilde{d}}_k^n=d_k+\frac{1}{2}〈{\left(z_k^n\right)}^2〉.}\hfill \end{array}$$

(51)

Finally, we calculate the approximate posterior distribution for the parameter of inverse variance $\beta$.

$$\begin{array}{cc}q\left(\beta \right)& \propto exp\left({〈lnp\left(\mathbf{Y},\mathbf{Z};\mathbf{\theta }\right)〉}_{q\left(\mathbf{Z}\right)/q\left(\beta \right)}\right)\propto exp\left(\sum _{k=1}^K{〈lnp\left({\mathbf{y}}_k|{\mathbf{z}}_c,{\mathbf{z}}_k,\beta \right)〉}_{q\left({\mathbf{z}}_c\right)q\left({\mathbf{z}}_k\right)}+lnp\left(\beta |e,f\right)\right)\hfill \\ & \propto exp\left({\displaystyle \frac{KM}{2}ln\beta -\frac{\beta }{2}\sum _{k=1}^K{〈{∥{\mathbf{y}}_k-{\mathbf{A}}_k\left({\mathbf{z}}_c+{\mathbf{z}}_k\right)∥}_2^2〉}_{q\left({\mathbf{z}}_c\right)q\left({\mathbf{z}}_k\right)}+\left(e-1\right)ln\beta -f\beta }\right)\hfill \\ & \propto Gamma\left(\beta ;\tilde{e},\tilde{f}\right).\hfill \end{array}$$

(52)

Thus $\beta$ is distributed as $Gamma(\beta ;\tilde{e},\tilde{f})$, where

$$\begin{array}{c}{\displaystyle \tilde{e}=e+\frac{KM}{2},}\hfill \\ {\displaystyle \tilde{f}=f+\frac{1}{2}\left(\sum _{k=1}^K{〈{∥{\mathbf{y}}_k-{\mathbf{A}}_k({\mathbf{z}}_c+{\mathbf{z}}_k)∥}_2^2〉}_{q\left({\mathbf{z}}_c\right)q\left({\mathbf{z}}_k\right)}\right).}\hfill \end{array}$$

(53)

The variational optimization proceeds by iteratively updating Equations (42), (45), (48), (50), and (52) until convergence occurs to hyperparameters $\mathbf{\theta }$. Finally, we can obtain the reconstructed signal by applying the maximum a posteriori estimation.

$$\begin{array}{c}\hfill {\widehat{\mathbf{\rho }}}_k=\underset{{\mathbf{z}}_c+{\mathbf{z}}_k}{argmax}p\left(\mathbf{Z}\right|\mathbf{Y};\mathbf{\theta })=\underset{{\mathbf{z}}_c}{argmax}q\left({\mathbf{z}}_c\right)+\underset{{\mathbf{z}}_k}{argmax}q\left({\mathbf{z}}_k\right)={\mathbf{\mu }}_c+{\mathbf{\mu }}_k,k=1,2,\cdots ,K.\end{array}$$

The proposed HVB-DCS algorithm for solving the JSM-1 reconstruction problem is summarized in Algorithm 1.

Algorithm 1: HVB-DCS Algorithm

Input: A set of measurement vectors $\{{\mathbf{y}}_1,{\mathbf{y}}_2,\cdots ,{\mathbf{y}}_K\}$ and corresponding measurement matrices

$\{{\mathbf{A}}_1,{\mathbf{A}}_2,\cdots ,{\mathbf{A}}_K\},k=1,\cdots ,K$.

Output: The reconstructed signal ${\widehat{\mathbf{\rho }}}_k={\mathbf{\mu }}_c+{\mathbf{\mu }}_k$, $k=1,\cdots ,K$.

Initialize the hyperparameters. Set the initial values of the variables (a, b, $\{c_k,k=1,\cdots ,K\}$,

$\{d_k,k=1,\cdots ,K\}$, e, f) as ${10}^{-6}$.

Compute the variational distribution for the common component.

Compute ${\mathsf{\Sigma }}_c={\left({\displaystyle \sum _{k=1}^K}\langle \beta \rangle {\mathbf{A}}_k^T{\mathbf{A}}_k+{\mathsf{\Lambda }}_{\langle {\mathbf{\alpha }}_c\rangle }\right)}^{-1}$, and ${\mathbf{\mu }}_c=\langle \beta \rangle {\mathsf{\Sigma }}_c{\displaystyle \sum _{k=1}^K}{\mathbf{A}}_k^T({\mathbf{y}}_k-{\mathbf{A}}_k{\mathbf{\mu }}_k).$

Compute the variational distribution for the prior of the common component.

Update $q\left({\alpha }_c^n\right)$, compute

${\displaystyle \tilde{a}=a+\frac{1}{2},{\tilde{b}}_n=b+\frac{1}{2}〈{\left(z_c^n\right)}^2〉.}$

Compute the variational distributions for the innovation components.

Compute ${\mathsf{\Sigma }}_k={\left(\langle \beta \rangle {\mathbf{A}}_k^T{\mathbf{A}}_k+{\mathsf{\Lambda }}_{\langle {\mathbf{\alpha }}_k\rangle }\right)}^{-1}$, and ${\mathbf{\mu }}_k=\langle \beta \rangle {\mathsf{\Sigma }}_k{\mathbf{A}}_k^T\left({\mathbf{y}}_k-{\mathbf{A}}_k{\mathbf{\mu }}_c\right)$.

Compute the variational distributions for the prior of the innovation components.

Update $q\left({\alpha }_k^n\right)$, compute

${\displaystyle {\tilde{c}}_k=c_k+\frac{1}{2},{\tilde{d}}_k^n=d_k+\frac{1}{2}〈{\left(z_k^n\right)}^2〉.}$

Compute the variational distribution for the prior of noise vector.

Update $q\left(\beta \right)$, compute

${\displaystyle \tilde{e}=e+\frac{KM}{2},}$

${\displaystyle \tilde{f}=f+\frac{1}{2}\left(\sum _{k=1}^K{〈{∥{\mathbf{y}}_k-{\mathbf{A}}_k({\mathbf{z}}_c+{\mathbf{z}}_k)∥}_2^2〉}_{q\left({\mathbf{z}}_c\right)q\left({\mathbf{z}}_k\right)}\right).}$

Iterate steps 2 , 3 , 4 , 5 and 6 until convergence occurs to hyperparameters.

Output ${\widehat{\mathbf{\rho }}}_k={\mathbf{\mu }}_c+{\mathbf{\mu }}_k$ for $k=1,\cdots ,K$.

5.2. Complexity Analysis

We present the computational complexity of the proposed HVB-DCS algorithm, as well as the comparison to the CS-based joint sparsity recovery algorithm and the RD algorithm. The computational complexity of the HVB-DCS algorithm at each iteration is dominated by the inversion operations on an $N\times N$ matrix in Equations (43) and (46), and the multiplications of a matrix and a vector in Equations (44) and (47). The computational complexities of the matrix inversion and matrix-vector multiplication are $\mathcal{O}\left(N^3\right)$ and $\mathcal{O}\left(N^2\right)$, respectively. By using matrix inversion lemma [39], the complexity of the matrix inversion can be reduced to $\mathcal{O}\left(M^3\right)$. Thus the overall computational complexity of the HVB-DCS algorithm is $\mathcal{O}\left(N_t\left(K+1\right)M^3\right)$, where $N_t$ is the total number of iterations. The small values of $N_t$ and M will result in acceptable computational complexity. In contrast, a large N will lead to high computational complexity, which makes the HVB-DCS algorithm impractical. To circumvent the high computational complexity problem, we can resort to a number of accelerating algorithms [40,41] to develop a more computationally efficient algorithm in the future work.

For the CS-based joint sparsity recovery algorithm, the dimensions of the sparse vector and the measurement vectors are $3N$ and $2M$, respectively. Thus the computational complexity of the CS-based algorithm is $\mathcal{O}\left(12N{M}^2\right)$. For the RD algorithm, the computational complexity is dominated by the fast Fourier transform (FFT) and inverse fast Fourier transform (IFFT), which have computational complexities of $\mathcal{O}\left(Mlo{g}_{2}M\right)$ and $\mathcal{O}(M/2lo{g}_{2}M)$, respectively. Thus, the overall computational complexity of the RD algorithm is $\mathcal{O}\left(K\left(3N/2+1\right)Mlo{g}_{2}M\right)$.

6. Simulations and Experiments

In this section, the performance of the proposed HVB-DCS algorithm is verified through both simulations and experiments. First, the effects of along-track and cross-track velocities on the SAR images reconstructed using the HVB-DCS algorithm are analyzed in Section 6.1. Secondly, the proposed HVB-DCS algorithm is compared with some classical GMTI algorithms, e.g., the RD [42] and CS [43] algorithms, based on point target simulation, in Section 6.2. Thirdly, the effects of problem size and data rate on complexity of the HVB-DCS algorithm are analyzed in Section 6.3. Finally, in Section 6.4, the proposed HVB-DCS algorithm is applied to the real background cluttered environment, where the real SAR data are collected by RADARSAT-1 satellite in the Vancouver region [42].

6.1. The Effects of Along-Track and Cross-Track Velocities on SAR Imaging

6.1.1. The Effect of Cross-Track Velocity

First, we consider the influence of cross-track velocity on SAR imaging and take channel one as an example. The parameters of the simulated SAR radar system are listed in

Table 1. SAR radar system parameters for simulation.

Parameter	Value
Wavelength (m)	0.03
Range bandwidth (MHz)	150
Platform height (m)	5000
Platform velocity (m/s)	150
Incidence angle (${}^{\circ }$)	45
PRF (Hz)	300
Aperture size (m)	2
Channel distance (m)	1

. Suppose that a moving target is at $(R_B,x_0)$, with a given cross-track velocity $v_r$ and along-track velocity $v_a=0$ m/s. The received signal of the moving target by channel one, after demodulation and range compression, which is presented in Equation (13), is rewritten as

$$\begin{array}{c}\hfill \begin{array}{c}\hfill s_1\left(t_m\right)={\rho }_0^{\left(1\right)}{\omega }_a\left(t_m-k\Delta \tau \right)exp\left(-j2\pi \frac{v^2}{\lambda R_B}{\left(t_m-k\Delta \tau -\frac{v_r{R}_B}{v^2}\right)}^2\right),\end{array}\end{array}$$

(54)

where $k\Delta \tau =x_0/v$. Since the CS-based SAR imaging algorithm applied in the azimuth direction is essentially searching the most matched atom from a redundant dictionary Equation (15) to $s_1\left(t_m\right)$, the degree of match is mainly affected by the exponential term in Equation (54), whereas the minor shift of the azimuth windowing function has little influence. Assuming that the most matched atom in the redundant dictionary ${\mathsf{\Phi }}_1$ is

$$\begin{array}{c}{s}_0(t_m-k^{\prime }\Delta \tau )={\omega }_a(t_m-k^{\prime }\Delta \tau )exp\left(-j2\pi \frac{v^2}{\lambda R_B}{(t_m-k^{\prime }\Delta \tau )}^2\right).\hfill \end{array}$$

(55)

When ${\displaystyle k^{\prime }\Delta \tau =k\Delta \tau +\frac{v_r{R}_B}{v^2}}$, $s_1\left(t_m\right)$ matches the atom most, resulting in ${\displaystyle k^{\prime }=k+\frac{v_r{R}_B}{v^2\Delta \tau }}$. Obviously, for a stationary clutter with $v_r=0$ m/s, the reflection coefficient of the stationary clutter locates at the ${(N/2+k)}^{th}$ pixel in the azimuth direction. For a moving target with $v_r\ne 0$, the cross-track velocity of the moving target leads to an azimuth deviation of ${\displaystyle \frac{v_r{R}_B}{v^2\Delta \tau }}$ pixels, compared with that of the stationary clutter, on the SAR image reconstructed using the HVB-DCS algorithm.

The influence of cross-track velocity is demonstrated in Figure 5. In Figure 5, a static scattering center and a moving scattering center with the cross-track velocity of $v_r=0.5$ m/s and the along-track velocity of $v_a=0$ m/s are both at $(7071,0)$ m. The amplitude of reflection coefficient of the moving target is half that of the stationary clutter. The moving target has a deviation of 47 pixels in the azimuth direction compared with that of the stationary clutter on SAR image, which completely conforms to ${\displaystyle \frac{v_r{R}_B}{v^2\Delta \tau }}$.

6.1.2. The Effect of Along-Track Velocity

For the analysis on along-track velocity, it is assumed that the moving target is at $(R_B,0)$, with a given along-track velocity $v_a$ and cross-track velocity of $v_r=0$ m/s. Thus the received signal of the moving target, after demodulation and range compression, can be written as

$$\begin{array}{c}\hfill \begin{array}{c}\hfill s_1\left(t_m\right)={\rho }_0^{\left(1\right)}{\omega }_a\left(t_m\right)exp\left(-j2\pi \frac{{\left(v-v_a\right)}^2}{\lambda R_B}{t}_m^2\right).\end{array}\end{array}$$

(56)

Similarly, for channel one, the received signal after range compression can be represented linearly by the atoms in the redundant dictionary ${\mathsf{\Phi }}_1$, as

$$\begin{array}{c}{s}_1\left(t_m\right)={\displaystyle \sum _{i=-N/2+1}^{N/2}}{C}_i^{\left(1\right)}{\omega }_a\left(t_m-i\Delta \tau \right)exp\left(-j2\pi \frac{v^2}{\lambda R_B}{\left(t_m-i\Delta \tau \right)}^2\right),\hfill \end{array}$$

(57)

where $C_i^{\left(1\right)}$ denotes the reconstructed reflection coefficient in channel one. Substituting Equation (56) into Equation (57), and dividing both sides by the item ${\displaystyle exp\left(-j2\pi \frac{v^2}{\lambda R_B}{t}_m^2\right)}$, results in

$$\begin{array}{c}{\rho }_0^{\left(1\right)}{\omega }_a\left(t_m\right)exp\left({\displaystyle -j2\pi \frac{\left(v_a^2-2v{v}_a\right)}{\lambda R_B}{t}_m^2}\right)={\displaystyle \sum _{i=-N/2+1}^{N/2}}{D}_i^{\left(1\right)}{\omega }_a\left(t_m-i\Delta \tau \right)exp\left(-ji\omega t_m\right),\hfill \end{array}$$

(58)

where ${\displaystyle \omega =\frac{-4\pi v^2}{\lambda R_B}\Delta \tau }$, and $D_i^{\left(1\right)}=C_i^{\left(1\right)}exp\left({\displaystyle -j\frac{2\pi v^2}{\lambda R_B}{(i\Delta \tau )}^2}\right)$. When the along-track velocity $v_a$ of the moving target equals zero, from Equation (58), we can see that the left side of the equation equals a constant. Thus on the right side of the equation, the most sparse solution is $D_i^{\left(1\right)}=0(i\ne 0)$ and $D_0^{\left(1\right)}={\rho }_0^{\left(1\right)}$. In general, the along-track velocity of the target is much less than that of the flying platform, resulting in a slow variation of the terms on the left side of Equation (58). Thus the left side of Equation (58) contains only a few low frequency components. With the increase of left side of Equation (58), the high frequency component increases slowly. $D_i^{\left(1\right)}$ can be approximated to the coefficients of expansion of the slowly varying function under different frequency components. Therefore, with the increase of the left side of the Equation (58), the nonzero value $D_i^{\left(1\right)}$ will spread to both sides, centering at $D_0^{\left(1\right)}$, which leads to the defocus of the moving target on the SAR image. The influence of along-track velocity can be shown in Figure 6, where $v_a$ is set as different values of 0 m/s, 1 m/s and 2 m/s. It can be seen from Figure 6 that as the along-track velocity increases, the azimuth defocus of the moving target on the SAR image becomes more serious.

The defocus of the moving target on the SAR image can be treated as a basis mismatch problem [44]. In this work, the measurement matrices are built based on the echo signal models of two channels (Equations (13) and (19) in Section 4.1), assuming that the along-track velocity $v_a$ equals zero. However, in practice, this assumption does not always hold true. In this case, the constructed measurement matrix does not exactly match the underlying true one, leading to the basis mismatch problem. A parameter ${\beta }_d$ is defined in [44] to quantitatively evaluate the difference between the assumed basis and the true one. The larger the absolute value of $v_a$ is, the higher the mismatch level of the constructed measurement matrix will be, resulting in a higher ${\beta }_d$. Moreover, the defocus effect can be evaluated by the reconstruction error of the sparse vector recovered with mismatched basis. Based on [44], an upper bound of the ${\ell }_1$-norm reconstruction error is defined as:

$${\parallel {\rho }_k-{\widehat{\rho }}_k\parallel }_1\le N{\beta }_d{\parallel {\rho }_k\parallel }_q,$$

(59)

where ${\rho }_k$ and ${\widehat{\rho }}_k$ denote the true and reconstructed sparse vector respectively, N is the dimension of the sparse vector, ${\parallel ·\parallel }_q$ represents the ${\ell }_q$-norm, and $1\le q\le \infty$. Thus, the azimuth defocus on the SAR image is bounded given $v_a$. In addition, please refer to [45,46] for more analyses of basis mismatch problem in CS-based SAR imaging and the methods to deal with it.

Furthermore, the effect of $v_a$ on the azimuth defocus is analyzed through simulation. The degree of azimuth defocus is represented by the relative reconstruction error, which is defined as

$$err\triangleq {\parallel {\rho }_k-{\widehat{\rho }}_k\parallel }_1/{\parallel {\rho }_k\parallel }_1.$$

(60)

Figure 7 shows the relative reconstruction error $err$ versus $v_a$. The result illustrates that the reconstruction error $err$ is low when the absolute value of $v_a$ is relatively small.

6.2. The Simulation on Point Targets

In this section, the performance of the HVB-DCS algorithm is verified by simulations on point targets. The nearest slant range $R_B$ is set as 7071 m. Four scattering centers are located at the range bin $R_B=7071$ m, including three static scattering centers, and a moving scattering center with cross-track velocity of $v_r=0.5$ m/s and along-track velocity of $v_a=0$ m/s. The three stationary targets are located at $-5$ m, 0 m, $+5$ m in the azimuth direction, and the moving target is located at 0 m in the azimuth direction. The amplitudes of reflection coefficient of three stationary targets are the same, which are twice that of the moving target. The parameters of the simulated SAR radar system are the same as those listed in Table 1.

In the field of SAR-GMTI, the proposed HVB-DCS algorithm is compared to the RD [42] and CS [43] algorithms. The flow charts of the RD based GMTI system, the CS-based GMTI system, and the HVB-DCS based GMTI system are shown in Figure 8, Figure 9 and Figure 10, respectively. For the RD based GMTI system and CS-based GMTI system, an additional DPCA [47] procedure is needed for moving-target indication, which is shown in both Figure 8 and Figure 9. DPCA can effectively suppress clutter and detect slow moving targets. In contrast, in the HVB-DCS based GMTI system, the moving targets (the innovation components) and the stationary targets (the common component) can be directly distinguished using the DCS framework. Thus the DPCA is no longer required in the HVB-DCS based GMTI system and removed (Figure 10). This significantly simplifies the architecture of the GMTI system, and saves the computing time.

Figure 8, Figure 9 and Figure 10 show that all the three algorithms, i.e., RD, CS and HVB-DCS algorithms, can effectively identify the moving target. For the RD algorithm, we use the oversampled raw data to image. The result of the RD algorithm is shown in Figure 11. As a widely used SAR imaging algorithm, the RD algorithm has very high running speed, but the on-board memory required by the sampling rate is large. For the CS algorithm, the simulated radar system works in the random sampling mode introduced in Section 4.1. In our experiment, we choose randomly half of transmitted pulses in azimuth and recover the complex-valued SAR image based on the CS algorithm. The result is shown in Figure 12, which shows that the image reconstructed by using the CS algorithm is more clear, and the side-lobes are effectively suppressed. However, the computational cost of the CS algorithm is relatively high. For the HVB-DCS algorithm, we further drop randomly $62.5\%$ of the transmitted pulses in azimuth to recover the SAR image. In our experiment, we reduce to about $37.5\%$ of the original samples in the azimuth direction and reconstruct the SAR image using the HVB-DCS algorithm. The result is shown in Figure 13, which shows that the HVB-DCS algorithm has a great suppression effect on stationary clutter and is capable of identifying the moving target accurately. The effects of both the data rate in azimuth, and noise/clutter level on the detection performance of moving targets are tested by simulations. The adopted metric is the improvement factor (IF), which is defined as IF = ${\mathrm{SCNR}}_{\mathrm{out}}$/${\mathrm{SCNR}}_{\mathrm{in}}$, where ${\mathrm{SCNR}}_{\mathrm{in}}$ and ${\mathrm{SCNR}}_{\mathrm{out}}$ are the signal-to-clutter-noise ratios of input signal and output signal, respectively. We test the HVB-DCS and CS algorithms to detect moving targets by varying the data rate in the azimuth direction, at different ${\mathrm{SCNR}}_{\mathrm{in}}$ levels, except that the performance of RD algorithm is tested at different ${\mathrm{SCNR}}_{\mathrm{in}}$ levels with data rate $100\%$. The IFs averaged over 100 Monte Carlo trials are shown in Figure 14. For the HVB-DCS algorithm, at the same ${\mathrm{SCNR}}_{\mathrm{in}}$ level, the IF first increases dramatically with the data rate in azimuth. When the data rate in the azimuth direction exceeds $37.5\%$, the IF grows slowly. For the CS algorithm, at the clutter/noise levels of ${\mathrm{SCNR}}_{\mathrm{in}}=-11$ dB and ${\mathrm{SCNR}}_{\mathrm{in}}=-6$ dB, the IF is negative infinite when the data rate in azimuth is lower than $37.5\%$. Furthermore, at the clutter/noise level of ${\mathrm{SCNR}}_{\mathrm{in}}=-2$ dB, the IFs are negative infinite when the data rate in azimuth is lower than $25\%$. With the data rate in azimuth approaching $37.5\%$, the IF increases to a value larger than zero. When the data rate in the azimuth direction increases to more than $50\%$, the IF gradually increases with the data rate in azimuth at the same ${\mathrm{SCNR}}_{\mathrm{in}}$ level. The main reason is that when the data rate in azimuth decreases, the reconstruction error increases gradually, thus affecting the performance of clutter suppression. Moreover, it is noted that the IFs of HVB-DCS algorithm with the data rate $37.5\%$ in azimuth are almost at the same level with those of CS algorithm with the data rate $50\%$ in azimuth. The IFs of RD algorithm (indicated by asterisks) are the lowest among the three algorithms. This verifies that the proposed HVB-DCS algorithm outperforms the CS and RD algorithms in terms of clutter suppression and moving target detection performance.

To examine the reconstruction accuracy of the proposed HVB-DCS algorithm, we introduce the reconstruction error, which is defined as

$${\displaystyle e_{rec}=\frac{{\sum }_{l=1}^L{∥{\widehat{x}}_l-x_l∥}_2}{{\sum }_{l=1}^L{∥x_l∥}_2},}$$

(61)

where L is the number of range bins, $x_l$ and ${\widehat{x}}_l$ denote the true and reconstructed sparse vector at the $l^{th}$ range bin, respectively. The reconstruction errors of the RD, CS and HVB-DCS algorithms, are calculated by varying data rates in azimuth, at different ${\mathrm{SCNR}}_{\mathrm{in}}$ levels. The results are shown in Figure 15. The reconstruction error of the HVB-DCS algorithm is lower than $0.2$ when the data rate is not less than $37.5\%$, whereas the reconstruction error of the CS algorithm is lower than $0.3$ when the data rate is not less than $50\%$. In addition, compared to the CS algorithm, the HVB-DCS algorithm has smaller reconstruction errors when the data rate exceeds $37.5\%$. The main reasons are due to the following two aspects. First, the DCS exploits both intra-signal and inter-signal correlations, which ensures a good correlation between the channels in the case of lower data rate. Secondly, the hierarchical structure is more suitable for complex models. Moreover, the higher the noise/clutter level is, the larger the reconstruction errors are, which shows that the detection performance degrades with the increase of the noise/clutter level.

6.3. Effects of Problem Size and Data Rate on Complexity of HVB-DCS

The goal of this section is to analyze the computational complexity of the HVB-DCS algorithm in terms of problem size and data rate using simulations. The computational complexity of the HVB-DCS algorithm is further compared with those of the RD and CS-based algorithms. First, we set the data rate as $37.5\%$ and $50\%$ for the HVB-DCS and CS-based algorithms respectively, and $100\%$ for the RD algorithm. Meanwhile, we vary the dimension of sparse vector N (also, the number of samples in the azimuth direction) in the range of $[100,900]$. We plot the runtime (averaged over 100 trials) of the HVB-DCS, RD and CS-based algorithms versus N in Figure 16a. The results show that the RD algorithm achieves the shortest runtime. When the data rate is at the level of $37.5\%$ and $N\le 500$, the proposed HVB-DCS algorithm takes almost the same runtime with the CS-based algorithm to perform moving target indication. When the data rate is at the level of $37.5\%$ and $N>500$, the proposed HVB-DCS algorithm takes longer runtime than the CS-based algorithm. When the data rate is at the level of $50\%$, the runtime of HVB-DCS is shorter than that of the CS-based algorithm in the range of $[100,900]$.

In addition, we examine the runtime of the HVB-DCS and CS-based algorithms versus data rate when setting N at different levels, i.e., 100, 500 and 900 respectively. The data rate varies in the range of $[12.5\%,100\%]$ for the two algorithms. The results are shown in Figure 16b. When $N=100$, the runtime of HVB-DCS is shorter than that of the CS-based algorithm at any data rate in the range of $[12.5\%,100\%]$. When $N=500$ and the data rate is larger than $37.5\%$, HVB-DCS takes shorter runtime than that of the CS-based algorithm, and the similar result occurs when $N=900$ and the data rate is larger than $50\%$.

6.4. Experiment on Real SAR Data

In order to verify the performance of the proposed HVB-DCS algorithm in practice, we carry out an experiment on data from a real background cluttered environment, which are collected by the satellite RADARSAT-1 in the Vancouver region [42]. The parameters of RADARSAT-1 satellite are shown in

Table 2. RADARSAT-1 parameters.

Parameter	Value
Slant range of scene center (km)	$988.65$
Beam squint angle (rad)	$0.0279$
Effective radar velocity (m/s)	7062
PRF (Hz)	$1256.98$
Sampling rate (MHz)	$32.317$
Range FM rate (MHz/$\mathsf{\mu }$s)	$0.72135$
Pulse duration ($\mathsf{\mu }$s)	$41.75$
Radar center frequency (MHz)	5300

.

In the experiment, the real SAR data from single-channel SAR system (RADARSAT-1) is used to simulate dual-channel raw data. The detailed procedures are as follows. According to the imaging geometry models for point targets as in Equations (13) and (19) in Section 4.1, there is a phase difference $\exp \left(-j\pi d^2/\left(2\lambda R_B\right)\right)$ and a time delay $d/\left(2\left(v-v_a\right)\right)$ between the received signals from two channels. First, all the real data from the first channel is delayed $\Delta t$ time, and then multiplied with $\exp \left(-j\pi d^2/\left(2\lambda R_B\right)\right)$, to produce the raw data from the second channel of the simulated dual-channel SAR system. Secondly, in order to simulate the channel imbalance, some random amplitudes and phase errors are added to the raw data of the second channel. Finally, the echo signals of nine targets (including five stationary targets and four moving targets moving) are added to the raw data of both channels. Based on the raw data from the simulated dual-channel SAR system, three algorithms, e.g., the RD, CS and HVB-DCS algorithms, are utilized to detect the moving targets.

Figure 17 shows the SAR image of the observation scene, which contains the static clutter, five stationary point targets and, four moving targets with the same cross-track velocity of $v_r=5$ m/s, and the same along-track velocity of $v_a=5$ m/s. The radar image covers an area of about $2.8$ km in azimuth and $1.9$ km in range. The corresponding parameters of nine targets are listed in

Table 3. Simulation parameters of five stationary targets and four moving targets.

No.	Azimuth Coordinate (m)	Nearest Slant Range (km)	Along-Track Velocity ${\mathit{v}}_{\mathit{a}}$ (m/s)	Across-Track Velocity ${\mathit{v}}_{\mathit{r}}$ (m/s)
1	$-500$	$0.8$	5	5
2	$-400$	$987.65$	0	0
3	$800$	$987.65$	0	0
4	$200$	$988.17$	0	0
5	$-1100$	$988.17$	5	5
6	$100$	$988.17$	5	5
7	$-500$	$988.73$	5	5
8	$-400$	$988.73$	0	0
9	$800$	$988.73$	0	0

. The four moving targets in full SAR RADARSAT-1 image are marked with the red rectangles in Figure 17. We use the oversampled raw data in azimuth for the RD algorithm, and $50\%$ of the oversampled raw data in azimuth for CS algorithm and $37.5\%$ of the oversampled raw data in azimuth for the HVB-DCS algorithm.

GMTI results (Figure 18) show that all of the three algorithms can distinguish the moving targets from the stationary clutter successfully. Moreover, the final ${\mathrm{SCNR}}_{\mathrm{out}}$ of the three algorithms are calculated to evaluate detection performance. The initial ${\mathrm{SCNR}}_{\mathrm{in}}$ equals $-69$ dB. The ${\mathrm{SCNR}}_{\mathrm{out}}\mathrm{s}$ of the SAR image recovered based on the RD, CS and HVB-DCS algorithms are $-13$ dB, $-12$ dB and $-5$ dB, respectively. Note that, there is an obvious increase in the $\mathrm{SCNR}$. The RD algorithm has the lowest ${\mathrm{SCNR}}_{\mathrm{out}}$. The SAR image of CS-based algorithm has a much lower ${\mathrm{SCNR}}_{\mathrm{out}}$ ($-12$ dB) compared with the proposed HVB-DCS algorithm ($-5$ dB). This is due to the reason that in the CS-based GMTI system, the mismatch between channels leads to the failure in energy cancelation of static targets and clutters, whereas the proposed HVB-DCS algorithm takes advantage of the correlation between two channels and suppresses the clutter and noise efficiently, even at a much lower data rate in the azimuth direction.

7. Conclusions

In this paper, we propose a novel HVB-DCS algorithm for GMTI using the dual-channel SAR system. The results from both the simulation on point targets and the experiment on data from the real background cluttered environment show that the proposed HVB-DCS algorithm can successfully detect multiple moving targets, while suppressing the clutter efficiently. The proposed algorithm achieves better detection performance for the GMTI system, in the metric of IF and reconstruction error, compared with the RD and CS algorithms.