## *Review* **Compressed Sensing Radar Imaging: Fundamentals, Challenges, and Advances**

#### **Jungang Yang, Tian Jin, Chao Xiao \* and Xiaotao Huang**

College of Electronic Science and Technology, National University of Defense Technology, Changsha 410073, China

**\*** Correspondence<EMAIL_ADDRESS>
Received: 3 June 2019; Accepted: 11 July 2019; Published: 13 July 2019

**Abstract:** In recent years, sparsity-driven regularization and compressed sensing (CS)-based radar imaging methods have attracted significant attention. This paper provides an introduction to the fundamental concepts of this area. In addition, we will describe both sparsity-driven regularization and CS-based radar imaging methods, along with other approaches in a unified mathematical framework. This will provide readers with a systematic overview of radar imaging theories and methods from a clear mathematical viewpoint. The methods presented in this paper include the minimum variance unbiased estimation, least squares (LS) estimation, Bayesian maximum a posteriori (MAP) estimation, matched filtering, regularization, and CS reconstruction. The characteristics of these methods and their connections are also analyzed. Sparsity-driven regularization and CS based radar imaging methods represent an active research area; there are still many unsolved or open problems, such as the sampling scheme, computational complexity, sparse representation, influence of clutter, and model error compensation. We will summarize the challenges as well as recent advances related to these issues.

**Keywords:** radar imaging; synthetic aperture radar; compressed sensing; sparse reconstruction; regularization

#### **1. Introduction**

Radar imaging technique goes back to at least the 1950s. In the past 60 years, it has been stimulated by hardware performance, imaging theories, and signal processing technologies. Figure 1 shows the developmental history of radar imaging methods.

**Figure 1.** Developmental history of radar imaging methods.

Since the development of radar imaging techniques, the main theory that has been used has always been matched filtering [1–3]. Matched filtering is a linear process; it has the advantages of simplicity and stability. However, the drawbacks of the matched filtering method are also obvious. Since it does not exploit any prior information concerning the expected targets, its performance is limited by the signal bandwidth. It also requires a dense sampling to record the signals, according to the Shannon–Nyquist sampling theorem. Thus, the matched filtering method places significant requirements on the measured data, but only produces results with limited performance. As higher and higher imaging performance is demanded, the matched filtering method will struggle to meet the requirements.

Apart from the matched filtering framework, from a more generic mathematical viewpoint, radar imaging can be viewed as an inverse problem [4–7], whereby a spatial map of the scene is recovered using the measurements of the scattered electric field. The radar observation process is a Fredholm integral (F-I) equation of the first kind [8]. Due to observation limitations, such as limited bandwidth and limited observation angles, this inverse problem is usually ill-posed [9,10]. The classic least squares (LS) estimation method cannot solve such ill-posed inverse problems efficiently. The matched filtering method can be viewed as using an approximation to eliminate the irreversible or unstable term in the LS solution. This approximation leads to limited resolution and side-lobes in the results. Thus, matched filtering methods typically provide an image that blurs the details of the scene. Using proper models for the targets, super-resolution methods can improve the resolution of the imaging result [11,12].

Besides using approximation, the ill-posed inverse problem can be solved by another approach, i.e., adding an extra constraint to the LS formula and yielding a stable solution. This approach is called regularization [8]. In order to make the solution after regularization closer to the true value, the additional constraint should represent appropriately some prior knowledge. The regularization approach can also be explained by the Bayesian maximum a posteriori (MAP) estimation theory [6,13,14], which uses prior knowledge in a probabilistic way.

In the radar imaging scenario, imposing sparsity is one possible form of prior knowledge [15]. The advantages of the sparsity-driven regularization methods include increased image quality and robustness to limitations in data quantity. Compressed sensing (CS) refers to the use of under-sampled measurements to obtain the coefficients of a sparse expansion [16–20].

This paper summarizes the fundamentals, challenges and recent advances of sparse regularization and CS-based radar imaging methods. Using a unified mathematical model, we derive the best estimator (i.e., the minimum variance unbiased estimator), the LS estimator, the Bayesian MAP estimator, matched filtering, regularization, and CS reconstructions of the scene. The characteristics of these methods and their connections are also analyzed. Finally, we present some key challenges and recent advances in this area. These include the sampling scheme, the computational complexity, the sparse representation, the influence of clutter, and the model error compensation.

#### **2. Mathematical Fundamentals of Radar Imaging**

#### *2.1. Radar Observation Model*

In the continuous signal domain, under the Born approximation, the radar observation process can be denoted as [4]

$$s(\mathbf{r}) = \int A(\mathbf{r}, \mathbf{r'}) g(\mathbf{r'}) d\mathbf{r'} + n \tag{1}$$

where *s*(**r**) denotes the observed data at the observation position of **r**, *g*(**r** ) denotes the reflectivity coefficient at **r** in the scene, *A*(**r**,**r** ) denotes the system response from **r** to **r**, and *n* denotes noise.

Assuming the system is shift invariant, Equation (1) can be rewritten as

$$s(\mathbf{r}) = \int A(\mathbf{r} - \mathbf{r'}) g(\mathbf{r'}) d\mathbf{r'} + n \tag{2}$$

It can be seen that the radar observation model is a convolution process. Equation (1) is a Fredholm integral (F-I) equation of the first kind [8]. From a mathematical viewpoint, radar imaging can be viewed as the solution of the F-I equation—i.e., we want to recover *g*(**r**) from the observed data *s*(**r**) using the observation equation. Unfortunately, according to the theory of integral equations, solving the F-I equation is usually an ill-posed problem [8].

In practice, since digitization is commonly used, the observed data are discrete. Based on Equation (1), the discrete observation model can be written as

$$\mathbf{s} = \mathbf{A}\mathbf{g} + \mathbf{n} \tag{3}$$

where **s** is stacked from the samples of *s*(**r**), **g** is stacked from the samples of *g*(**r** ), **A** is formed from samples of *A*(**r**,**r** ), and **n** is the observation noise vector.

#### *2.2. Best Linear Unbiased Estimate and Least Squares Estimate of the Scene*

From the observation model shown in (3), radar imaging can be viewed as an estimation problem, in which the scene **g** is estimated based on the observed data **s** in a noisy environment. According to estimation theories, the minimum variance unbiased estimate is the "best" estimate in terms of estimation square error. From Equation (3), it can be seen that when the radar observation model is linear, the minimum variance unbiased estimate is the best linear unbiased estimate [13]—i.e., the expression of the best estimate of the scene is

$$\stackrel{\star}{\mathbf{g}} = \left(\mathbf{A}^{H}\mathbf{C}^{-1}\mathbf{A}\right)^{-1}\mathbf{A}^{H}\mathbf{C}^{-1}\mathbf{s} \tag{4}$$

where **C** is the covariance matrix of the noise term (**C** = *E* & **nn***H*' ).

In practice, a more tractable approach is LS estimation, which can be denoted as

$$\hat{\mathbf{g}} = \operatorname\*{argmin}\_{\mathbf{g}} \|\mathbf{s} \cdot \mathbf{A} \mathbf{g}\|\_2^2 \tag{5}$$

Therefore, the LS estimate of the scene is

$$\hat{\mathbf{g}} = \left(\mathbf{A}^{H}\mathbf{A}\right)^{-1}\mathbf{A}^{H}\mathbf{s} \tag{6}$$

If **n** is white Gaussian noise, we have **C** = σ2**I**, where **I** is the identity matrix. Under such condition, Equations (4) and (6) are the same. Therefore, the LS estimate will equal to the best estimate in white Gaussian noise [13].

If we want to use Equation (6) to calculate the best estimate of the scene, a prerequisite is that (**A***H***A**) is invertible. However, in practice, this prerequisite is usually not satisfied, as discussed below. We assume that the size of **A** is *M* × *N*, where *M* denotes the number of measurements and *N* denotes the number of unknown grid points. Then, the size of (**A***H***A**) is *N* × *N*.

One case is that *M* < *N*, i.e., the number of measurements is less than the unknown variables. CS is a typical example of this case. In such a case, rank(**A***H***A**) = rank(**A**) ≤ *M* < *N*, i.e., (**A***H***A**) is irreversible.

In the above case, it can be seen that due to limited number of measurements, (**A***H***A**) is irreversible. Is it possible to make (**A***H***A**) invertible by increasing the number of measurements (i.e., make *M* > *N*. As mentioned previously, due to physical limitations, such as limited bandwidth and limited observation angles, if we take more measurements, the interval between the adjacent measurements will be smaller. Thus, the coherence between the adjacent columns in **A** will increase. Consequently, (**A***H***A**) −1 will probably be ill-conditioned.

In summary, the LS solution usually contains irreversible or ill-posed terms. This problem is inherent, and is derived from the property of the F-I equation of the first kind [8].

#### *2.3. Matched Filtering Method*

Examining Equation (6), it can be seen that the irreversible or ill-posed term is (**A***H***A**) −1 . We can multiply (**A***H***A**) in the left side of Equation (6) to eliminate (**A***H***A**) −1 . In this way, we can avoid explicitly calculating the nonexistent or unstable term (**A***H***A**) −1 . This leads to the matched filtering method, which can be denoted as

$$
\stackrel{\star}{\mathbf{g}}\_{\text{MF}} = (\mathbf{A}^H \mathbf{A}) \stackrel{\star}{\mathbf{g}} = \mathbf{A}^H \mathbf{s} \tag{7}
$$

Equation (7) can be viewed as multiplying the best estimate of the scene with (**A***H***A**). The matrix (**A***H***A**) is the autocorrelation of the system response, which usually has a sinc pulse shape [1,21]. The matched filtering result can be viewed as the convolution of the best estimate of the scene and the sinc function. A point target will be spread, and side-lobes will also appear in the matched filtering result [21]. This implies that the matched filtering method can only provide an image that blurs the details of the scene. The matched filtering method has a limited resolution, which depends on the autocorrelation of the system response [1].

Figure 2 shows an example of the matched filtering method. Six point targets are set in the scene. It can be seen that the matched filtering result is the convolution of the targets and the autocorrelation of the system response. As a result, an idea point target is spread into a sinc waveform. Consequently, targets will interfere with each other, and two closely spaced targets may not be resolved in the matched filtering result.

Equation (7) is the original form of the matched filtering equation. In practice, in order to reduce the computational cost and make it more convenient for implementation, some transformations and approximations are usually adopted for Equation (7). Equation (7) can represent many widely used imaging algorithms, such as backprojection algorithms, range Doppler algorithms, chirp scaling algorithms, and ωK algorithms [1].

**Figure 2.** Matched filtering example. Two closely spaced targets cannot be resolved.

#### *2.4. Regularization Method*

Examining the LS formula (Equation (5)), it can be seen that it only relies on the observed data. In order to make the ill-posed inverse problem become well-posed, we can add an extra constraint to the LS formula [8–10]. This leads to the regularization method, which can be denoted as

$$\hat{\mathbf{g}} = \operatorname\*{argmin}\_{\mathbf{g}} \left\{ \|\mathbf{s} \cdot \mathbf{A} \mathbf{g}\|\_{2}^{2} + \lambda L(\mathbf{g}) \right\} \tag{8}$$

where λ is the regularization parameter and *L*(**g**) is the added penalty function. In order to make the solution of Equation (8) closer to the true value, *L*(**g**) should represent appropriate prior knowledge for the problem.

A typical choice of *L*(**g**) is

$$L(\mathbf{g}) = \|\mathbf{g}\|\_p^p \tag{9}$$

where ·*p* denotes the *p*-norm, i.e.,

$$\|\mathbf{g}\|\_{p} = \begin{cases} \left(\sum\_{i=1}^{N} \left|\mathbf{g}\_{i}\right|^{p}\right)^{1/p} & p > 0\\ \text{Number of nonzero elements in } \mathbf{g} & p = 0 \end{cases} \tag{10}$$

Then, Equation (8) can be rewritten as

$$\hat{\mathbf{g}} = \operatorname\*{argmin}\_{\mathbf{g}} \left\{ \|\mathbf{s} \cdot \mathbf{A} \mathbf{g}\|\_{2}^{2} + \lambda \|\mathbf{g}\|\_{\mathcal{V}}^{\mathcal{V}} \right\} \tag{11}$$

The choice of *p* can control the result of the regularization method. If we want to enforce sparsity in the result, we should choose *p* in the range 0 ≤ *p* ≤ 1 [16,17]. For *p* = 1, Equation (11) can be compared to the Lasso solution of the CS type methods [16]. Equation (11) can be solved by gradient search algorithms, such as the Newton iteration [22].

#### *2.5. Bayesian Maximum a Posteriori Estimation*

It should be noted that in Equation (11), the added constraint term λ**g** *p p* represents prior knowledge [17,23]. Another prior knowledge-based estimation method is Bayes theory. The main idea behind the Bayesian estimation framework is to account explicitly for the errors, and also for incomplete prior knowledge. Assuming that the noise **n** in Equation (3) is white and Gaussian, we have

$$p(\mathbf{n}) \propto \exp\left\{-\frac{1}{2\sigma^2} \|\mathbf{n}\|\_2^2\right\} \tag{12}$$

where σ2 is the noise variance. Then we obtain the expression of likelihood

$$p(\mathbf{s}|\mathbf{g}) \propto \exp\left\{-\frac{1}{2\sigma^2} \|\mathbf{s} - \mathbf{g}\|\_2^2\right\} \tag{13}$$

We assume that the scene has a prior probability density function, as

$$p(\mathbf{g}) \propto \exp\{-\alpha \|\mathbf{g}\|\_{p}^{p}\} \tag{14}$$

If 0 ≤ *p* ≤ 1, the magnitude of the scene is more likely to concentrate around zero, which implies that the scene is sparse. For a review on sparsity enforcing priors for the Bayesian estimation approach, the reader can refer to [6].

Using the prior probability density of **g** shown in (14), and according to the Bayes rule, we obtain

$$p(\mathbf{g}|\mathbf{s}) = \frac{p(\mathbf{s}|\mathbf{g})p(\mathbf{g})}{p(\mathbf{s})} \propto \frac{1}{p(\mathbf{s})} \exp\left\{-\frac{1}{2\sigma^2} \left\|\mathbf{s} - \mathbf{g}\right\|\_2^2 - \alpha \left\|\mathbf{g}\right\|\_p^p\right\} \tag{15}$$

Then the MAP estimate can be obtained easily as

$$\hat{\mathbf{g}} = \arg\max\_{\mathbf{g}} p(\mathbf{g}|\mathbf{s}) = \arg\min\_{\mathbf{g}} \|\mathbf{s} - \mathbf{g}\|\_2^2 + 2\sigma^2 \alpha \|\mathbf{g}\|\_p^p \tag{16}$$

Comparing Equations (11) and (16), it can be seen that when λ = 2σ2α, these two equations are equivalent, i.e., the regularization method is equivalent to Bayesian MAP estimation.

#### *2.6. Compressed Sensing Method*

For the observation model shown in Equation (3), if the scene (i.e., **g**) is sparse, according to CS theory, it can be stably reconstructed using reduced data samples. The reconstruction method can be written as [16,17]

$$\hat{\mathbf{g}} = \operatorname\*{argmin}\_{\mathbf{g}} \|\mathbf{g}\|\_{0} \text{ s.t. } \|\mathbf{s} \cdot \mathbf{A} \mathbf{g}\|\_{2}^{2} < \varepsilon \tag{17}$$

where s.t. means subject to and ε denotes the allowed data error in the reconstruction process.

Equation (17) is NP-hard and computationally difficult to solve [17]. Matching pursuit is an approximate method for obtaining an 0 sparse solution. In CS theory, a more tractable approach is taking the 1-norm instead of the 0-norm, which is called the 1 relaxation:

$$\hat{\mathbf{g}} = \operatorname\*{argmin}\_{\mathbf{g}} \|\mathbf{g}\|\_1 \text{ s.t. } \|\mathbf{s} \cdot \mathbf{A} \mathbf{g}\|\_2^2 < \varepsilon \tag{18}$$

If **g** is sparse and **A** satisfies some specific conditions, Equations (18) and (17) will have the same solution, and this solution is the exact or approximate recovery of **g** [16,17]. Equation (18) can be solved using convex programming, which is more tractable than the original 0-norm minimum problem. Unlike the matched filtering method, CS method does not have an exact or pre-defined resolution, since it is a non-linear method. Generally, the resolution capability of the CS method is much better than the matched filtering method if the targets are sparse.

Figure 3 shows an example of compressed sensing. The simulated scene is the same as the matched filtering example shown in Figure 2. Only 1/20 signal samples are used for the CS reconstruction. It can be seen that the two closely spaced targets are well resolved. This implies that the CS method can obtain better results using less data than the matched filtering method. The reason is that prior information concerning signal sparsity is utilized in the CS model.

Equation (18) is a constrained optimization problem. According to the Lagrange theory, it can be transformed into an unconstrained optimization problem, which will have the same form as Equation (11). For appropriate choices of λ and *p* = 1, Equations (11) and (18) will be equivalent [16,17]. This implies that CS is a special case of the regularization method.

**Figure 3.** Compressed sensing example; closely spaced targets are well resolved.

#### *2.7. Summary of Radar Imaging Methods*

The above subsections introduced the LS estimator, matched filtering, regularization methods, Bayesian MAP estimation, and the CS method. In this subsection, we will summarize these methods and analyze their connections.

Table 1 lists the main characteristics and describes some connections between these imaging methods. The LS estimation only relies on the observed data, and cannot solve the ill-posed radar imaging problem efficiently. The matched filtering method can be viewed as using an approximation to avoid the ill-posed term in the LS solution. The regularization method, Bayesian MAP estimation, and the CS method exploit prior knowledge concerning the targets in addition to the observed data, and they are equivalent in some cases.

Table 1 also shows the equivalent geometric illustration for each method in R2. The observation equation can only confine the solution to a hyperplane (which becomes a line in R2), but cannot reliably produce a certain solution [17,23]. The other methods aim at obtaining a stable solution close to the true value, using some modifications that represent prior knowledge concerning the targets.

Figure 4 shows the block diagram and the relationship of the radar imaging methods. All of the radar imaging methods can be divided into two branches. The first branch does not use the prior information of the targets or scene, and it leads to the linear imaging methods; the most typical and widely used one in this branch is matched filtering. Another branch uses the prior information of the targets or scene. This leads to the non-linear methods. The most recently developed methods, including regularization methods, Bayesian methods, and CS methods belong to this branch.

**Table 1.** Characteristics and connections of radar imaging methods.

**Figure 4.** Block diagram and relationship of the radar imaging methods.

#### **3. Challenges and Advances in Compressed Sensing-Based Radar Imaging**

The use of regularization methods in radar imaging goes back at least to the year 2000 [21,24]. Since the CS theory was proposed in 2006, it has been explored for a wide range of radar [25–33] and radar imaging applications [4,34–38], including synthetic aperture radar (SAR) [39–42], inverse SAR (ISAR) [43–45], tomographic SAR [46–51], three-dimensional (3D) SAR [52–54], SAR ground moving target indication (SAR/GMTI) [55–61], ground penetrating radar (GPR) [62–64], and through-the-wall radar (TWR) [65–67]. In this paper, we will focus on two-dimensional (2D) imaging radar systems, i.e., SAR, GPR, and TWR.

After several years of development, although many interesting ideas have been presented in this area, there still exist a number of challenges, both in theory and practice [68]. The state of the art in this area has not yet reached the stage of practical application. We will present some challenges as well as recent advances in this part of the paper.

#### *3.1. Sampling Scheme*

CS usually involves random under-sampling [16,17]. A widely used waveform in traditional radar imaging is the linear frequency modulated (LFM) waveform. If we adopt the LFM waveform in CS-based radar imaging, a random sampling analog to digital (A/D) converter is needed, which is not easily realized in practice. This will require extra hardware components, which means that LFM waveforms are not ideally suited for CS.

Recently, many researchers have found that the stepped frequency waveform is much more suitable for CS than the LFM waveform [35,62,63,66,69]. Sparse and discrete frequencies are more convenient for hardware implementation. For a CS-based radar imaging system, a stepped frequency waveform may be the preferred choice. In practical application, a set of adjustable pseudorandom numbers can be generated to select the frequency points in the stepped frequencies. In this way, randomly generated frequencies, i.e., random and sparse measurement, can be realized, and the CS-based imaging model can be implemented.

Figures 5 and 6 show an example for CS-based stepped frequency radar imaging. The main equipment in the experimental system is a vector network analyzer (VNA). The experiment is carried out in a non-reflective microwave chamber. Five targets in the scene are shown in Figure 5. Figure 6a shows the backprojection result, using the fully sampled data (81 azimuth measurements × 2001 frequencies). Figure 6b shows the CS reconstruction result using under-sampled data (27 azimuth measurements × 128 frequencies). Considering the aspects of resolution and sidelobe levels, the CS reconstruction result is even better than the backprojection result, although it uses less sampled data. The reason is that prior information concerning signal sparsity is used in the CS model, while the backprojection method uses no prior information.

**Figure 5.** Experimental scene for CS-based stepped frequency radar imaging. (**a**) Five reflectors in the microwave chamber. (**b**) Transmitter and receiver antennas.

**Figure 6.** (**a**) Backprojection result of full data (81 azimuth measurements × 2001 frequencies). (**b**) CS result of under-sampled data (27 azimuth measurements × 128 frequencies).

#### *3.2. Computational Complexity*

In the regularization or CS model for a 2D radar imaging system, the 2D observed data and the 2D scene grid are both stacked into column vectors. This will lead to a huge size measurement matrix. For example, the original fully sampled data are 2048 × 2540 points (azimuth × range); if a 512 × 512 pixel image is reconstructed from a reduced sampling data consist of 256 × 256 points. Then the size of the matrix **A** is 65,536 × 262,144. Since regularization or CS reconstruction is a non-linear process, such a large measurement matrix will result in a huge computational burden for image reconstruction. In addition, the total memory to access the measurement matrix is 128 gigabytes (assuming float point and complex numbers are used). This is a too much memory space for normal desktop computers. Considering that data size is usually larger than the above example in practice, it is difficult for conventional methods to reconstruct a moderate-size scene by using normal computers.

A common idea for reducing computational complexity and memory occupancy is to split big data into sets of small data [70]. Based on this thought, a segmented reconstruction method for CS based SAR imaging has been proposed [71]. In this method, the whole scene is split into a set of small subscenes. Since the computational complexity is non-linear to the data size, the reconstruction time can be reduced significantly. The sensing matrices for the method proposed in [71] are much smaller than those for the conventional method. Therefore, the method also needs much less memory. Due to the short reconstruction time and lower memory requirement of the method proposed in [71],

reconstructing a moderate-size scene in a short time is no longer a difficult task. The processing steps of the segmented reconstruction method are shown in Figure 7.

**Figure 7.** Processing steps of the segmented reconstruction method for CS-based synthetic aperture radar (SAR) imaging (taken from [71]).

Figures 8 and 9 show an example of the segmented reconstruction method [71]. Figure 8 shows the experimental scene of an airborne SAR system, which contains six trihedral reflectors. Figure 9a shows the conventional CS reconstruction result, where the reconstruction time is 44,032 s (12 h 14 min). The whole scene is split into five segments, and Figure 9b shows the segmented reconstruction result, where the reconstruction time is now reduced to 1498 s (25 min). It can be seen that, using the segmented reconstruction method, the reconstruction time is significantly reduced, while the reconstruction precision is nearly the same.

(**a**) Sight A of the scene. (**b**) Sight B of the scene.

**Figure 8.** Trihedral reflectors in the scene. Trihedral reflectors 1–4 are large, and trihedral reflectors 5 and 6 are small (taken from [71]).

**Figure 9.** (**a**) Conventional CS reconstruction result (reconstruction time = 44,032 s). (**b**) Segmented reconstruction result (reconstruction time = 1498 s) (taken from [71]).

#### *3.3. Sparsity and Sparse Representation*

Sparsity of the scene is an essential requirement for sparse regularization or CS methods. For an SAR scene, an extended scene is usually not sparse in itself (not sparse in the canonical basis), except for the case of a few dominant scatterers in a low reflective background [35]. Therefore, a sparse representation is needed to use a sparsity-driven method.

CS-based optical imaging has successfully used sparse representations [72]. However, radar imaging involves complex-valued quantities; the raw data and the imaging result are both complex-valued. Since the phase of the scene are potentially random, it is very difficult to find a transform basis to sparsify a complex-valued and extended scene [73,74].

Structured dictionaries and dictionary learning ideas are proposed in [75] and [76], respectively. An alternative approach is to handle the magnitude and phase separately [41]. Although the phase of the scene is potentially random, the magnitude of the scene usually has better sparse characteristics. However, this approach has a much higher computational complexity than standard CS reconstruction. Another method investigates physical scattering behavior [4,77]. For example, a car can be represented as the superposition of responses from plate and dihedral shapes.

Figure 10 shows a simulation example for an extended and complex-valued scene. There are two extended objects in the scene, one of which has a round shape while the other has a rectangular shape. Both the two objects have random phases associated with them. It can be seen that the DCT (Discrete Cosine Transform) results of the magnitude are sparse.

Figure 11a shows the result of matched filtering. Since the random phase leads to speckle, it can be seen that although the scene has a smooth shape, the matched filtering result has obvious fluctuation. Figure 11b shows the result of conventional CS reconstruction without sparse representation. The reconstruction algorithm is SPGL1 [78]. Since the scene is not sparse in the canonical basis, the reconstruction is not accurate. Figure 11c shows the result of the method using a magnitude sparse representation [41]; it can be seen that the reconstruction result is much better than Figure 11a,b. Figure 11d shows the result of the method using the improved magnitude sparse representation method proposed in [79]. In the proposed method, besides the sparsity, the real-valued information of the magnitude and the coefficient distribution of the sparse representation are also utilized. It can be seen that both the shape and speckle are further improved.

**Figure 10.** (**a**) Magnitude of the scene, (**b**) phase of the scene, (**c**) DCT result of the magnitude (taken from [79]).

**Figure 11.** Simulation results: (**a**) matched filtering result, (**b**) conventional CS reconstruction result without sparse representation, (**c**) result of the method with magnitude sparse representation, and (**d**) result of the method with improved magnitude sparse representation (taken from [79]).

Figure 12 shows the real data results. The raw data is acquired by an airborne SAR system. Figure 12 contains a scene of farmland with trellises. The reflectivity from the trellises is very strong. From the real data result, it can be seen that CS with the improved magnitude sparse representation method can produce an image with less speckle and clearer edges of different regions than the previous methods.

**Figure 12.** Real data reconstruction results (scene of farmland with trellises): (**a**) matched filtering result (full data), (**b**) conventional CS reconstruction result without sparse representation, (**c**) result of CS with magnitude sparse representation, and (**d**) result of CS with improved magnitude sparse representation (taken from [79]).

#### *3.4. Influence of Clutter*

Another practical case is when the targets of interest are sparse, but there also exists clutter in the scene. Clutter arises from reflections within the scene, so the image may no longer be sparse if significant clutter returns are present. Typical examples include GPR and TWR imaging. The interesting targets, such as landmines and humans, are usually sparse, but they are often buried in the ground surface clutter and wall clutter.

Some methods have been proposed to remove the ground surface clutter and wall clutter for downward-looking GPR and TWR [64,65]. These methods are effective in cases when the clutter is concentrated in a fixed range cell or limited to several range cells.

Another scenario is TWR/SAR imaging of moving targets. A sparsity-driven change detection method is proposed in [67]. The stationary targets and clutter are removed via change detection, and then CS reconstruction is applied to the resulting sparse scene. In [55], a SAR/GMTI method using distributed CS is proposed, which can cope with the non-sparse stationary clutter.

A more difficult case is when both the targets and clutter are stationary, and the clutter is distributed over the whole scene. Forward-looking GPR may fall into this category. Figure 13 shows a real data example for this case. In such a scenario, shrubs and rocks above the ground surface may cause strong azimuth clutter. Short range clutter is usually also strong, due to the large grazing angle and short range. Besides the strong clutter far away from the target (landmine), there is also ground surface clutter around the target. In [68], an idea is proposed to build a model in which the clutter is also taken into account as a norm in the objective function. In [80], the forward-looking clutter is suppressed in two steps. In the first step, the strong clutter outside of the reconstruction region is suppressed first. In the second step, the clutter in the reconstruction region is suppressed by selecting a proper β, which represents the ratio of the non-zeros area in the reconstructed scene. The reconstruction results are shown in Figure 14.

**Figure 13.** Real data example for the clutter problem in forward-looking GPR (backprojection result using full-sampled data). Taken from [80].

**Figure 14.** Reconstruction results in clutter environment with different parameters (taken from [80]).

#### *3.5. Model Error Compensation*

In the regularization or CS methods, we usually assume that the model is exact. However, in practice, the model may also contain errors. For example, imperfect knowledge of the observation position will lead to errors in the measurement matrix. This effect resembles motion errors that arise in traditional airborne SAR imaging. Figure 15 shows the geometry of the observation position errors or motion errors in SAR.

Several methods have been proposed to deal with model errors in CS-based or sparsity-driven radar imaging. A phase error correction method for sparsity-driven SAR imaging is proposed in [81]. An autofocus method for compressively sampled SAR is proposed in [82]. This method can correct phase errors in the reconstruction process. Both the methods proposed in [81,82] deal with phase errors in the observed data, or approximately treat the observation position-induced model errors as phase errors in the observed data. In [83], the platform position errors are investigated and compensated. That method considers the azimuth offset errors and also uses some approximations.

**Figure 15.** Geometry of the observation position errors in SAR. (Taken from [84]).

In [84], a model error compensation method is proposed. An iterative algorithm cycles through steps of target reconstruction, and observation position error estimation and compensation are used. This method can estimate the observation position error exactly, while only relying on the observed data.

Figure 16 shows a real data result using the method proposed in [84]. The data set used in this figure is the same as that used for Figure 9. In the data acquisition process, the airplane is expected to fly along a straight line. However, due to the air current's influence, the trajectory of the airplane may slightly deviate from the expected one. As a result, the observation position data inevitably contain some errors.

**Figure 16.** Observation position error compensation for airborne SAR data. (**a**) Result without observation position error compensation. (**b**) Result with observation position error compensation (taken from [84]).

Figure 16a shows the original CS reconstruction result. Since the observation position errors are not compensated, it can be seen that the targets are somewhat defocused. Figure 16b shows the corresponding CS reconstruction result with compensation for observation position error. It can be seen that the focusing quality is improved using the method proposed in [84]. The peak of the targets has an increase of about 20%, and the sidelobes are also significantly reduced.

#### **4. Conclusions**

In radar imaging area, there are many relevant techniques and methods, such as matched filtering, the range Doppler algorithm, the chirp scaling algorithm, the ωK algorithm, regularized methods, and CS methods. These techniques and methods are quite different in their forms. This paper tries to understand these techniques and methods in a unified mathematical framework.

Based on theoretical analysis, it can be seen that sparsity-driven regularization or CS-based radar imaging methods have potentially significant advantages. However, although many interesting ideas have been presented, very few of them have been verified with real data. There are still many unsolved or open problems in this area. In the issues discussed in this paper, the sampling scheme, fast reconstruction strategy, and model error problems are basically solved. However, issues concerning the sparsity or sparse representation of a complex and extended scene are still not completely solved. Strong clutter may break the sparsity of a scene, while sparse representation methods for an extended scene are currently not perfect. The state of the art in these areas has not yet reached the stage of practical application, and further investigations are needed in the future.

**Funding:** This work was supported in part by the National Natural Science Foundation of China under Grant 61401474, and in part by the Hunan Provincial Natural Science Foundation under Grant 2016JJ3025.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Compressive Sensing-Based Bandwidth Stitching for Multichannel Microwave Radars**

#### **Paul Berry 1,\*, Ngoc Hung Nguyen 2 and Hai-Tan Tran 1**

Received: 18 November 2019; Accepted: 19 January 2020; Published: 24 January 2020

**Abstract:** The problem of obtaining high range resolution (HRR) profiles for non-cooperative target recognition by coherently combining data from narrowband radars was investigated using sparse reconstruction techniques. If the radars concerned operate within different frequency bands, then this process increases the overall effective bandwidth and consequently enhances resolution. The case of unknown range offsets occurring between the radars' range profiles due to incorrect temporal and spatial synchronisation between the radars was considered, and the use of both pruned orthogonal matching pursuit and refined *l*1-norm regularisation solvers was explored to estimate the offsets between the radars' channels so as to attain the necessary coherence for combining their data. The proposed techniques were demonstrated and compared using simulated radar data.

**Keywords:** radar signal processing techniques; radar imaging; multiband processing; compressive sensing; sparse reconstruction; bandwidth stitching

#### **1. Introduction**

The construction of high range resolution profiles (HRRP) of targets is a precursor to feature extraction for automatic target recognition (ATR), and normally requires the employment of a high-bandwidth waveform following detection by a lower resolution radar mode. Examples of recent papers in the non-cooperative target recognition (NCTR) literature focusing on feature extraction for ATR following HRRP construction are [1–3]. This paper considers the problem of HRRP construction, but using low resolution radars operating in different frequency bands for the purpose of combining their signals to achieve a higher resolution, and examines the problem of their data not being mutually coherent.

The ability to acquire high resolution range profiles of targets has improved over time as hardware capability has developed, with higher resolution being achieved by increasing the time-bandwidth product. In the early approaches, for narrowband radars with very limited instantaneous bandwidth, stepped-frequency waveforms were used with a single I,Q (that is, baseband quadrature) signal sample received after each frequency step. The set of samples is effectively used for the Fourier transform of the slant range profile, enabling the range response to be obtained simply by implementing an inverse Fourier transform (see, e.g., [4]). The greater the frequency range, the higher the range resolution, but the downside is that the burst of pulses can be so long that a scatterer may migrate between range cells, causing smearing of the range profile, and therefore requiring range compensation. An example of a recent paper involving the use of stepped-frequency waveforms is [5], and recent papers which have investigated the effects of target motion and aspect sensitivity are [6,7].

An advance on the stepped-frequency approach is to increase the time-bandwidth product using stretch processing, whereby a wideband LFM waveform is transmitted, and pulse-compression is achieved in hardware by mixing the received signal with an extended replica of the transmitted waveform. A point target at a particular range will manifest itself as a single frequency which is proportional to its range. The result, which is digitally sampled in time in order to facilitate the identification of the frequency components, is therefore a superposition of discrete frequencies, each corresponding to a point target at a different range. An application of the inverse Fourier transform again recovers the range profile (see, e.g., [4,8,9]).

A spectral analysis technique to improve range resolution was proposed in [10,11] based on autoregressive linear prediction. The main idea is to combine the mutually coherent signals received from multiple waveforms transmitted sequentially or concurrently, which have widely separate carrier frequencies or may even occupy entirely different frequency bands. Viewed in the spectral domain, the received signals from individual waveforms may be seen to occupy discrete wavebands which are separate or contiguous. If contiguous, then they can potentially be coherently combined to synthesise the signals that would have been received from a single wider bandwidth waveform in the manner presented in [12]. If separate, then presumably this coherent combination of signals would still be feasible, as would be the interpolation of the frequency response in the gaps between the bands under the a priori assumption that the signals are returned from discrete scatterers using spectral estimation techniques (see, e.g., [13]). Alternatively, the signals from different frequency bands can be jointly processed without explicitly filling the gaps between the bands. Since no new synthetic frequency band is actually constructed, this approach can be referred to as bandwidth stitching to distinguish it from bandwidth interpolation and extrapolation. The main challenge of this approach lies in the presence of phase errors in different frequency bands resulting from post-processed motion compensation which is often carried out separately for each frequency band. In this paper, we focus on the problem of bandwidth stitching for radar high-resolution range profiling and explore the use of sparse reconstruction to deal with the phase error problem.

It is convenient to formulate these ideas in the spectral domain, within which point scatterers appear as discrete sinusoids and which are amenable to analysis by spectral estimation techniques such as autoregression, as presented in [10,11]. Compressive sensing and sparse reconstruction, however, provide for the possibility of alternative signal representations, potentially allowing for greater flexibility and discriminating between signals of physical origin and receiver noise [13–20]. Instances of non-sinusoidal signals are waveforms in fast-time and signals returned from rotating objects when the angle of rotation is large. Compressive sensing and sparse reconstruction can also handle the situation corresponding to data being non-uniformly sampled in time or space, such as non-uniform PRF (pulse repetition frequency) waveforms and random sparse arrays.

Compressive sensing and sparse reconstruction were exploited in [17,19,21,22] to address the problem of gaps in the data both in slow-time and in frequency for inverse synthetic aperture radar (ISAR) imaging. However, these works assumed that the data were coherent across different sub-bands and that there were no model uncertainties. The work [23] took account of the possible lack of mutual coherence between the radars operating on the different sub-bands arising from incorrect timing synchronisation, or, equivalently, errors in antenna phase's centre-relative locations. This is achieved by fitting an ultra-wideband all-pole signal model to the mutually-coherent sub-bands, which is then used for bandwidth interpolation and extrapolation prior to recovering the range profile by means of an inverse Fourier transform. This paper, however, proposes the use of compressive sensing and sparse reconstruction to deal with the non-coherence problem between different sub-bands.

To address the sub-band non-coherence problem, two different approaches were explored: (i) greedy pursuit and (ii) *l*1-norm regularisation. In the first approach, pruned orthogonal matching pursuit (POMP) [24], which was originally developed for micro-Doppler parameter estimation, is adopted to deal with the dictionary mismatch which is due to the phase errors in each sub-band resulting from the motion-compensation post processing. The main idea is to parameterise the dictionary as a function of the phase errors and to construct multiple realisations of the dictionary. A selective learning process is then used to discard the dictionaries which correspond to incorrect

values of phase errors. Since a straight application of the POMP algorithm to the problem under consideration would have been computationally expensive, we first applied the POMP algorithm pairwise to sub-bands in order to estimate the phase errors, and then utilised the conventional OMP algorithm to determine the range profile based on the estimated phase error values. In the second approach, an *l*1-norm regularisation problem can be solved jointly both for the range profile vector and the phase errors [25]. The work [25] offers two solutions for joint synthetic aperture radar (SAR) imaging and phase error correction. The first solution is not applicable to the problem under consideration because a constant phase error in each sub-band is assumed. On the other hand, the second solution considers general arbitrary phase errors, and thus can be applied to our problem. Since the phase errors within each sub-band are a linear function of the range error coming from post-processing motion compensation, we also present refined variants of the second solution of [25] to take into account this underlying structure of the phase error.

The paper is organised as follows. Section 2 formulates the problem of bandwidth stitching for HRRP in the presence of phase errors. The POMP algorithm is applied in Section 3 to the bandwidth stitching problem under consideration. Section 4 presents *l*1-norm regularisation solvers. Numerical performance comparisons are provided in Section 5 and conclusions are drawn in Section 6.

#### **2. Problem Formulation**

Consider a multistatic radar system consisting of *M* radar channels on different and distinct frequency sub-bands, approximately co-located, and illuminating a common target, such that their radar lines of sight (LoS) coincide but their range profiles are out of alignment. Each channel can individually produce a one dimensional range profile of the target, but with relatively coarse resolution. The bandwidth stitching problem can be briefly stated as follows: for the *M* generally non-coherent channels, the aim is to coherently combine, or "stitch," the channels together so that they can effectively produce a single range profile with resolution corresponding to the combined overall signal bandwidth.

Let *fm*,*n* (*n* = 1, ... , *Nm*) denote the *n*th frequency bin of the *m*th channel (*m* = 1, ... , *M*). Here, *Nm* is the number of frequency bins in the *m*th channel. We base the formulation on the point-scatterer model and assume that the target can be defined as consisting of *K* scattering centres at local line-of-sight coordinates *xk* (or local "down ranges") and having complex-valued reflectivity coefficients *αk*, which are also assumed frequency-independent. The down-converted, pulse-compressed, motion-compensated signals received in each channel, in the frequency domain, can be written as

$$\mathbf{S}\_{\rm m} = [\dots, \mathbf{S}\_{m, \nu}, \dots]\_{n=1, \dots, N\_{\rm m}}^{T} \tag{1}$$

where superscript *T* denotes the transpose operation, and

$$S\_{m,n} = |A(f\_{m,n})|^2 \exp\left\{-\frac{4\pi jf\_{m,n}}{c} \Delta R\_m \right\} \sum\_{k=1}^{K} a\_k \exp\left\{-\frac{4\pi jf\_{m,n}}{c} x\_k \right\}.\tag{2}$$

Here, *A*(*fn*,*m*) represents the transmit radar waveform, the squared amplitude resulting from pulse compression processing; constant *c* denotes the speed of light; and Δ*Rm* accounts for the range errors in the motion-compensation processing. Bandwidth stitching in this context amounts to estimating these phase errors as accurately as possible.

The signals *Sm* in (1) can be rewritten in a more compact form as

$$\mathbf{S}\_{\mathfrak{m}} = \Lambda\_{\mathfrak{m}} \mathbf{F}\_{\mathfrak{m}}^{\dagger} \mathfrak{a}^{\dagger},\tag{3}$$

where

$$\Lambda\_{\mathcal{W}} = \begin{array}{c} \text{diag}\left\{ \dots, \exp\left\{ -\frac{4\pi j f\_{m,n}}{c} \Delta R\_m \right\}, \dots \right\}\_{n=1,\dots,N\_m} \end{array} \tag{4}$$

$$\begin{array}{rcl} F\_m^\dagger &=& [\dots, F\_{m,k}, \dots]\_{k=1,\dots,K} \end{array} \tag{5}$$

$$F\_{m,k}^{\dagger} = \left[ \dots, \exp\left\{ -\frac{4\pi j f\_{m,n}}{c} x\_k \right\}, \dots \right]\_{N=1,\dots,N\_m}^{T} \tag{6}$$

$$\mathfrak{a}^{\dagger} \quad = \begin{array}{c} \mathfrak{l} \\ [\ldots, \mathfrak{a}\_{k'} \ldots \mathfrak{l}]\_{k=1,\ldots,K}^{T} \end{array} \tag{7}$$

Here, "diag" denotes a diagonal matrix; **Λ***m* is referred to as the phase error matrix, of dimension *Nm* × *Nm*; *F*† *m* and *α*† are respectively, dimensions *Nm* × *K* and *K* ×1; *Sm* is a column vector of dimension *Nm* × 1; and the dagger symbol † refers to the *K* actual scatterers on the target.

To apply the sparse representation techniques of compressive sensing, we discretise the target's local range coordinate *x* using a regularly-spaced range grid {*xl*} for *l* = 1, ... , *Lx*, with *Lx K*, and construct the *Nm* × *Lx* dictionary matrices

$$F\_{\rm ll} = [\dots, F\_{\rm m,l}, \dots]\_{l=1,\dots,l\_{\rm x}} \tag{8}$$

where

$$F\_{\mathbf{m},l} = \left[ \dots, \exp\left\{ -\frac{4\pi j f\_{\mathbf{m},n}}{c} \mathbf{x}\_l \right\}, \dots \right]\_{n=1,\dots,N\_m}^T \tag{9}$$

Are the "atoms" of dictionary *Fm* in the frequency domain. The corresponding range profile vector

$$\mathfrak{a} = [\ldots, \mathfrak{a}\_l, \ldots]\_{l=1, \ldots, L\_x}^T \tag{10}$$

Spans over the range grid {*xl*}. The received signal *Sm* can also be written as

$$S\_m = \Lambda\_m F\_m \mathfrak{a}.\tag{11}$$

Since the target usually contains only a small number of dominant scattering centres relative to the total number of range resolution cells, the range profile *α* can be considered sparse (i.e., containing a small number of non-zero elements).

In the presence of unknown noise, (11) becomes

$$
\mathcal{S}\_m = \Lambda\_m F\_m \mathfrak{a} + \mathfrak{n}\_m. \tag{12}
$$

where *nm* is the additive noise for channel *m*. Stacking up the individual channel signals *S*˜ *m*, *m* = 1, . . . , *M*, gives

$$
\vec{S} = \Lambda F \mathfrak{a} + \mathfrak{n},
\tag{13}
$$

where

$$\bar{\mathbf{S}}^{\top} = [\dots, \bar{\mathbf{S}}\_{m}^{T}, \dots]\_{m=1, \dots, M}^{T} \tag{14}$$

$$\mathbf{A}\_{\;\;\;\;\prime} = \text{diag}\{\ldots, \mathbf{A}\_{\text{W}\_{\prime}}, \ldots\}\_{\:\!\;\;\text{m} = 1, \ldots, \text{M}}\tag{15}$$

$$F\_{\ldots} = [\ldots, F\_{m'}^{T}, \ldots]\_{m=1, \ldots, M}^{T} \tag{16}$$

$$\mathfrak{n}\_{\cdot} = [\ldots, \mathfrak{n}\_{m'}^T, \ldots]\_{m=1, \ldots, M}^T. \tag{17}$$

Note that *S*˜ and *n* are column vectors of size (∑*m Nm*) × 1; diagonal phase error matrix **Λ** is of size (∑*m Nm*) × (∑*m Nm*); dictionary matrix *F* is (∑*m Nm*) × *Lx*; and the range profile *α* is again a column vector of size *Lx* × 1. Stacking the received signals amounts to a vertical stacking of the dictionary matrices from all channels and a diagonal concatenation of the corresponding phase error matrices. The stacking of multiple channels in this manner can improve the estimation accuracy for *α*, as will be demonstrated later in the paper.

A *problem statement* can thus be expressed as follows: given *S*˜ as the measured signal,

$$\text{find } \mathfrak{a} \text{ and } \Lambda, \text{ subject to } \left\{ \begin{array}{l} \mathcal{S} \approx \Lambda F \mathfrak{a}, \\ \mathfrak{a} \text{ is sparse} \end{array} \right. \tag{18}$$

The estimation of *α* over {*xl*} is the process of range profiling, giving the main desired output, whereas the estimation of the phase error matrix **Λ** is really only a necessary intermediate result; it is a function of Δ*R*1, Δ*R*2, ... , and Δ*RM* (recall that Δ*Rm* is the range estimation error resulting from the motion-compensation process for channel *m*). Furthermore, since these errors arise from a lack of precise knowledge of the relative locations of the radar channels' phase centres and are small relative to a range resolution cell, we may assume, without loss of generality, that Δ*R*1 = 0.

#### **3. Greedy Pursuit Solutions**

In this section, we adopt the pruned OMP (POMP) technique, which was originally proposed for micro-Doppler parameter estimation, [24], for the problem of bandwidth stitching for range profiling. We start with the simplest case of two channels and then generalise it to the multiple channel case.

#### *3.1. The Two-Channel Case*

For this case, the signal model in (13) can be expressed as

$$\mathcal{S} = \Lambda(\Delta \mathcal{R}\_2) F \mathfrak{a} + \mathfrak{n}, \tag{19}$$

where **Λ**(Δ*R*2) is a function of the single unknown relative range error Δ*R*2,

$$\mathbf{A}(\Delta \mathbb{R}\_2) = \text{diag}\{I\_{\mathbb{N}\_1}, \mathbf{A}\_2(\Delta \mathbb{R}\_2)\},\tag{20}$$

with

$$\Lambda\_2(\Delta R\_2) = \text{diag}\left\{ \dots, \exp\left\{ -\frac{4\pi j f\_{2,n}}{c} \Delta R\_2 \right\}, \dots \right\}\_{n=1,\dots,N\_2} \tag{21}$$

In addition to the sparse range profile vector *α*, Δ*R*2 is the only additional unknown parameter to be estimated. Let us rewrite (19) as

$$
\vec{S} = \Phi(\Delta R\_2)\,\mathfrak{a} + \mathfrak{n} \tag{22}
$$

where

$$
\Phi(\Delta R\_2) = \mathbf{A}(\Delta R\_2) \, F.
$$

In this form, the problem can be viewed as a joint sparse reconstruction and parameter estimation problem with the parametric dictionary **Φ**(Δ*R*2) itself a function of the parameter Δ*R*2. This can be considered as a special dictionary learning problem where the objective is to solve simultaneously for both the sparse solution of *α* and the range error Δ*R*2.

To solve this problem, we adopt the POMP technique [24], which embeds a pruning operation into the iterative process of OMP. The main idea of POMP is to construct multiple realisations of the dictionary **Φ** based on a number *L*Δ of candidate values of Δ*R*2; the OMP algorithm is applied to each dictionary realisation to find the atom which correlates most strongly with the current residual for that dictionary, and to recompute the residual with that atom's contribution to the residual removed. To overcome possible excessive computations arising from outlier candidate values of Δ*R*2, a pruning operation is performed to exclude the half of the dictionaries which yield the largest residual errors, until a single dictionary realisation remains. The OMP iterations for the remaining dictionary are continued until a termination criterion is satisfied. The candidate value of Δ*R*2 corresponding to this

dictionary gives the final estimate of Δ*R*2. The basic POMP algorithm is summarised in Table 1 and its computational cost is shown in Appendix C to be of the order of *L*Δ*NmLx*.

**Table 1.** The pruned orthogonal matching pursuit (POMP) algorithm (*M* = 2).

#### INPUT:

#### PROCEDURE:

	- set the initial indexes of active dictionaries to Θ1 = {1, . . . , *L*Δ};
	- set the corresponding residual vectors to *r*1 = ··· = *rL*Δ = *S*˜;
	- set the initial support Λ to ∅, the null set;
	- **for** every *l* ∈ Θ*i*, perform OMP as follows

$$\begin{array}{c} \text{- } \mathsf{Identity} \\ \begin{array}{c} c\_{l} = \boldsymbol{\Phi}\_{l}^{H} r\_{l} \\ \boldsymbol{j}\_{l} = \text{arg}\,\mathsf{max}\_{j} \,|\,\mathsf{c}\_{j}| \\ \mathsf{- } \mathsf{Merge} \,\mathsf{supp}\,\mathsf{rots} \\ \begin{array}{c} \Lambda\_{l} = \Lambda\_{l} \cup \boldsymbol{j}\_{l} \\ \mathsf{Update}^{\*} \end{array} \\ \mathsf{- } \mathsf{Update}^{H} \colon \\ \begin{array}{c} \mathsf{\hat{\mathsf{\hat{\boldsymbol{I}}}}\_{l} \Lambda\_{l} = \left(\boldsymbol{\Phi}\_{l,\Lambda\_{l}}^{H} \boldsymbol{\Phi}\_{l,\Lambda\_{l}}\right)^{-1} \boldsymbol{\Phi}\_{\boldsymbol{I},\Lambda\_{l}}^{H} \boldsymbol{\mathcal{S}} \\ \mathsf{r}\_{l} = \boldsymbol{\mathcal{S}} - \boldsymbol{\Phi}\_{l,\Lambda\_{l}}^{\boldsymbol{\Lambda}} \boldsymbol{\hat{\mathsf{\hat{\boldsymbol{I}}}}\_{l} \boldsymbol{\Lambda}\_{l} \end{array} \end{array}$$

**if** |Θ*i*| > 1

Remove indices of |Θ*i*|/2 candidate dictionaries that correspond to |Θ*i*|/2 largest residual errors from |Θ*i*|

**end if**

**end for**

OUTPUT:

∗ **Φ***l*,Λ*l* consists of the columns of **Φ***l* with indices belonging to Λ*l* and *α*ˆ Λ*l* consists of the elements of *α*ˆ *l* with indices belonging to Λ*l*.

*3.2. The Multi-Channel Case*

For this case, the noisy signal model (19) becomes

$$\mathcal{S} = \Phi(\Delta R\_{2\prime}, \dots, \Delta R\_M)\mathfrak{a} + \mathfrak{n} \tag{23}$$

where

$$\Phi(\Delta R\_2, \dots, \Delta R\_M) = \text{diag}\left\{ I\_{\mathbb{N}\_1}, \Lambda\_2(\Delta R\_2), \dots, \Lambda\_M(\Delta R\_M) \right\} \text{ F.} \tag{24}$$

Here, the dictionary matrix is a function of the (*M* − 1) unknowns Δ*R*2, Δ*R*3, . . . and Δ*RM*.

The POMP algorithm could be extended to multiple channels by computing candidate dictionaries based on a multi-dimensional grid of candidate values for Δ*R*2, Δ*R*3, ... , and Δ*RM*. The grid would consist of a total of (*M* − 1) dimensions, where the *m*th dimension corresponds to the unknown range error Δ*Rm*+1 of the (*m* + 1) channel. Note that only a one-dimensional grid for Δ*R*2 is required for the case of two channels. However, the cardinality of the dictionary set is exponentially dependent on the number of available channels; i.e., *V*(*M*−1), where *V* denotes the number of grid points in each parameter dimension. As a result, although this extension would be simple and straightforward, it is computationally expensive.

To alleviate this computational burden, we instead apply the POMP algorithm pairwise to channels in order to estimate the range errors Δ*R*2, ... , Δ*RM* relative to the first channel, and then utilise the conventional OMP algorithm to determine the range profile vector *α* based on these estimated values of Δ*R*2,..., Δ*RM*. The procedure is summarised as below:

**STEP 1:** Estimation of range errors.

	- Calculate input signal:

$$\vec{S}\_{1m} = [\vec{S}\_1^T, \vec{S}\_m^T]^T \tag{25a}$$

$$F\_{1m} = [F\_1^T, F\_m^T]^T. \tag{25b}$$

$$\Phi\_{1m,l} = \text{diag}\{I\_{N\_1}, \Lambda\_m(\Delta R\_m^{(l)})\} F\_{1m}. \tag{25c}$$

End for.

**STEP 2:** Estimation of range profile vector.

• Compute signal and dictionary.

$$\mathcal{S} = [\dots, \mathcal{S}\_{m'}^T \dots \mathcal{I}]\_{m-1, \dots, M}^T \tag{26}$$

$$F = [\dots, \mathbf{F}\_m^T, \dots]\_{m=1, \dots, M}^T \tag{27}$$

$$\mathbf{A} = \text{diag}\{\dots, \mathbf{A}\_m(\Delta R\_m), \dots\}\_{m=1,\dots,M} \tag{28}$$

$$
\Phi = \Lambda F.
\tag{29}
$$

• Estimate *α* using OMP given *S*˜ and **Φ**.

The computational cost of the general POMP algorithm is shown in Appendix C to be of the order of (*M* − 1)*L*Δ*NmLx*.

#### **4.** *L***1-Norm Regularisation Approach**

The sparse reconstruction problem (18) can be solved via the following *l*1 regularised optimisation:

$$\min\_{\mathbf{a},\mathbf{A}} \left\{ \|\mathbf{S} - \mathbf{A}\mathbf{F}\mathbf{a}\|\|\_{2}^{2} + \mu \|\mathbf{a}\|\|\_{1} \right\},\tag{30}$$

where *μ* is a regularisation parameter. It should be emphasized that this is not a conventional *l*1 regularisation formulation because of the unknown phase error matrix **Λ** resulting from the estimation error of the motion-compensation phase. Therefore, **Λ** must be jointly estimated with *α*:

$$\{\hat{\mathfrak{a}}, \hat{\Lambda}\} = \operatorname\*{arg\,min}\_{\mathfrak{a}, \Lambda} \left\{ \|\mathbf{S} - \Lambda \mathbf{F} \mathfrak{a}\|\|\_{2}^{2} + \mu \|\mathfrak{a}\|\|\_{1} \right\}.\tag{31}$$

Two solutions for this joint estimation problem were presented in [25]. The first solution assumes that the phase error matrix for the *m*th sub-band is modelled as

$$\Lambda\_{\mathfrak{M}} = \exp\{j\phi\_{\mathfrak{M}}\} \ I\_{\mathcal{N}\_{\mathfrak{M}} \times \mathcal{N}\_{\mathfrak{M}}}.\tag{32}$$

In other words, the phase errors for different frequency bins of a particular sub-band are identical. However, this assumption is invalid in the problem under consideration because the phase error is a function of frequency and thus has different values for different frequency bins. Therefore, that solution is not applicable in this case. The second solution considers a general phase error matrix

$$\Lambda\_{\mathfrak{M}} = \text{diag}\{\dots, \exp\{j\phi\_{\mathfrak{m},\mathfrak{n}}\}, \dots\}\_{\mathfrak{n} = 1, \dots, N\_{\mathfrak{m}}} \tag{33}$$

where the phase errors *φm*,*n* can be arbitrary. Although this solution can be used, it does not exploit the underlying structure of the phase errors; i.e., *φm*,*n* = −4*π fm*,*n c* Δ*Rm*. In what follows, we will also present other refined versions, building on the second solution of [25], while exploiting prior knowledge of the structure of the phase error.

The *l*1 norm can be approximated as [26–29]:

$$||\mathfrak{a}||\_1 \approx \sum\_{l=1}^{L} \left( |\mathfrak{a}\_l|^2 + \delta \right)^{1/2} \tag{34}$$

In order to overcome the nondifferentiably of the *l*1 norm at the origin. Here, *δ* is a small non-negative parameter. Using this approximation, the minimisation problem in (31) becomes

$$\{\hat{\mathbf{a}},\hat{\mathbf{A}}\} = \operatorname\*{arg\,min}\_{\mathbf{a},\mathbf{A}} \left\{ \left\| \mathbf{S} - \mathbf{A}\mathbf{F}\mathbf{a} \right\|\_{2}^{2} + \lambda \sum\_{l=1}^{L} \left( \left| \alpha\_{l} \right|^{2} + \delta \right)^{1/2} \right\}.\tag{35}$$

The solution of (35) tends to the solution of (31) as *δ* approaches zero. Therefore, a small value of *δ* should be used to ensure the validity of this approximation. The quasi-Newton approach can be adopted to solve the modified *l*1 regularised optimisation (35), as below.

The gradient of the objective function of (35) is given by

$$\nabla(\mathfrak{a}) = H(\mathfrak{a})\mathfrak{a} - 2\,\,\mathcal{F}^H \boldsymbol{\Lambda}^H \mathbf{S},\tag{36}$$

where the superscript *H* denotes the Hermitian transpose operation. Here, *H* is the Hessian matrix given by

$$H(\mathfrak{a}) = 2\,\mathrm{F}^{H}\Lambda^{H}\Lambda\mathrm{F} + \lambda\mathcal{W}(\mathfrak{a}) = 2\,\mathrm{F}^{H}\mathrm{F} + \lambda\mathcal{W}(\mathfrak{a}),\tag{37}$$

where

$$\mathcal{W}(\mathfrak{a}) = \text{diag}\left\{ \dots, \left( |a\_l|^2 + \delta \right)^{-1/2}, \dots \right\}. \tag{38}$$

Since the Hessian matrix is a function of the unknown *α*, the minimisation (35) is solved iteratively. Given the estimates *α*ˆ(*i*) and **Λ**ˆ (*i*) from a previous iteration *i*, the new solutions at iteration *i* + 1 are obtained in the following two steps.

1. Calculate ˆ*α*(*i* + 1) by setting ∇(*α*) = **0** given *H*(*α*ˆ(*i*)) and **Λ**ˆ (*i*):

$$\begin{split} \hat{\mathfrak{a}}(i+1) &= 2\left(H(\hat{\mathfrak{a}}(i))\right)^{-1} \mathbf{F}^{H}(\hat{\mathbf{A}}(i))^{H} \mathbf{S} \\ &= \left(\mathbf{F}^{H}\mathbf{F} + \frac{1}{2}\lambda\mathbf{W}(\hat{\mathfrak{a}}(i))\right)^{-1} \mathbf{F}^{H}(\hat{\mathbf{A}}(i))^{H} \mathbf{S}. \end{split} \tag{39}$$

2. Calculate **Λ**ˆ (*i* + 1) given ˆ*α*(*i* + 1). The phase error matrix **Λ**ˆ (*i* + 1) is obtained by solving:

$$\hat{\mathbf{A}}(i+1) = \underset{\mathbf{A}}{\arg\min} \|\mathbf{S} - \mathbf{A}\mathbf{F}\mathbf{\hat{a}}(i+1)\|\_{2}^{2} \tag{40}$$

or equivalently

$$\hat{\mathbf{A}}\_{m}(i+1) = \underset{\mathbf{A}\_{m}}{\text{arg min}} \|\mathbf{S}\_{m} - \mathbf{A}\_{m}\mathbf{F}\_{m}\mathbf{a}(i+1)\|\_{2}^{2} \tag{41}$$

for *m* = 2, . . . , *M*. Note that **Λ**1 = *IN*1 (since Δ*R*1 = 0); thus, no estimation is required for **Λ**1.

The algorithm may be halted when the objective function falls below a threshold, or when a maximum number of iterations is reached, or when the relative change in the objective function falls below a threshold.

Various methods for calculating the phase error matrix **Λ**ˆ (*i* + 1) in Step 2 are given in the following sections.

#### *4.1. Unstructured Approach*

Ignoring the underlying structure of the phase errors, i.e., *φm*,*n* = −4*π fm*,*n c* Δ*Rm*, **Λ**ˆ *m* can be considered as a diagonal matrix with arbitrary elements *φm*,*n*:

$$\mathbf{A}\_{\mathcal{W}} = \text{diag}\{\dots, \exp\{j\phi\_{\mathcal{W},\mathcal{W}}\}, \dots\}\_{n=1,\dots,N\_{\mathcal{W}}}.\tag{42}$$

Therefore, *φm*,*n* can be estimated as [25]

$$\phi\_{m,n}(i+1) = \tan^{-1} \frac{\circledast \{ S\_{m,n} \hat{Y}\_{m,n}^\*(i+1) \}}{\Re \{ S\_{m,n} \hat{Y}\_{m,n}^\*(i+1) \}} \tag{43}$$

where {·} and {·} denote operations to extract the imaginary and real parts of a complex number, and tan−1 stands for a four-quadrant arctangent operation. Here, *Y*ˆ *m*,*n*(*i* + 1) is the *n*th element of *Y*ˆ *m*(*i* + 1) which is defined as *Y*ˆ *m*(*i* + 1) = *Fmα***ˆ**(*i* + 1). As a result, we obtain:

$$\dot{\mathbf{A}}\_{m}(i+1) = \text{diag}\{\dots, \exp\{j\hat{\phi}\_{m,n}(i+1)\}, \dots\}\_{n=1,\dots,N\_{m}}.\tag{44}$$

#### *4.2. Gauss–Newton Approach*

Taking into account the underlying structure of the phase errors, **Λ***m*(Δ*Rm*) is a function of Δ*Rm*, and the minimisation (41) can be re-expressed as

$$
\overline{\Delta}\overline{R}\_{\mathfrak{M}}(i+1) = \underset{\Delta R\_{\mathfrak{m}}}{\text{arg min}} \left\| \mathbf{S}\_{\mathfrak{m}} - \mathbf{A}\_{\mathfrak{m}}(\Delta R\_{\mathfrak{m}}) F\_{\mathfrak{m}} \mathfrak{A}(i+1) \right\|\_{2}^{2}.\tag{45}
$$

By letting *em* = *Sm* − **Λ***m*(Δ*Rm*)*Fmα***ˆ**(*i* + 1), we have

$$\mathbf{e}\_{m} = [\dots, \mathbf{e}\_{m,\nu\nu}, \dots]\_{m=1,\dots,N\_m}^T \tag{46}$$

where

$$\varepsilon\_{\rm m,n} = S\_{\rm m,n} - \mathcal{U}\_{\rm m,n} \exp\left\{ -\frac{4\pi j f\_{\rm m,n}}{c} \Delta R\_m \right\} \sum\_{l=1}^{L} \hbar\_l (i+1) \exp\left\{ -\frac{4\pi j f\_{\rm m,n}}{c} \mathbf{x}\_l \right\} \tag{47}$$

and

$$\mathcal{U}\_{m,n} = \exp\{j\phi\_{m,n}\}.\tag{48}$$

Here *α*ˆ*l*(*i* + 1) is the *l* th element of *α*ˆ(*i* + 1). As we are estimating real quantities, it is more convenient to reformulate the problem as the minimisation of a real function in order to apply the Gauss–Newton. The details of the Gauss–Newton algorithm for updating Δ 2*Rm*(*i* + 1) are given in Appendix A.

Using Δ 2*Rm*(*i* + 1), we obtain the phase error matrix as

$$\hat{\Lambda}\_{m}(i+1) = \text{diag}\left\{ \dots, \exp\left\{ -\frac{4\pi j \hat{\lambda} \hat{R}\_{m}(i+1)}{c} f\_{m,n} \right\}, \dots \right\}\_{n=1,\dots,N\_{m}}.\tag{49}$$

#### *4.3. Linear Regression-Based Approach*

By noting that *φm*,*n* = −4*π*Δ*Rm c fm*,*n*, the gradient −4*π*Δ*Rm c* can be calculated via a linear least squares estimator using *φ*ˆ*m*,*n* obtained from (43) [30]. Specifically, we have

$$\left[ -\frac{4\pi\widehat{\Delta R}\_m(i+1)}{c}, \widehat{\Phi}\_m^\dagger(i+1) \right]^T = (A\_m^T A\_m)^{-1} A\_m^T b\_m(i+1) \tag{50}$$

where

$$\mathbf{A}\_{\mathfrak{m}} = [\ldots, \mathbf{A}\_{m, \mathfrak{n}} \ldots \mathbf{I}]\_{n=1, \ldots, \aleph\_{m} \prime}^{T} \text{ with } \mathbf{A}\_{m, \mathfrak{n}} = [f\_{m, \mathfrak{n}}, \mathbf{1}] \tag{51a}$$

$$\mathbf{b}\_{\rm m}(i+1) = [\dots, \hat{\phi}\_{\rm m,n}^{\rm unwarppped}(i+1), \dots]\_{\rm n-1, \dots, N\_{\rm m}}^{T}. \tag{51b}$$

Note that *φ*ˆunwrapped *m*,*n* (*i* + 1) is the unwrapped version of *φ*ˆ*m*,*n*(*i* + 1) and *φ*ˆ† *m*(*i* + 1) is an estimate for the initial phase *φ*† *m*(*i* + 1) which results from the unwrapping process. From (50), Δ 2*Rm*(*i* + 1) is obtained and can then be used for computing **Λ**ˆ *m*(*i* + 1) as in (49).

#### *4.4. Differenced-Phase-Based Approach*

Subtracting the estimated phase errors of two successive frequency bins, we obtain:

$$-\frac{4\pi (f\_{m,n+1} - f\_{m,n})}{c} \Delta R\_{\text{ll}} = \phi\_{\text{m,n+1}} - \phi\_{\text{m,n}}.\tag{52}$$

Therefore, Δ*Rm* can be estimated as [30]

$$
\hat{\Delta R}\_m(i+1) = -\frac{c}{4\Delta f\_m(N\_m - 1)} \sum\_{n=1}^{N\_m - 1} \Delta \phi\_{m,n}(i+1) \tag{53}
$$

where

$$
\Delta\phi\_{m,n}(i+1) = \tan^{-1}\frac{\sin(\phi\_{m,n+1}(i+1) - \phi\_{m,n}(i+1))}{\cos(\hat{\phi}\_{m,n+1}(i+1) - \hat{\phi}\_{m,n}(i+1))}.\tag{54}$$

Note that the four-quadrant arctangent has been used here to handle the phase wrapping. An estimate of the phase error matrix **Λ**ˆ *m*(*i* + 1) is now obtained as in (49) using Δ 2*Rm*(*i* + 1) in (53).

#### **5. Simulation and Discussion**

Numerical simulations are presented in this section to evaluate the performance of the methods described in previous sections.

#### *5.1. Scenario 1: Two Sub-Bands*

We consider a synthetic scenario with two sub-bands at carrier frequencies of *f*1 = 6 GHz and *f*2 = 8 GHz, each having a bandwidth of *B* = 300 MHz and 64 frequency steps (i.e., *N* = *N*1 = *N*2 = 64). The range profile is discretised over a grid with a length of (*N* − 1)*c*/(2*B*) = 31.5 m and a grid step of ΔGrid = *c*/(10*B*) = 0.1 m. We consider a far-field target consisting of six point scatterers which are aligned with the grid. Figure 1 plots the true range profile of the target. We set Δ*R*2 = 2.78ΔGrid for the case of existing phase errors. The signal-to-noise ratio is set to 10 dB.

Figure 2 compares the reconstructed range profiles obtained by the conventional OMP algorithm without and with the presence of phase errors. The OMP is terminated when the signal residual reaches the noise level or after 15 iterations have been carried out. We observe that OMP successfully reconstructs the range profile of the target by correctly identifying the scatterers of the target with accurate range and coefficient estimates when no phase errors exist. However, OMP provides unsatisfactory results in the presence of phase errors, where the reconstructed image is observed

as exhibiting many spurious scatterers. Similar observations are obtained for the results obtained by the conventional *l*1-norm regularised optimisation solver (without phase error correction), as shown in Figure 3. Here, we set *δ* = 10−5 and *λ* = 0.001 max |*cl*| where *cl* is the *l* th element of *c* = *FHS*˜. The *l*1-norm regularised optimisation solver is stopped if the relative change in the *l*2-norm of the range profile vector *α* falls below 10−5 or after it reaches 500 iterations. The performance degradation of these conventional sparse reconstruction techniques is not unexpected, since they were not originally developed to cope with dictionary mismatch arising from the presence of the phase errors.

**Figure 1.** True range profile of synthetic target under consideration.

(**a**) Without the presence of phase errors.

(**b**) With the presence of phase errors.

**Figure 2.** Performance of conventional OMP.

(**b**) With the presence of phase errors.

**Figure 3.** Performance of conventional *l*1-norm regularised optimisation solver.

Figure 4 shows the reconstructed range profile obtained by the POMP. POMP constructs candidate dictionaries based on a grid of Δ*R*2 with a grid step size of ΔGrid/100. The same stopping criteria of OMP is used for POMP. We observe that POMP produces a range profile which is almost identical to the ground truth, thereby demonstrating the effectiveness of POMP in terms of dealing with the phase errors between different sub-bands.

**Figure 4.** Performance of POMP in the presence of phase errors.

Figure 5 shows the results obtained by different *l*1-norm regularised optimisation solvers with phase error correction, as presented in Section 4. The same parameters and stopping criteria of the conventional *l*1-norm regularised optimisation solver as described above are used in the simulations. Although these algorithms exhibit some improvements over the conventional *l*1-norm regularisation (i.e., without phase error correction), they provide poorer results compared to that of POMP. Specifically, the peaks of the reconstructed range profiles obtained by these algorithms only appear close to but not exactly at the true scatterer positions. In addition, the magnitudes of the peaks are much smaller than the ground truth values.

(**d**) With structured DP-based error correction

**Figure 5.** Performance of *l*1-norm regularised optimisation solver with phase error correction.

The inferior performances of these algorithms can be explained by noting that the *l*1-norm regularised optimisation in (31) is a nonconvex problem due to the phase error matrix **Λ**. Figure 6 plots the objective function of (31) as a function of Δ*R*2 assuming that *α* is perfectly known. We observe that this objective function has many local maxima and minima, confirming the non-convexity of the *l*1-norm regularised optimisation in (31). The reason for this non-convexity is explained in Appendix B. Due to this nonconvexity, the iterative solvers presented in Section 4 are prone to converge to local minima; thus, limiting the effectiveness of this approach.

**Figure 6.** Illustration of the nonconvexity of the *l*1-norm regularised optimisation problem (31). The objective function of (31) is plotted against Δ*R*2 assuming that *α* is perfectly known.

The performance of the OMP, POMP, and *l*1-norm regularised optimisation methods are now compared using the earth mover's distance (EMD) between the true and reconstructed range profiles. EMD [31] is a metric estimating the distance between two distributions or equivalently the minimal amount of work required to transform one distribution to the other. Figures 7 and 8 show the EMD performance of the OMP, POMP, and *l*1-norm regularised optimisation methods, averaged over 100 Monte Carlo runs, versus different levels of SNR (signal-to-noise ratio) and phase error, respectively. It is observed that the POMP method yields the smallest EMD values amongst all algorithms considered. Since a smaller value of EMD corresponds to a higher level of similarity between the true and reconstructed range profiles, this observation indicates that the reconstructed range profile obtained by POMP is closer to the ground-truth range profiles than those obtained from the OMP and *l*1-norm regularised optimisation methods. This verifies the performance advantage of the POMP method from a statistical point of view.

**Figure 7.** Earth mover's distance (EMD) performance of the OMP, POMP, and *l*1-norm regularised optimisation methods versus various of SNRs (Δ*R*2 = 2.78ΔGrid).

**Figure 8.** EMD performance of the OMP, POMP, and *l*1-norm regularised optimisation methods versus various levels of phase error (SNR = 10 dB).

#### *5.2. Scenario 2: Four Sub-Bands*

We now consider another scenario with four sub-bands at carrier frequencies of *f*1 = 6 GHz, *f*2 = 8 GHz, *f*3 = 10 GHz, and *f*4 = 12 GHz, each having a bandwidth of *B* = 300 MHz and 64 frequency steps (i.e., *N* = *N*1 = *N*2 = *N*3 = *N*4 = 64). The range errors are set to Δ*R*2 = 2.78ΔGrid, Δ*R*3 = 1.33ΔGrid, and Δ*R*4 = 3.69ΔGrid. Other simulation parameters and the true range profile of the target remain unchanged as in the previous simulation example.

Figure 9 compares the reconstructed range profiles obtained by the conventional OMP algorithm and the POMP algorithm presented in Section 3.2. OMP results in an unsatisfactorily reconstructed range profile with many spurious peaks, as expected, because it ignores the phase errors between different sub-bands. In contrast, the POMP is capable of reconstructing the true range profile with a high accuracy thanks to the use of dictionary learning with a pruning process. Note that, given the inferior performance of the *l*1-norm regularised optimisation approach compared to POMP, as demonstrated in the previous simulation scenario, this approach is excluded from the comparison here.

(**b**) Range profile reconstructed by POMP

**Figure 9.** Performance comparison between OMP and POMP for Simulation Scenario 2 (with four sub-bands).

#### **6. Conclusions**

This paper explores the use of the POMP algorithm and *l*1-norm regularisation solvers for the problem of sparsity-driven HRRP with bandwidth stitching in the presence of phase errors. We observe that the *l*1-norm regularisation solvers do not provide significant performance improvement over the conventional sparse reconstruction algorithms due to the nonconvexity of *l*1-norm regularised optimisation when phase errors exists. In contrast, POMP is observed to be capable of effectively dealing with the phase errors and thus be able to reconstruct the range profile of the target with high accuracy. Simulation results show a significant performance improvement by POMP over OMP and the conventional and refined *l*1-norm regularisation. In future work, we propose using experimental data for a more general scenario where the true scatterers constituting the target are located in off-grid positions with respect to the dictionary grid, and the true range errors have off-grid values. We shall also consider the more general case of frequency-dependence of scatterer RCS (radar cross-section). A potential approach for this is to exploit the framework of spectral compressive sensing [32,33].

**Author Contributions:** Conceptualisation, P.B., N.H.N., and H.-T.T.; methodology, P.B., N.H.N., and H.-T.T.; simulation, N.H.N.; validation, P.B., N.H.N., and H.-T.T.; writing—original draft preparation, N.H.N.; writing—review and editing, P.B., N.H.N., and H.-T.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Defence Science and Technology Group, Australia.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Appendix A. Derivation of Gauss–Newton Algorithm for Estimation of Δ** 1*Rm***(***i* **+ 1)**

We rewrite *em*,*n* as

$$\mathcal{L}\_{m,n} = S\_{m,n} - \exp\left\{-\frac{4\pi j f\_{m,n}}{c} \Delta R\_m \right\} Z\_{m,n} \tag{A1}$$

with

$$Z\_{m,n} = \mathcal{U}\_{m,n} \sum\_{l=1}^{L} \mathfrak{A}\_{l}(i+1) \exp\left\{-\frac{4\pi j f\_{m,n}}{c} \mathbf{x}\_{l}\right\}.\tag{A2}$$

By noting that

$$\begin{split} &\mathbb{E}\exp\left\{-\frac{4\pi f f\_{m,n}}{c}\Delta R\_{m}\right\}Z\_{m,n} \\ &= Z\_{m,n}^{R}\cos\left\{\frac{4\pi f\_{m,n}}{c}\Delta R\_{m}\right\} + Z\_{m,n}^{I}\sin\left\{\frac{4\pi f\_{m,n}}{c}\Delta R\_{m}\right\} \\ &\quad + j\left(Z\_{m,n}^{I}\cos\left\{\frac{4\pi f\_{m,n}}{c}\Delta R\_{m}\right\} - Z\_{m,n}^{R}\sin\left\{\frac{4\pi f\_{m,n}}{c}\Delta R\_{m}\right\}\right) \end{split} \tag{A3}$$

where *ZR m*,*n* and *ZI m*,*n* are the real and imaginary components of *Zm*,*n*, we can decouple and stack the real and imaginary components of *em* to form a real-valued vector as

$$\boldsymbol{\mathfrak{e}}\_{\rm{m}} = \left[ \left( \boldsymbol{\mathfrak{e}}\_{\rm{m}}^{\rm{R}} \right)^{T}, \left( \boldsymbol{\mathfrak{e}}\_{\rm{m}}^{\rm{I}} \right)^{T} \right]^{T} \tag{A4}$$

where

$$\boldsymbol{\sigma}\_{m}^{R} = \left[ \dots , \boldsymbol{\varepsilon}\_{m,n}^{R}, \dots \right]\_{n=1,\dots,N\_{m}}^{T} \tag{A5}$$

$$\mathbf{e}\_{\mathfrak{m}}^{\mathrm{I}} = \left[ \dots , \mathbf{e}\_{\mathfrak{m}, \mathfrak{m}}^{\mathrm{I}} \dots \right]\_{\mathfrak{m} = 1, \dots, N\_{\mathfrak{m}}}^{\mathrm{T}} \tag{A6}$$

and

$$\varepsilon\_{m,n}^R = Z\_{m,n}^R \cos\left\{\frac{4\pi f\_{m,n}}{c} \Delta R\_m \right\} + Z\_{m,n}^I \sin\left\{\frac{4\pi f\_{m,n}}{c} \Delta R\_m \right\} \tag{A7}$$

$$\epsilon\_{m,n}^{l} = Z\_{m,n}^{l} \cos\left\{\frac{4\pi f\_{m,n}}{c} \Delta R\_m \right\} - Z\_{m,n}^{R} \sin\left\{\frac{4\pi f\_{m,n}}{c} \Delta R\_m \right\}.\tag{A8}$$

Using (46)–(A8), the minimization (45) becomes

$$
\overrightarrow{\Delta R}\_{m}(i+1) = \underset{\overrightarrow{\Delta R}\_{m}}{\text{arg min}} \|\mathfrak{e}(\Delta R\_{m})\|\_{2}^{2}.\tag{A9}
$$

By adopting the Gauss–Newton algorithm, Δ 2*Rm*(*i* + 1) can be computed from Δ 2*Rm*(*i*) as

$$
\widehat{\Delta R}\_{\mathfrak{M}}(i+1) = \widehat{\Delta R}\_{\mathfrak{M}}(i) - \left(f\_{\mathfrak{m}}^T(i)f\_{\mathfrak{m}}(i)\right)^{-1} f\_{\mathfrak{m}}^T(i) \mathfrak{e}(\widehat{\Delta R}\_{\mathfrak{m}}(i)) \tag{A10}
$$

where (Δ 2*Rm*(*i*)) is an estimated version of computed from Δ 2*Rm*(*i*) and *Jm*(*i*) is the Jacobian of (Δ*Rm*) with respect to Δ*Rm* evaluated at Δ*Rm* = Δ 2*Rm*(*i*).

The expression for the Jacobian *Jm* of (Δ*Rm*) with respect to Δ*Rm* is given by

$$J\_m = \left[ \left( J\_m^R \right)^T, \left( J\_m^I \right)^T \right]^T \tag{A11}$$

where

$$J\_m^R = \left[ \dots, J\_{m,n}^R, \dots \right]\_{n=1,\dots,N\_m}^T \tag{A12}$$

$$J\_m^I = \begin{bmatrix} \dots \dots \, J\_{m,n'}^I \dots \end{bmatrix}\_{n=1,\dots,N\_m}^T \tag{A13}$$

and

$$J\_{m,n}^R = \frac{4\pi f\_{m,n}}{c} \left( -Z\_{m,n}^R \sin\left\{\frac{4\pi f\_{m,n}}{c} \Delta R\_m \right\} + Z\_{m,n}^I \cos\left\{\frac{4\pi f\_{m,n}}{c} \Delta R\_m \right\} \right) \tag{A14}$$

$$\left[J\right]\_{m,n}^{l} = -\frac{4\pi f\_{m,n}}{c} \left(Z\_{m,n}^{l} \sin\left\{\frac{4\pi f\_{m,n}}{c} \Delta R\_m \right\} + Z\_{m,n}^{R} \cos\left\{\frac{4\pi f\_{m,n}}{c} \Delta R\_m \right\} \right). \tag{A15}$$

#### **Appendix B. Analysis of the Nonconvexity of the** *l***1-Norm Regularised Optimisation Problem** (31)

From Equation (31)

$$\{\hat{\mathbf{a}}, \hat{\mathbf{A}}\} = \operatorname\*{arg\,min}\_{\mathbf{a}, \mathbf{A}} \left\{ \|\mathbf{S} - \mathbf{A} \mathbf{F} \mathbf{a}\|\|\_{2}^{2} + \mu \|\mathbf{a}\|\|\_{1} \right\},\tag{A16}$$

or

$$\{\hat{\mathfrak{a}}, \hat{\mathsf{A}}\} = \underset{\mathfrak{a}}{\arg\min} \arg\min \left\{ \|\mathfrak{S} - \mathsf{A}\mathsf{F}\mathfrak{a}\|\_{2}^{2} + \mu \|\mathfrak{a}\|\_{1} \right\}.\tag{A17}$$

$$\|\|\mathbf{S} - \Lambda \mathbf{F}\mathbf{a}\|\|\_{2}^{2} = \left\{\mathbf{S} - \Lambda \mathbf{F}\mathbf{a}\right\}^{H} \left\{\mathbf{S} - \Lambda \mathbf{F}\mathbf{a}\right\} \tag{A18}$$

$$\|\|\mathbf{S} - \mathbf{A}\mathbf{F}\mathbf{a}\|\|\_{2}^{2} = \mathbf{S}^{H}\mathbf{S} - 2\Re\left[\left(\mathbf{F}^{H}\boldsymbol{\Lambda}^{H}\mathbf{S}\right)^{H}\mathbf{a}\right] + \mathbf{a}^{H}\mathbf{F}^{H}\mathbf{F}\mathbf{a}.\tag{A19}$$

Only the term (*FH***Λ***HS*)*Hα* is dependant on Δ*R*2, so the objective function is minimised with respect to Δ*R*2 when the correlation between *FH***Λ***HS* and *α* is maximised. Now

$$F^H \Lambda^H \mathbf{S} = F\_1^H \mathbf{S}\_1 + F\_2^H \Lambda\_2^H \mathbf{S}\_2 \tag{A20}$$

And the terms *FH* 1 *S*1 and *FH* 2 **Λ***H* 2 *S*2 represent the conventional pulse compression (i.e., transformation from frequency domain to range domain) for each of the radars with a range offset

Δ*R*2. These are non-sparse and low resolution due to the oversampling in range. So *FH***Λ***H*(Δ*R*2)*S* represents the linear superposition of the two conventional range profiles. When these range profiles from the two radars are correctly aligned in range, they will better correlate with the true range profile *α*. Also if the range offset Δ*R*2 is such that a scatterer for one radar is superimposed upon a different scatterer for the other radar, a local minimum will occur. That explains Figure 6 and the reason for its non-convexity.

#### **Appendix C. Analysis of Operation Count for POMP**

Consider first the case of *M* = 2. Let the size of the grid for the parameter Δ*R*2 be *L*Δ = 2*N*Δ , so that a dictionary is constructed for each of these 2*N*Δ values of Δ*R*2. OMP is implemented with the number of dictionaries halved at each stage; hence, the name "pruned" OMP. At stage *k* of OMP, the number of complex multiplications and divisions for a single dictionary is denoted *Comp*(*k*). Here *k* = 1, ... , *N*Δ + 1 with 2*N*Δ−*k*+1 dictionaries considered at stage *k*. The purpose of this section is to estimate the dependence of the computational cost of POMP on the size of the dictionary *Lx* and the grid size *L*Δ for the parameter Δ*R*2.

With reference to Table 1, the significant costs for POMP are associated with the Identify and Update steps. At each stage of OMP for a given dictionary, the Identify step performs atom/residual correlations which require ∼ *O*(*NmLx*) complex multiplications. The Update step performs a linear least squares estimation requiring Gaussian elimination which, at stage *k* of OMP, has an operation count ∼ *O*(*k*3).

Due to the halving of the number of dictionaries at each stage, the total operational count required until only one dictionary is left (although more OMP steps may be required for that dictionary until the residual is sufficiently small) is therefore of order

$$\sum\_{k=1}^{N\_{\Lambda}+1} \left( N\_{\text{H}} L\_{\text{x}} + k^3 \right) 2^{N\_{\Lambda}-k+1} \tag{A21}$$

or

$$(2^{N\_{\Lambda}+1}-1)N\_{\text{ll}}L\_{\text{x}}+(N\_{\Lambda}+1)^3+2^{N\_{\Lambda}+1}\sum\_{k=1}^{N\_{\Lambda}}k^3z^k\tag{A22}$$

with *z* = 1 2 .

The finite sum ∑*n k*=1 *k*3*zk* is referred to as a low-order polylogarithm, for which a formula may be derived [34]. This formula can be shown to have a leading term of order *n*3*zn*+3 so that the overall operation count for POMP is ∼ *O*((2*N*Δ+1 − 1)*NmLx* + *N*3 Δ). As the size of the grid for Δ*R*2 is *L*Δ = 2*N*Δ , *N*Δ = log *L*Δ log 2 , and the operation count in terms of *L*Δ is ∼ *O L*Δ*NmLx* + log *L*Δ log 2 3 . We see that to leading order the computational cost is proportional to the size of the grid for Δ*R*2.

For the case of general *M* this cost is multiplied by (*M* − 1).

#### **References**

© 2020 Commonwealth of Australia. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Compressive Sensing for Tomographic Imaging of a Target with a Narrowband Bistatic Radar**

**Ngoc Hung Nguyen 1,***∗***, Paul Berry 2 and Hai-Tan Tran 2**

Received: 12 November 2019; Accepted: 9 December 2019; Published: 13 December 2019

**Abstract:** This paper introduces a new approach to bistatic radar tomographic imaging based on the concept of compressive sensing and sparse reconstruction. The field of compressive sensing has established a mathematical framework which guarantees sparse solutions for under-determined linear inverse problems. In this paper, we present a new formulation for the bistatic radar tomography problem based on sparse inversion, moving away from the conventional *k*-space tomography approach. The proposed sparse inversion approach allows high-quality images of the target to be obtained from limited narrowband radar data. In particular, we exploit the use of the parameter-refined orthogonal matching pursuit (PROMP) algorithm to obtain a sparse solution for the sparse-based tomography formulation. A key important feature of the PROMP algorithm is that it is capable of tackling the dictionary mismatch problem arising from off-grid scatterers by perturbing the dictionary atoms and allowing them to go off the grid. Performance evaluation studies involving both simulated and real data are presented to demonstrate the performance advantage of the proposed sparsity-based tomography method over the conventional *k*-space tomography method.

**Keywords:** radar tomography; compressive sensing; sparse reconstruction; bistatic radar; radar imaging; parameter-refined orthogonal matching pursuit (PROMP); orthogonal matching pursuit (OMP); *k*-space tomography; narrowband radar; off-grid compressive sensing

#### **1. Introduction**

Radar imaging has received much attention for several decades, having a wide range of applications in both civilian and military domains [1–3]. In principle, to obtain high-resolution radar images, a wide bandwidth of radar waveform is required for a fine resolution in the range direction, while a large antenna aperture is required for a fine resolution in the cross-range direction. To overcome the physical constraints of the radar aperture size, a synthesized aperture with a much larger size can be formed by exploiting the relative motion between the radar and target. This is, in fact, the main idea behind the synthetic aperture radar (SAR) and inverse SAR (ISAR) [1]. In recent years, there has been an increasing demand on the radio-frequency (RF) electromagnetic spectrum due to rapid advances in radar and communications, and the radar has to compete for spectrums with many different services, including radio and television broadcasting, communications, and radio-navigation [4]. As a result, the constraints on spectrum availability may present severe limits on signal bandwidth, prompting the need for high-resolution imaging techniques using narrowband radars.

As a consequence, Doppler tomography has been considered for narrowband radar imaging [5–12], which is also called "Doppler radar tomography" (DRT). The main idea of DRT is to utilize the information given by the Doppler frequencies induced from the relative radar-target rotational motion to construct an image of the target, which can be conveniently formulated in the slow-time *k*-space [8]. Imaging can also be formulated in the more traditional fast-time *k*-space, the support for which is created by sweeping out the complex samples of the received signal in the angular direction. For a particular transmit frequency, the complex samples for all available aspect angles form a circular arc in the spatial frequency space. The traditional range-Doppler ISAR imaging can be considered as a special case of this fast-time *k*-space technique when a wideband signal is available and the total rotation angle is small enough that the support region can be approximated as being rectangular. The inversion process for image formation has evolved from traditional tools, such as filtered back projection, to the more modern non-uniform fast Fourier transform (NUFFT). The *k*-space radar tomography was also considered in bistatic settings [13,14]. The bistatic radar offers several advantages over a monostatic radar, including higher performance against stealth targets, less vulnerability to jamming, and its covertness.

The main objective of this paper is to present a new tomographic imaging technique for a narrowband bistatic radar based on the framework of compressive sensing and sparse reconstruction. The field of compressive sensing has established a mathematical framework which guarantees sparse solutions for underdetermined linear inverse problems that occur across numerous engineering and mathematical science fields. In particular, this framework has found applications in various radar imaging problems, ranging from moving target indication, ISAR imaging, coherence imaging, multichannel imaging, micro-Doppler imaging, to through-the-wall radar imaging (see, e.g., [15–26]). The key contributions of this paper are summarized as follows.

The remainder of this paper is organized as follows. Section 2 describes the signal model for bistatic radar tomography. Section 3 formulates bistatic radar tomographic imaging as a sparse inversion problem. Section 4 derives the sparse solution based on the PROMP algorithm. Performance studies with simulated and real data are presented in Section 5. The paper ends in Section 6 with some concluding remarks.

#### **2. Signal Model**

Figure 1 shows the geometry for the problem of bistatic radar tomographic imaging under consideration. The transmitter Tx and receiver Rx are located in the far field of the target of interest. The bistatic angle between the transmitter and receiver with respect to the target is denoted as *β*. The transmitter is narrowband, transmitting a continuous waveform at a single frequency *f* (i.e., the wavelength is *λ* = *c*/ *f* , where *c* is the speed of signal propagation). A local target coordinate frame T (*x*1, *x*2, *x*3), which is fixed and rotated with the target, is chosen as the reference frame. Here, both the transmitter and receiver lie on the image plane X (*x*1, *x*2) of the target, and the origin of the frame T is placed at the target rotation centre. It is also assumed that the target rotational speed Ω is constant over the coherent processing interval (CPI) and is accurately estimated a priori.

**Figure 1.** The radar-target geometry of the considered bistatic radar tomographic imaging problem.

The receiver takes one complex-valued scattered signal sample for each rotation angle *θn* = Ω*tn* of the target (with respect to the axis *x*2 of the frame T ) at time *tn*. The expression of the scattered signal sample collected at time *tn* is given by [13,14]

$$s(t\_{\mathbb{H}}) = \int\_{\mathcal{X}\_1} \int\_{\mathcal{X}\_2} \sigma(\mathbf{x}) \exp\left\{-j\frac{2\pi}{c}f\left[\mathbf{R}(t\_{\mathbb{H}}) + 2\cos\left(\frac{1}{2}|\beta|\right)\mathbf{x} \cdot \mathbf{u}(t\_{\mathbb{H}})\right]\right\} d\mathbf{x}\_1 d\mathbf{x}\_2. \tag{1}$$

where *x* = [*x*1, *x*2] *T*, *R* is the total bistatic range between the target centre (or focus point) and the transmitter and receiver, and *σ*(*x*) is the scatterer reflectivity distribution projected onto the image plane. Note that the total bistatic range *R* is, in general, a function of time *tn* because of the target translational motion. In (1), *u* = [*u*1, *u*2] *T* denotes the unit vector along the bisector of the bistatic angle *β*. It is also noted that the radar tomographic imaging problem under interest is considered in the rotating local frame T , and the signal model (1) has already taken into account the rotational motion of the target by the rotation of the unit vector *u*(*tn*) relative to T .

Since *x*1 and *x*2 in the image domain are discrete variables in practical radar imaging applications, the reflectivity function *σ*(*x*) is commonly discretized over the image plane X (*x*1, *x*2) onto a grid of points *xm* for *m* ∈ {1, . . . , *M*}, as

$$\sigma(\mathbf{x}) = \sum\_{m=1}^{M} \sigma\_{\mathfrak{M}} \delta(\mathbf{x} - \mathbf{x}\_{\mathfrak{M}}).\tag{2}$$

Substituting (2) into (1) and assuming that translational motion compensation is accurately accomplished by additional pre-processing, we obtain

$$s(t\_{\rm ll}) = \sum\_{m=1}^{M} \sigma\_{\rm m} \exp\left\{-j\frac{4\pi}{c}f\cos\left(\frac{1}{2}|\beta|\right)\mathbf{x}\_{\rm m} \cdot \mathbf{u}(t\_{\rm tr})\right\}.\tag{3}$$

A compact vector form of (3) is given by

$$\mathbf{s} = \Phi \boldsymbol{\sigma},\tag{4}$$

where

$$\mathbf{s} = \begin{bmatrix} s(t\_1), \dots, s(t\_N) \end{bmatrix}^T \tag{5}$$

$$\sigma = [\sigma\_1, \dots, \sigma\_M]^T \tag{6}$$

and

$$\Phi = [\Phi(\mathbf{x}\_1), \dots, \Phi(\mathbf{x}\_M)] \tag{7a}$$

$$\boldsymbol{\Phi}(\mathbf{x}\_{m}) = \begin{bmatrix} \boldsymbol{\Phi}(\mathbf{x}\_{m}, t\_{1}), \dots, \boldsymbol{\Phi}(\mathbf{x}\_{m}, t\_{N}) \end{bmatrix}^{T} \tag{7b}$$

$$\phi(\mathbf{x}\_{\mathrm{ll}}, t\_{\mathrm{ll}}) = \exp\left\{-j\frac{4\pi}{c}f\cos\left(\frac{1}{2}|\beta|\right)\mathbf{x}\_{\mathrm{ll}} \cdot \boldsymbol{\omega}(t\_{\mathrm{ll}})\right\}.\tag{7c}$$

In practice, where noise is presented, the radar received signal becomes

$$
\ddot{s} = s + \mathfrak{u} = \Phi \sigma + \mathfrak{u},\tag{8}
$$

where *n* = [*n*(*t*1),..., *n*(*tN*)]*T*, with *n*(*tn*) denoting the complex-valued noise term at time *tn*.

#### **3. Sparse Inversion Formulation of Bistatic Radar Tomography**

The objective of any target imaging problem is to construct a spatial reflectivity map of the target from the backscattered radar signal. Specifically, the ultimate objective is to estimate the unknown reflection vector *σ* from the noisy received signal vector *s*˜ by solving (8). Since the number of signal samples received, as often is the case in practice, is much smaller than the number of grid points in the reflectivity map (i.e., *N M*), solving (8) is essentially an underdetermined linear inverse problem which requires additional regularization constraints to obtain meaningful solutions.

Typical target images captured by microwave radar signals has been known in the literature to contain a few dominant scattering centers (see, e.g., [15–17,19–23]). As a result, the reflection vector *σ* only has a small number of non-zero elements, thus enjoying a sparse characteristic. Such a sparse characteristic of *σ* can be utilized as a regularizing constraint to solve the underdetermined inverse problem (8), that is,

$$\text{A fixed sparse } \sigma \text{ such that } \mathfrak{s} \approx \Phi \sigma. \tag{9}$$

This sparse inversion problem can be effectively solved under the compressive sensing framework using sparse reconstruction algorithms. Note that in the compressive sensing context, the matrix **Φ** is commonly referred to as the dictionary, and the columns of **Φ** are called the atoms, each representing the theoretical scattered signal component of a hypothetical scatterer residing on a grid point of the reflectivity map.

Compressive sensing and sparse reconstruction have been extensively studied in the last two decades, with various techniques proposed. Comprehensive surveys of the state-of-the-art on this topic can be found in [22,23,30,31]. The objective of this paper is to apply the sparse reconstruction approach to the bistatic radar tomographic imaging problem.

The main challenge for this work is that the true scatterers constituting the target do not coincide exactly with the grid which is used to construct the dictionary, leading to dictionary mismatch problems which in turn significantly degrade the performance of conventional sparse reconstruction techniques [27,28]. Several methods have been presented in the literature to address the off-grid dictionary mismatch problems based on the ideas of joint-sparse recovery [32], dictionary perturbation [33,34], sparse Bayesian learning [35,36], and parameter perturbation [29,37]. In this paper, we will exploit the use of the PROMP method [29,37], that is, a parameter perturbation method, to solve the sparse inversion problem (9). The main motivations of using PROMP are twofold. Firstly, PROMP is capable of tackling the dictionary mismatch problem by perturbing the dictionary atoms and allowing them to go off the grid. Secondly, PROMP is computationally efficient and thus suitable for real-time operation because it belongs to the greedy pursuit family which identifies the support of the solution in an iterative manner based on the level of correlation between the input data and the dictionary atoms.

#### **4. Parameter-Refined Orthogonal Matching Pursuit**

Table 1 summarizes the overall structure of the PROMP algorithm. As a variant of the greedy pursuit technique, PROMP solves the sparse inversion problem (9) by identifying the support of *σ* in an iterative greedy manner. In particular, it starts with an empty support set Λ[0] = ∅ and sets the signal *s*˜ as the initial signal residual *r*[0] . Like other greedy techniques, one column of **Φ** (corresponding to one atom of the dictionary) that produces the largest correlation with the current signal residual *r*[*i*] is chosen and added to the support set Λ[*i*] in each iteration. However, a unique feature of PROMP is that it allows the dictionary atoms to go off the grid by perturbing their parameters, thus it can overcome the off-grid dictionary mismatch problem. Specifically, the updated step of PROMP, as different to that of other greedy pursuit techniques, not only estimates the coefficients *σ*[*i*] *k* but also determines the positions *x* [*i*] *k* = *x* [*i*] 1,*k*, *x* [*i*] 2,*k T* , *k* = 1, ... , *i*, of the scatterers associated with the current support set Λ[*i*] via the least-square sense as

$$\begin{aligned} \left\{ \left. \boldsymbol{\sigma}\_{k}^{[i]}, \dot{\boldsymbol{\mathfrak{x}}}\_{k}^{[i]} \right\} \right\}\_{k=1,\ldots,i} &= \arg\min \left\| \bar{\mathbf{s}} - \sum\_{k=1}^{i} \boldsymbol{\sigma}\_{k}^{[i]} \boldsymbol{\Phi} \left( \mathbf{x}\_{k}^{[i]} \right) \right\|\_{2} \\ \text{subject to} & \\ \left\| \left. \dot{\boldsymbol{\mathfrak{x}}}\_{1,k}^{[i]} - \dot{\boldsymbol{\mathfrak{x}}}\_{1,k}^{[i]} \right\| \leq \zeta \text{ and } \left\| \left. \dot{\boldsymbol{\mathfrak{x}}}\_{2,k}^{[i]} - \dot{\boldsymbol{\mathfrak{x}}}\_{2,k}^{[i]} \right\| \leq \zeta \text{.} \end{aligned} \tag{10} $$

2,*k*

subject to

1,*k*

which is in fact a nonlinear least-square (NLS) estimation problem. Here,  $\bar{x}\_k^{[i]} = \left[\bar{x}\_{1,k}^{[i]}, \bar{x}\_{2,k}^{[i]}\right]^T$ ,  $k = 1, \dots, k$  are the positions of the dictionary atoms in the current support set  $\Lambda^{[i]}$ . Note that the constraint in (10) ensures the position estimate for each scatterer starting within the vicinity of the corresponding dictionary atom. A nominal resolution of  $\lambda/2$  can be used to set the value of  $\tilde{\iota}$ .

#### INPUT: ˜*s*, **Φ**.

	- Identify:

$$\begin{array}{l} \mathbf{c}^{[i]} = \boldsymbol{\Phi}^{H} \boldsymbol{r}^{[i-1]}\\ \boldsymbol{j}^{[i]} = \operatorname{arg} \max\_{\boldsymbol{j}} \left| \mathbf{c}\_{\boldsymbol{j}}^{[i]} \right| \end{array}$$

quit the iteration **if** *j* [*i*] ∈ Λ[*i*−1]

Λ[*i*] = Λ[*i*−1] ∪ *j* [*i*]

$$\left\{\boldsymbol{\mathfrak{d}}\_{k}^{[i]},\boldsymbol{\mathfrak{x}}\_{k}^{[i]}\right\}\_{k=1,\ldots,i} = \arg\min \left\|\boldsymbol{\mathfrak{s}} - \sum\_{k=1}^{i} \boldsymbol{\sigma}\_{k}^{[i]} \boldsymbol{\Phi} (\boldsymbol{\mathfrak{x}}\_{k}^{[i]})\right\|\_{2}$$

$$\begin{aligned} \text{subject to} \\ \left\|\boldsymbol{\mathfrak{x}}\_{1,k}^{[i]} - \boldsymbol{\mathfrak{x}}\_{1,k}^{[i]}\right\| \leq \zeta \text{ and } \left\|\boldsymbol{\mathfrak{x}}\_{2,k}^{[i]} - \boldsymbol{\mathfrak{x}}\_{2,k}^{[i]}\right\| \leq \zeta \end{aligned}$$

$$r^{[i]} = \mathfrak{s} - \sum\_{k=1}^{i} \partial\_k^{[i]} \Phi \left( \mathfrak{x}\_k^{[i]} \right)$$

**end for**.

*σ*ˆ [*i*] *k* , ˆ*x* [*i*] *k*

OUTPUT: &

 $k = 1, \dots, i$
 $\frac{\text{The superscript } ^H \text{ stands for the \"Herrnition\"intersection\" .}}{}$

Since the NLS problem in (10) does not admit a closed-form solution, in what follows we will derive an iterative solution based on the Gauss–Newton (GN) approach [38]. Note that the scatterer reflection coefficient is a complex-valued variable, while the scatterer position is a real-valued variable. For the sake of convenience, we transform (10) into a NLS problem purely in the real-valued domain. In particular, we re-express the cost function in (10) as

$$\left\|\bar{\mathbf{s}} - \sum\_{k=1}^{i} \sigma\_{k}^{[i]} \boldsymbol{\Phi} \left(\mathbf{x}\_{k}^{[i]}\right)\right\|\_{2} = \left\|\begin{bmatrix} \text{Real}\{\bar{\mathbf{s}}\} \\ \text{Imag}\{\bar{\mathbf{s}}\} \end{bmatrix} - \sum\_{k=1}^{i} \begin{bmatrix} \text{Real}\left\{\sigma\_{k}^{[i]} \boldsymbol{\Phi} \left(\mathbf{x}\_{k}^{[i]}\right)\right\} \\ \text{Imag}\left\{\sigma\_{k}^{[i]} \boldsymbol{\Phi} \left(\mathbf{x}\_{k}^{[i]}\right)\right\} \end{bmatrix}\right\|\_{2} \tag{11}$$

where explicit expressions of Real ! *σ*[*i*] *k φ x* [*i*] *k* " and Imag ! *σ*[*i*] *k φ x* [*i*] *k* " are given by

$$\text{Real}\left\{\sigma\_{k}^{[i]}\Phi\left(\mathbf{x}\_{k}^{[i]}\right)\right\} = \left[\ldots, \left(\text{Real}\{\sigma\_{k}^{[i]}\}\cos\theta\_{k}^{[i]}(t\_{\text{\textquotedblleft}}) - \text{Imag}\{\sigma\_{k}^{[i]}\}\sin\theta\_{k}^{[i]}(t\_{\text{\textquotedblright}})\right), \ldots\right]\_{n=1,\ldots,N'}^{T} \tag{12a}$$

$$\text{Imag}\left\{\boldsymbol{v}\_{k}^{[i]}\boldsymbol{\Phi}\left(\mathbf{x}\_{k}^{[i]}\right)\right\} = \left[\dots\right,\left(\text{Real}\{\boldsymbol{v}\_{k}^{[i]}\}\sin\theta\_{k}^{[i]}(t\_{n}) + \text{Imag}\{\boldsymbol{v}\_{k}^{[i]}\}\cos\theta\_{k}^{[i]}(t\_{n})\right),\dots\right]\_{n=1,\dots,N}^{T}\tag{12b}$$

and

$$\theta\_k^{[i]}(t\_\hbar) = -\frac{4\pi f}{c}\cos\left(\frac{|\beta|}{2}\right)\left(x\_{1,k}^{[i]}u\_1(t\_\hbar) + x\_{2,k}^{[i]}u\_2(t\_\hbar)\right). \tag{13}$$

Now we define

$$\bar{z} = \begin{bmatrix} \text{Real}\{\bar{\mathfrak{s}}\} \\ \text{Imag}\{\bar{\mathfrak{s}}\} \end{bmatrix}, \quad z = \sum\_{k=1}^{i} \begin{bmatrix} \text{Real}\left\{\sigma\_{k}^{[i]} \boldsymbol{\Phi}\left(\mathbf{x}\_{k}^{[i]}\right)\right\} \\ \text{Imag}\left\{\sigma\_{k}^{[i]} \boldsymbol{\Phi}\left(\mathbf{x}\_{k}^{[i]}\right)\right\} \end{bmatrix} \tag{14}$$

and

$$\mathfrak{g}^{[i]} = \left[\mathfrak{f}\_1^{[i]T}, \dots, \mathfrak{f}\_i^{[i]T}\right]^T \\ \text{with } \mathfrak{f}\_k^{[i]} = \left[\sigma\_{\mathbb{R},k'}^{[i]} \sigma\_{I,k'}^{[i]} \mathbf{x}\_{1,k'}^{[i]} \mathbf{x}\_{2,k}^{[i]}\right]^T. \tag{15}$$

Here, *σ*[*i*] *R*,*k* = Real{*σ*[*i*] *k* }, *σ*[*i*] *I*,*k* = Imag{*σ*[*i*] *k* }. Noting that *z* is a function of *ξ*[*i*] , (10) is equivalent to

$$\mathfrak{E}^{[i]} = \underset{\mathfrak{E}^{[i]}}{\arg\min} \left\| \bar{z} - z\left(\mathfrak{E}^{[i]}\right) \right\|\_{2}.\tag{16}$$

This is a NLS problem solely in the real-valued domain, and its solution can be obtained via the following GN iteration [38]

$$\mathfrak{f}^{[i]}(h+1) = \mathfrak{f}^{[i]}(h) + \left(\Gamma^T(h)\Gamma(h)\right)^{-1}\Gamma^T(h)\left(\overline{z} - z\left(\mathfrak{f}^{[i]}(h)\right)\right) \tag{17}$$

for *h* = 0, 1, ... , where **Γ**(*h*) = **Γ** ˆ *ξ*[*i*] (*h*) is the Jacobian matrix of *z* with respect to *ξ*[*i*] evaluated at *ξ*[*i*] = ˆ *ξ*[*i*] (*h*) and *z* ˆ *ξ* is an estimate of *z* calculated at *ξ*[*i*] = ˆ *ξ*[*i*] (*h*).

The expression of the Jacobian matrix **Γ**(*ξ*[*i*] ) is given by

$$\Gamma = \left[ \dots, \Gamma\_k, \dots \right]\_{k=1,\dots,j \atop (\Delta) \dots (\Delta)} \tag{18a} \tag{18a}$$

$$
\Gamma\_k = \begin{bmatrix}
\Gamma\_k^{(1)} & \Gamma\_k^{(2)} & \Gamma\_k^{(3)} & \Gamma\_k^{(4)} \\
\Gamma\_k^{(5)} & \Gamma\_k^{(6)} & \Gamma\_k^{(7)} & \Gamma\_k^{(8)}
\end{bmatrix} \tag{18b}
$$

$$\Gamma\_k^{(1)} = \begin{bmatrix} \dots \dots \cos \theta\_k^{[i]}(t\_n), \dots \end{bmatrix}\_{n=1,\dots,N}^T \tag{18c}$$

$$\begin{aligned} \mathbf{T}\_k^{(5)} &= \left[ \dots , \sin \theta\_k^{[i]} (t\_n) , \dots \right]\_{n=1, \dots, N}^T \\ \mathbf{T}\_k^{(2)} &= \left[ \dots , -\sin \theta\_k^{[i]} (t\_n) , \dots \right]\_{n=1 \quad \dots \quad \dots}^T \end{aligned} \tag{18d}$$

$$\begin{aligned} \boldsymbol{\Gamma}\_{k}^{(2)} &= \left[ \dots , -\sin \theta\_{k}^{[\boldsymbol{i}]} (\boldsymbol{t}\_{\boldsymbol{n}}) , \dots \right]\_{\boldsymbol{n}=1,\dots,N}^{T} \\ \boldsymbol{\Gamma}\_{k}^{(6)} &= \left[ \dots , \cos \theta\_{k}^{[\boldsymbol{i}]} (\boldsymbol{t}\_{\boldsymbol{n}}) , \dots \right]\_{\boldsymbol{n}=1,\dots,N}^{T} \end{aligned} \tag{18}$$

$$\mathbf{T}\_{k}^{(3)} = \left[ \dots , \frac{4\pi f}{c} \cos\left(\frac{|\beta|}{2}\right) u\_1(t\_n) \left( \sigma\_{k,k}^{[i]} \sin\theta\_k^{[i]}(t\_n) + \sigma\_{l,k}^{[i]} \cos\theta\_k^{[i]}(t\_n) \right) \dots \right]\_{n=1,\dots,N}^{T} \tag{18g}$$

$$\Gamma\_{k}^{(7)} = \left[ \dots, \frac{4\pi f}{c} \cos\left(\frac{|\beta|}{2}\right) u\_1(t\_n) \left( -\sigma\_{R,k}^{[i]} \cos\theta\_k^{[i]}(t\_n) + \sigma\_{I,k}^{[i]} \sin\theta\_k^{[i]}(t\_{\text{il}}) \right), \dots \right]\_{n=1,\dots,N}^{T} \tag{18h}$$

$$\Gamma\_k^{(4)} = \left[ \dots, \frac{4\pi f}{c} \cos\left(\frac{|\beta|}{2}\right) u\_2(t\_n) \left( \sigma\_{R,k}^{[i]} \sin\theta\_k^{[i]}(t\_n) + \sigma\_{I,k}^{[i]} \cos\theta\_k^{[i]}(t\_n) \right), \dots \right]\_{n=1,\dots,N}^T \tag{18i}$$

$$\mathbf{T}\_{k}^{(8)} = \left[ \dots , \frac{4\pi f}{c} \cos\left(\frac{|\beta|}{2}\right) \boldsymbol{u}\_{2}(t\_{n}) \left( -\sigma\_{\mathbf{R},k}^{[i]} \cos\theta\_{k}^{[i]}(t\_{n}) + \sigma\_{I,k}^{[i]} \sin\theta\_{k}^{[i]}(t\_{n}) \right) , \dots \right]\_{n=1,\dots,N}^{T} . \tag{18}$$

The following decision logic is then applied at the end of each GN iteration to ensure the constraint in (10) is met:

$$\begin{split} \text{if } \left( \hat{\mathfrak{x}}\_{1,k}^{[i]}(h+1) - \bar{\mathfrak{x}}\_{1,k}^{[i]} \right) &\stackrel{\text{>}}{\gtrless} \pm \zeta \text{ set } \hat{\mathfrak{x}}\_{1,k}^{[i]}(h+1) = \bar{\mathfrak{x}}\_{1,k}^{[i]} \pm \zeta, \\ \text{if } \left( \hat{\mathfrak{x}}\_{2,k}^{[i]}(h+1) - \bar{\mathfrak{x}}\_{2,k}^{[i]} \right) &\stackrel{\text{>}}{\lesssim} \pm \zeta \text{ set } \hat{\mathfrak{x}}\_{2,k}^{[i]}(h+1) = \bar{\mathfrak{x}}\_{2,k}^{[i]} \pm \zeta. \end{split} \tag{19}$$

When the constraint is in effect, a re-estimation of the reflection coefficients is performed as

$$\left[\mathfrak{d}\_1^{[i]}(h+1), \dots, \mathfrak{d}\_i^{[i]}(h+1)\right]^T = \left(\mathbf{Y}^H\mathbf{Y}\right)^{-1}\mathbf{Y}^H\mathbf{\tilde{s}}\tag{20}$$

with **Υ** = *φ*(*x*ˆ [*i*] 1 (*h* + 1)), ... , *φ*(*x*ˆ [*i*] *i* (*h* + 1)) . The GN iteration can be halted after a fixed number of iterations or if the *l*2 norm of the updating term falls below a given threshold.

The GN iteration is initialized to the solution ˆ *ξ*[*i*−1] obtained from the previous PROMP iteration *i* − 1 and the newly selected atom:

$$\mathfrak{z}^{[i]}(0) = [\mathfrak{z}^{[i-1]T}, \text{Real}\{\boldsymbol{\vartheta}\_{j[i]}\}, \text{Imag}\{\boldsymbol{\vartheta}\_{j[i]}\}, \bar{\mathfrak{x}}\_{j[i]}^T]^T \tag{21}$$

where *x*¯*j* [*i*] is the position of the newly selected atom at index *j* [*i*] within the dictionary and *σ*¯*j* [*i*] = *φH*(*x*¯*j* [*i*])*φ*(*x*¯*j* [*i*]) −1 *φH*(*x*¯*j* [*i*])*r*[*i*−1] is the corresponding initial coefficient estimate for this atom.

#### **5. Results**

In this section, we demonstrate the performance superiority of the proposed sparsity-based tomography method based on the PROMP algorithm over the conventional *k*-space tomography method via results using both simulated and real data. The result comparison also includes the performance of the OMP algorithm to illustrate the off-grid dictionary mismatch problem of the sparsity-based tomography formulation and to verify the effectiveness of PROMP in dealing with this issue.

#### *5.1. Results with Simulated Data*

We consider two synthetic targets with two and eight scatterers, respectively, as depicted in Figure 2. In this simulation, the target rotational speed is set to Ω = 37.70 rad/s and the signal frequency is set to *f* = 9.96 GHz. The constraint value *ζ* is set to *ζ* = *λ*/2 for PROMP. The dictionary matrix is constructed using a regularly spaced grid in Cartesian coordinates with the grid step size of Δ = *λ*/5 on each *x*- and *y*-axis. The sampling frequency at the receiver is set to 2.16 kHz. PROMP and OMP iterations are halted if the signal residual reaches the noise level.

**Figure 2.** Ground truth images of two synthetic targets under consideration.

Figure 3 compares the reconstructed images for synthetic target 1 obtained by the *k*-space, OMP and PROMP algorithms for various noise levels. Here, the number of data samples is *N* = 360 (i.e., the CPI approximately being one full rotation cycle of the target). Note that the SNR is defined by SNR = 20 log10( *s* / *n* ). We observe that the *k*-space technique produces images with two main peaks corresponding to the true target scatterers. However, along with these two main peaks, the images obtained by the *k*-space technique also contain other sidelobes with many spurious peaks. Specifically, the higher the noise is, the poorer the performance of the *k*-space technique (i.e., yielding a larger numbers of spurious peaks). In contrast, such a problem associated with spurious peaks does not appear in the OMP and PROMP images, thus demonstrating the performance advantage of the sparsity-based tomography approach over the conventional *k*-space approach. It is observed that one of the scatterers is split into multiple peaks in the OMP images. This observation can be

explained by the fact the OMP solution relies on the fixed dictionary which is built based on a grid of atoms while the true scatterers of the target do not coincide with this dictionary grid, thereby demonstrating the dictionary mismatch problem. On the other hand, by perturbing the dictionary atoms and allowing them to go off the grid, the PROMP algorithm can effectively overcome the dictionary mismatch problem by exhibiting a clean image with only two peaks corresponding to the true scatterers. More importantly, the locations of the peaks in the PROMP images almost exactly match the locations of the true scatterers.

**Figure 3.** Reconstructed images obtained by the *k*-space, OMP and PROMP algorithms for synthetic target 1.

Figure 4 shows the results for synthetic target 2. This is a more challenging scenario because target 2 contains much more scatterers than target 1. We observer that the conventional *k*-space method is struggling to produce reliable image results because of the interaction between the sidelobes of different main peaks, especially in large noise scenarios. Such an interaction leads to some strong spurious peaks which have similar magnitudes to the correct peaks that correspond to the true scatterers, thus making the resulting images severely distorted. On the other hand, compared to the *k*-space method, OMP results in much more satisfactory images. However, the OMP performance is significantly affected by the dictionary mismatch problem arising from off-grid scatterers. As a result, the true scatterers are split into multiple peaks in the OMP images, and some spurious peaks also appear. In contrast, PROMP produces clean and clear images which are almost identical to the ground truth target image, even at large noise levels.

**Figure 4.** Reconstructed images obtained by the *k*-space, OMP and PROMP algorithms for synthetic target 2.

To further demonstrate the superior performance of PROMP, the reconstructed images obtained from less data samples (i.e., with CPI = 2/3 and 1/3 target rotation cycle) are shown in Figure 5. With a limited number of data samples, the *k*-space method results in unsatisfactory images with incorrect peaks, while OMP and PROMP are observed to retain their good performance. In addition,

similar to the observations in Figure 4, PROMP outperforms OMP and provides a better image of the target thanks to its ability to deal with off-grid scatterers.

**Figure 5.** Reconstructed images obtained by the *k*-space, OMP and PROMP algorithms for synthetic target 2, given less data samples.

Note that, to satisfy the sparsity condition, the number of dominant scatterers constituting the target must be sufficiently small compared to the number of grid points on the reflectivity map. The required sparsity level in general depends on several factors, including the number of data samples, the noise level, as well as the level of coherence between the atoms of the dictionary. In compressive sensing, the restricted isometry property and the mutual incoherence property establish theoretical connections between those factors required for the effectiveness of sparse reconstruction [23]. However, these analytical metrics are overly-conservative and do not reflect the average performance which is often of interest from the practical point of view [39].

We now compare the performance of the *k*-space, OMP, and PROMP methods using the earth mover's distance (EMD) between the true and reconstructed images. EMD [40] is a widely-used metric to compare the similarity between different images. In principle, EMD is an estimate of the distance between two distributions which is equivalent to the minimal amount of work required for one distribution to be transformed to the other [40]. Figure 6 shows the EMD performance of the *k*-space, OMP, and PROMP methods, averaged from 1000 Monte Carlo runs, against various levels of SNR for the synthetic target 2 and CPI = 1 target rotation cycle. We observe that the PROMP method exhibits an EMD much smaller than those of the *k*-space and OMP methods. This indicates that, from a statistical point of view, the image obtained by PROMP is much closer to the ground-truth image than those obtained by the *k*-space and OMP methods, thus verifying the performance superiority of the PROMP method.

**Figure 6.** EMD performance of the *k*-space, OMP and PROMP algorithms versus SNR for synthetic target 2 and CPI = 1 target rotation cycle.

Table 2 compares the runtimes of the *k*-space, OMP, and PROMP algorithm for the image results shown in Figure 4b. For a fair comparison, all methods were implemented in MATLAB on the same Intel Core i7 3.40 GHz CPU with 16 GHz RAM. We observe that the *k*-space method is much slower than the OMP and PROMP methods. The reason for this is that the *k*-space method requires the non-uniform fast Fourier transform to be performed, thus being computationally more demanding compared to OMP and PROMP. On the other hand, the OMP and PROMP methods are computationally fast, thanks to the fact that they belong to the greedy pursuit family. In addition, each iteration of PROMP is about 7.3 times slower than each iteration of OMP because PROMP incorporates a NLS solver, rather than a linear least-squares solver as in OMP. However, the overall timerun of PROMP is only 2.4 times slower than that of OMP because PROMP requires many fewer iterations than OMP to make the signal residual reach the noise level.

**Table 2.** Complexity comparison .

 All methods are implemented in MATLAB on an Intel Core i7 3.40 GHz CPU with 16 GHz RAM, for the image results shown in Figure 4b.

#### *5.2. Results with Real Data*

The experimental data used in this paper was collected in the Mumma Radar Laboratory at the University of Dayton, Ohio, USA. Although this paper focuses on narrowband tomographic imaging, a wideband waveform at X-band with stepped frequency pulses over 101 regular frequency steps from 8 GHz to 12 GHz was used in the experiment. The aim of using wideband data was only for the removal of extraneous clutter components existing in the lab environment, as described in [13]. After that, only the measured data from one discrete frequency is actually used for algorithm performance evaluation for the problem of narrowband tomographic imaging under consideration.

Figure 7 shows photos of the experimental system configuration. The experimental setup involves transmitting and receiving horn antennas mounted on separate robotic arms in a controlled laboratory environment with radar-absorbing material to reduce the radar reflections from the floor and walls. These robotic arms could be oriented and positioned with high precision. During the experiment, the antennas were kept stationary while the target was rotated over 1◦ steps through 360◦, where the stepped-frequency waveform was transmitted and sampled (one sample for each frequency). Recall from above that only one frequency sample set was used for tomographic imaging purposes. The experimental target was comprised of two vertical metallic rods with a 19 cm separation to emulate two point scatterers which rotate around a vertical pedestal. Figure 8 shows the wideband *k*-space tomographic image obtained from monostatic data with all 101 frequency steps. This will be used as a reference benchmark for the performance evaluation of narrowband bistatic imaging presented in this section.

**Figure 7.** Photos of experimental system configuration: (**a**) an antenna mounted on a robotic arm, and (**b**) two vertical metallic rods secured to a rotating pedestal.

G%

**Figure 8.** A benchmark experimental target image using wideband signals received by the monostatic receiver for all 101 frequency steps between 8 GHz and 12 GHz.

Figure 9 shows the experimental narrowband images obtained by the *k*-space, OMP, and PROMP algorithms using a narrowband signal received by the bistatic receiver with *β* = 86◦ at a single frequency of 8.8 GHz for various values of CPI. Since the SNR is unknown, the OMP and PROMP iterations are halted when the change in the signal residual norm falls below 1% of the input signal norm. Compared to the benchmark image in Figure 8, the *k*-space method produces images with a much lower quality when only a narrowband signal from a single frequency is available. We observe that the images obtained by the *k*-space method are distorted with numerous spurious sibelode peaks. In particular, the number of spurious sibelode peaks increases significantly for shorter CPI. Moreover, the main lodes corresponding to the true scatterers are also spread when the CPI is reduced. In contrast, the PROMP images contain two clear peaks at the locations very close to the peaks of the reference image in Figure 8, even when the CPI is reduced to one third of the target rotation cycle. This observation demonstrates the performance superiority of the sparsity-based tomography approach over the conventional *k*-space tomography approach. Figure 9 also shows the images obtained by OMP to illustrate the dictionary mismatch problem associated with off-grid scatterers, where each true scatterer is split into multiple peaks; thus, verifying the effectiveness of PROMP in terms of tackling the dictionary mismatch problem.

**Figure 9.** Experimental narrowband images obtained by the *k*-space, OMP, and PROMP algorithms using a narrowband signal received by the bistatic receiver with *β* = 86◦ at a single frequency of 8.8 GHz, for various values of CPI.

Figure 10 shows the experimental results where the data is downsampled. Here, we observe a similar relative performance comparison to Figures 3–5 and 9, once again confirming the performance advantages of the proposed sparsity-based tomographic imaging method based on the PROMP algorithm.

**Figure 10.** Experimental narrowband images obtained by the *k*-space, OMP, and PROMP algorithms for downsampled data (CPI = 1 target rotation cycle).

#### **6. Conclusions**

In this paper, we have proposed a new sparsity-based bistatic radar tomographic imaging method exploiting the use of the PROMP algorithm. A new formulation for radar tomography building on the framework of compressive sensing and sparse reconstruction was presented, moving away from conventional *k*-space tomography which is prone to sidelobe responses and their interference. The PROMP algorithm was adopted to obtain a sparse solution for the resulting sparsity-based tomography formulation. By perturbing the dictionary atoms and allowing the estimated scatterers to go off the grid, PROMP is capable of tackling the dictionary mismatch problem arising from off-grid scatterers. The performance

advantages of the proposed sparsity-based tomography method over the conventional *k*-space tomography method were demonstrated via numerical studies involving both simulated and real data.

**Author Contributions:** Conceptualization, N.H.N., H.-T.T. and P.B.; Methodology, N.H.N., H.-T.T. and P.B.; Simulation, N.H.N.; Validation, N.H.N., H.-T.T. and P.B.; Writing—Original Draft Preparation, N.H.N.; Writing—Review & Editing, N.H.N., H.-T.T. and P.B.

**Funding:** This research was funded by Defence Science and Technology Group, Australia.

**Acknowledgments:** Preliminary results of this work were reported in [41], where the feasibility and benefits of using compressive sensing and sparse reconstruction for narrowband bistatic radar tomographic imaging were demonstrated. This paper presents detailed algorithmic design and implementation as well as experimental evaluation of the PROMP-based tomographic imaging method.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 Commonwealth of Australia. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Two-Dimensional Augmented State–Space Approach with Applications to Sparse Representation of Radar Signatures**

#### **Kejiang Wu and Xiaojian Xu \***

School of Electronics and Information Engineering, Beihang University, Beijing 100191, China<EMAIL_ADDRESS>
**\*** Correspondence<EMAIL_ADDRESS>Tel.: +86-10-8231-6065

Received: 16 September 2019; Accepted: 21 October 2019; Published: 24 October 2019