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

In this work, we focus on sparse representation of two-dimensional (2-D) radar signatures for man-made targets. Based on the damped exponential (DE) model, a 2-D augmented state–space approach (ASSA) is proposed to estimate the parameters of scattering centers on complex man-made targets, i.e., the complex amplitudes and the poles in down-range and aspect dimensions. An augmented state–space approach is developed for pole estimation of down-range dimension. Multiple-range search strategy, which applies one-dimensional (1-D) state–space approach (SSA) to the 1-D data for each down-range cell, is used to alleviate the pole-pairing problem occurring in previous algorithms. Effectiveness of the proposed approach is verified by the numerical and measured inverse synthetic aperture radar (ISAR) data.


Introduction
Sparse representation of two-dimensional (2-D) radar signatures has been widely used in many applications, such as super-resolution radar imaging, data compression, and target identification [1][2][3][4]. 2-D radar signatures can be reconstructed with fewer data, where a set of parameters including the locations, amplitudes, and damping factors are used to represent the returned signals from spatial distributed scattering centers. Moreover, 2-D ultra-wideband images of radar targets can be obtained through parameter interpolation or extrapolation [5].
Over the years, several model-based spectral estimation approaches have been developed for sparse representation of 2-D radar signatures, in which the common idea of these approaches is the development of a parametric model based on electromagnetic scattering mechanisms. In this way, a set of parameters may be found to represent the original signatures. Please note that most sparse representation approaches are considered to have super-resolution capabilities because they are supposed to estimate the parameters of scattering centers that cannot be distinguished by standard processing [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. Typically, these approaches include the amplitude and phase estimation of a sinusoid (APES) algorithm [6], fast Fourier transform (FFT)-based technique (CLEAN) [7,8], the compressed sensing (CS)-based procedure [9][10][11] and the subspace-based approach [12][13][14][15][16][17][18][19][20][21]. The APES algorithm could provide accurate estimation of the complex amplitudes whereas it does not estimate the locations of scattering centers. The 2-D CLEAN technique is a computationally efficient procedure which uses the undamped exponential model and deconvolution algorithm for optimization [7,22]. The main process of the CS-based procedure is to represent the signatures using an over-complete dictionary and the corresponding coefficients [9,10]. The performance of this approach depends on the initial coefficients and the threshold used in terminating iterations.  (1) . Pole denotes a transfer function of the scattering center collected on target. Usually, the scattering centers are characterized via pairing the complex amplitudes and the poles in down-range and aspect dimensions.
While these subspace-based algorithms are demonstrated to be useful for the signatures of targets or simulated scattering centers, there are still some technique challenges, especially for sparse representation of the wideband radar signatures collected on complex man-made targets. First, the pairing schemes used by many existing algorithms cannot provide correct pole pairs in certain circumstances. i.e., MEMP, ACMP and 2-D system realization meet with the pole-pairing problem when there are numerous repeated poles in either the down-range or aspect dimension [17,20,23]. Second, the model order, which directly affects the result of parameter estimation such as 2-D MUSIC, MEMP, and 2-D ESPRIT et al., is difficult to be determined in the measurement environment [8]. Third, the compromises between the computational complexity and accuracy should be considered in the application of these algorithms [12]; for instance, the 2-D TLS Prony has smaller computational complexity than MEMP but with loss in accuracy [15].
In this paper, we focus on developing an approach for sparse representation of 2-D radar signatures collected on man-made targets. An augmented state-space approach is proposed for pole estimation of the down-range dimension. Multiple-range search strategy is then applied to estimate the pairing poles and corresponding amplitudes along the aspect dimension. Compared to the existing methods [14,17,18], the advantages of the 2-D augmented state-space approach (ASSA) are: (1) Computational complexity of the algorithm is much reduced since the newly defined Hankel matrix and several time-saving operations are adopted, whereas the pole-estimation accuracy is still at the same level; (2) The pole-pairing problem can be alleviated because all the poles are adaptively paired by using the multiple-range search strategy; And (3) an eigenvalue sequences transform algorithm is proposed, which could provide fast model order selection.
The remainder of the paper is organized as follows. Section 2 briefly presents the damped exponential models for 1-D and 2-D signals. In Section 3, a two-step procedure for 2-D ASSA is developed. Results for numerical and measured inverse synthetic aperture radar (ISAR) data representation are demonstrated to validate the effectiveness of the proposed procedure in Section 4. We conclude the paper in Section 5. Appendices are given to show more mathematical details of the proposed approach.

1-D Damped Exponential Model
As described in DE model [22,24], the 1-D radar signatures y( f n ) can be expressed as a summation of K scattering centers corrupted with noise w(n).
where n = 1,2, . . . ,N and N denotes the number of pulses, A k is the complex amplitude of the k-th scattering center; β k is the damping factor with respect to frequency; r k denotes the relative range; The parameter f n denotes the radar frequency f n = f c + (n − 1 − N2)∆ f where f c is the center frequency and N2 = ceil(N/2) denotes the smallest integer less than or equal to N/2; a k represents the amplitude of the k-th scattering center in pole form, the pole p k = exp[(β k + j2πr k /c)∆ f ] represents the transfer function of the k-th scattering center; c = 3 × 10 8 m/s is the propagation velocity. It worth noting that the data models used in this paper are all considered in a stepped frequency radar [25,26]. According to the discrete-time control theory and auto-regressive moving average (ARMA) model, Piou and Naishadham proposed a one-dimensional state-space approach (1-D SSA) which use a state-space description to the 1-D radar signatures in (1) [22].
where x(n) ∈ C K×1 is the state vector, y(n) denotes the signal sequence y( f n ), u(n) is the input vector, As described in [22], 1-D SSA could precisely estimate the state matrices A, B, and C. Once these three state matrices are computed, the model parameters in (1) can be estimated by using the eigen-decomposition technique.

2-D Damped Exponential Model
As a 2-D extension of 1-D DE model, the radar signatures obtained from different aspect angles are considered to be the summation of a finite number of dispersive scattering centers [13][14][15][16][17][18]. Typically, it is applicable to modeling the 2-D radar signatures with small aspect ranges [27].
a k s m k p n k + w(m, n) in matrix notation, the 2-D radar signatures in (5) can be expressed as: where m = 1, . . . , M and n = 1, . . . , N; {r 1k ,r 2k } give the relative locations of the k-th scattering center; {β 1k , β 2k } characterize the frequency and aspect dependence of scattering; θ m = θ 0 + (m − 1)∆θ denotes the m-th aspect angle where θ 0 is the starting angle and ∆θ represents the angle interval; w(m, n) is the Gaussian noise with zero-mean; y(m, n) denotes the signal sequences y(θ m , f n ); {s k , p k } refer to poles of the down-range and aspect dimension.

Two-Dimensional Augmented State-Space Approach
From the 2-D DE model in (5), the vector Y(n), which represents the n-th column of the noiseless signature matrix Y = Y − W (where W is the noise matrix), can be decomposed as.
where K denotes the number of scattering centers, l K indicates the column vector of ones with length K. Actually, we often meet the one-to-multiple matching situation, i.e., one pole p k is corresponding to more than one poles s k1 , s k2 , . . .. These repeated poles can be merged and (9) could be rewritten as.
where K 1 represents the number of non-repeated poles in the down-range dimension. l 1 , l 2 , . . . l K 1 denote the number of repeated poles, respectively, for the corresponding poles p 1 , p 2 , . . . p K 1 .
We define the matrix P as.
where * denotes the Hermitian operator; Q is a unitary matrix which satisfies Q * Q = QQ * = E, E is the identity matrix.
Two constant matrices, i.e., S and D, are defined to simplify the expression in (10).
Thus, the noiseless signatures Y can be written as.
Comparing (15) with (4), it can be found that these two equations have extremely similar structures. Thus, the parameters a k , s k and p k can be estimated by 2-D augmentation of 1-D SSA, which is called 2-D augmentation state-space approach (2-D ASSA).
Here, the 2-D ASSA consists of two steps. An augmented state-space approach is applied to estimating the pole p k in down-range dimension. Then multiple-range search strategy is used for estimating the matching pole s k as well as the corresponding amplitude a k . Details are as follows.

Pole Estimation of the Down-Range Dimension
First, we introduce a single augmented Hankel matrix H of size M(N − L + 1) × L by analogy of the enhanced Hankel matrix in 1-D SSA: where each element Y(n), first mentioned in (9), represents the n-th column of the matrix Y; L is the step size of the correlation window, which is heuristically set to be L=N2 [22]. Next, do the singular value decomposition (SVD) of the Hankel matrix.
where U sn ,U n ,V * sn and V * n are unitary matrices; sn and n are diagonal matrices; The matrices with subscript 'sn' refer to the components of signal space and those matrices with subscript 'n' denote the components of noise space; The rank of the noiseless matrix U sn sn V * sn is known as the model order which has been widely used in [14,22,27].
As proved in [14,17], the model order in the down-range dimension is equal to K 1 and should be less than the number of columns or rows of H at least. Thus, the number of non-repeated poles K 1 should satisfy the following condition.
Model order selection is an inevitable problem in modeling the 2-D signatures [14,22,27]. A series of eigenvalue-based criteria, such as Akaike information criterion (AIC), minimum description length (MDL) criterion and the minimum eigenvalue (MEV) criterion, have been proposed for solving this problem [28,29]. Those criteria are useful in model order selection but with time consuming process [30].
Here an eigenvalue sequences transform algorithm is proposed for estimating the model order due to its low computational complexity. This algorithm can provide fast estimation but with loss in robustness to noise. More details of this algorithm are presented in Appendix A.
Based on the linear systems theory [31] and the matrix expression in (14) and (16), the noiseless matrix H can be further factorized as: where Ω is the (N − L + 1)M × K 1 observability matrix which can be further expressed by the matrices S and P.
and Γ is the K 1 × L controllability matrix which can be expressed by the matrices P and D.
Considering that computational complexity of (19)-(21) is enlarged observably for large data sets, operations (22)-(24) are used for alternative steps which have lower computational load but with minimal calculation error.
where Z = 2 sn . As an analogy of the open-loop matrix in 1-D SSA, the augmented open-loop matrix P can be derived from the observability matrix by using least square.
Or it can be computed by the controllability matrix Γ.
where Ω r f is the first (N − L)M rows of Ω and Ω rl denotes the last (N − L)M rows of Ω; Γ c f represents the first (L − 1) columns of Γ and Γ cl is the last (L − 1) columns of Γ. From (20) and (21), these matrices can be rewritten as More details of the derivation of the matrices S, P, and D are listed in Appendix B. According to (11), the vector [ p 1 p 2 . . . p K 1 ] is obtained by performing the SVD.
where diag(Λ) denotes the diagonal element of matrix Λ.

Pole Estimation of the Aspect Dimension
In this step, multiple-range search, which applies one-dimensional (1-D) state-space approach (SSA) to the 1-D data for each down-range cell, is used for pole estimation of the aspect dimension and pole adaptive pairing. Details are as follows.
We introduce the Vandermonde sub-matrix O.
From (5), the noiseless signatures Y can be factorized as a product of the pairing matrix G and O, where each column of G is associated with each row of O.
Thus, the pairing matrix G in (34) can be calculated by using least square, i.e.
The k-th column of the pairing matrix G is presented as where k = 1, 2, . . . , K 1 . Using (36), each column of the pairing matrix is constructed in a same structure to 1-D DE model defined in (1). This indicates that the pole s kt and a kt , which correspond to the k-th pole p k in (32), can be solved by using 1-D SSA [22] to each column of the pairing matrix G. For example, for the first column of the pairing matrix G corresponds to p 1 , the parameter matrices (C, A, B) can be obtained by using 1-D SSA in (4).
The matching pole s 1t and the corresponding amplitude a 1t are computed as.
Thus, the pairing matrix G can be searched column by column until all poles s kt and the corresponding amplitudes a kt are estimated. No extra pairing scheme is required because the poles [s k1 , s k2 , . . . , s kl k ] and the amplitudes [a k1 , a k2 , . . . , a kl k ] have already been adaptively paired to the poles p k in (34). In addition, considering that the model order may be overestimated sometime, the searched pole pairs are usually checked again according to their amplitudes. Finally, locations and damping factors of all the scattering points can be obtained from (41) and (42).

Results and Discussion
In this Section, three examples are presented to demonstrate the usefulness of the proposed procedure, i.e., the numerical signatures with 14 point scattering centers, the numerical signatures of a sphere tipped cone-cylinder-frustum combination model, the measured ISAR data for an aircraft model. The results obtained by 2-D ESPRIT are used for comparison since it is one of the very few techniques which have been used in real radar applications [17,32].

Numerical Signatures with Point Scattering Centers
In this example, the noisy signatures composed of 14-point scattering centers are considered to be as follows:  As shown in Figure 1, these scattering points form a missile-like shape in Cartesian coordinates. It contains 13 pole pairs which have repeated poles in either the down-range or aspect dimensions. The noisy signatures in space domain by Fourier transform-based imaging algorithm (add Taylor window) are displayed in Figure 2a,b. As can be seen, the positions and decay rates of these scattering points are consistent with the corresponding ranges and damping factors. To sparse representation of these noisy signatures, the parameter for 2-D ASSA is chosen to be L = 20; the parameters for 2-D ESPRIT are set as: (P, Q) = (20, 20), β= 0.8, which are suggested in [17]. For each signal to noise ratio (SNR), the number of Monte Carlo simulation is 200. A series of range estimation results in different SNRs are presented in Figure 3. Please note that the model order in 2-D ESPRIT is pre-specified because the singular values-based criterion [28,29] cannot be used when there are repeated poles in either the down-range or aspect dimensions. Because of different pairing strategies, K 1 and K in 2-D ASSA could be estimated by the eigenvalue sequences transform algorithm (Appendix A) when SNR > 10 dB. However, the numbers K 1 and K in 2-D ASSA should be pre-specified or use the other criterion [30] when the noise level is higher (SNR ≤ 10 dB).
As can be seen in Figure 3a,b, the positions of scattering points estimated by 2-D ESPRIT are generally according to the preset poles in Figure 1, except for some missing or incorrect points (shown by the black rectangular). The possible reason of this problem is that the pairing procedure in 2-D ESPRIT may provide incorrect pole pairs when there are repeated poles in either the down-range or aspect dimension, i.e., six pole pairs (s 1 , p 1 ), (s 1 , p 2 ), (s 1 , p 3 ), (s 2 , p 1 ), (s 2 , p 2 ), (s 2 , p 3 ). In contrast, the results estimated by 2-D ASSA showed higher accuracy than 2-D ESPRIT for different pairing strategies.
The estimation accuracy of cross/down ranges and damping factors are displayed in Figure 4. Root mean square error (RMSE) is used as the evaluating indicator which is defined in (44). As we can see, estimation accuracy of these two algorithms are basically at the same level although the size of Hankel matrix used by 2-D ASSA is smaller than the block-Hankel matrix defined by 2-D ESPRIT.    δ RMSE = 10 log 10 1 where M c = 2000 denotes the number of Monte Carlo runs, X est and X real are the esimated and real parameters.

Numerical Signatures of Computer-Aided Design (CAD) Model
The numerical signatures are obtained using method of moment (MOM), where the target is from a computer-aided design (CAD) model of a sphere tipped cone-cylinder-frustum combination, (shown in Figure 6). The data was calculated from 8-12 GHz in 10 MHz frequency step size, view angle ranging from −5 to 5 deg with an increment of 0.25 deg. Figure 7a presents 2-D radar image processed using 2-D FFT. As it can be seen, strong scattering points can be observed at differential discontinuities, such as the base-edge, the body groove, and nosetip. To sparsely represent the numerical data, the parameters of this example are set as follows. For 2-D ASSA, L is set as N2. For 2-D ESPRIT, (P, Q) = (M2, N2), where M2 = ceil(M/2) denotes the smallest integer less than or equal to M/2, and β = 0.8.
Location estimation of key scattering points for these two algorithms are shown in Figure 7b-e. As can be seen, when the number of scattering centers is set to be 14, all key scattering points are accurately estimated by these two algorithms except for some minor differences. However, when K is set to be 18, 2-D ESPRIT encountered the pole-pairing problem and could not provide the right estimation. Relative reconstruction error (RRE) δ RRE (defined in (45)) of these two algorithms are shown in Figure 8. We can see that the RRE of 2-D ASSA has been falling when K increased from 10 to 30, whereas the RRE of 2-D ESPEIT stopped falling after K = 14. The result shows that 2-D ASSA tends to be more robust in pole-pairing for different numbers of scattering centers.
where Y recon denotes the signatures reconstructed by the estimated parameters. Y real is the original signatures.    Figure 9a displays the sub-band data image (8-10 GHz, from −5 to 5 deg) which was extracted from the full-band data. Compared to the full-band data image in Figure 7a, scattering centers of the base-edge and the body groove cannot be distinguished when the signal bandwidth is only 2 GHz. Figure 9b presents the scattering points estimated by 2-D ASSA and Figure 9c shows the full-band data image extrapolated by these extracted pole pairs. As can be seen, scattering centers of the base-edge and the body groove are distinguished from the extracted scattering points. The full-band data image extrapolated by 2-D ASSA is in keeping with the main scattering centers distributed in the original full-band data image (Figure 7a). It is clear from these figures that 2-D ASSA results in higher-resolution images than traditional 2-D FFT imagery.

Measured ISAR Signatures
In the third example, the measured ISAR signatures are acquired in an indoor test range where a mocked aircraft model is used. The data was collected using S-band radar with center frequency 3 GHz and 1.5 GHz bandwidth. The range of aspect angle is from −4 • to 4 • with an interval of 0.1 • . Figure 10a shows the ISAR image of the aircraft model generated by 2-D FFT. As it can be seen, the scattering distribution of this model is more complex than the previous model, i.e., many scattering points are densely distributed in the fuselage. The estimation parameter for 2-D ASSA is set to be L = N2; For 2-D ESPRIT, (P, Q) = (M2, N2) and β = 0.8.  Location estimation of main scattering centers for measured data are shown in Figure 10b-e. As we can see, both 2-D ESPRIT and 2-D ASSA could estimate the right positions of strong scattering centers. The scattering points estimated by 2-D ASSA are more densely distributed around the strong scattering centers than 2-D ESPRIT's for different pole-pairing strategies. A few relatively weak scattering points, i.e., the red point in rectangle box, are wrongly estimated by 2-D ESPRIT when the number of scattering points is set to be 115, and more wrongly estimated pole pairs can be seen in K > 115 which are not displayed in Figure 10. Figure 11a-h display a set of reconstructed 2-D images by the estimated pole pairs. Please note that the number of scattering points K in 2-D ESPRIT should be equal to the model order whereas K in 2-D ASSA could be larger than the model order K 1 (the number of non-repeated poles in the down-range dimension) when there are one-to-multiple pairing poles. That is the reason K could be set to 607 which is already exceed the limitation of model order in [17]. As we can see from the figures, all those two algorithms can reconstruct those strong scattering centers of target which are in consistent with the estimated poles in Figure 10a-d. As displayed in Figure 11f-h, more and more relatively weak scattering centers have been reconstructed by 2-D ASSA with the increasing number of K and the reconstruction result with compression ratio of 91.04% is extremely similar to the original data in Figure 10a.
Evaluation of the reconstructed results by 2-D ASSA are listed in Table 2. Compression ratio (CR) ε CR and image similarity degree (ISD) are defined in (46) and (47). The RRE δ RRE is used for evaluating the reconstructed result in frequency-domain, whereas the ISD γ ISD is for the result in image-domain. From the table, we can see that the numbers of scattering points in 2-D ASSA are negatively correlated with RRE (or ISD). In contrast, 2-D ESPRIT performs higher RRE than 2-D ASSA when the number of scattering centers remains the same. Please note that part of evaluation results by 2-D ESPRIT are abnormal because 2-D ESPRIT cannot do the right reconstruction when K > 307. The possible reason is that for 2-D ESPRIT, all amplitudes factors a 1 a 2 . . . a K (defined in (3)) are estimated together by using least square technique after all the pole pairs are confirmed. Thus, numerous incorrect pairing poles estimated by 2-D ESPRIT will lead to entirely wrong reconstruction. In comparison, 2-D ASSA performs better robustness for the measurement data.
where n P = 3 represents parameter number of each paired pole (a k , s k , p k ), K is the number of scattering points, M × N denotes size of the original signatures. I recon , I real represent the reconstructed and original 2-D images-domain data.  Table 3. All the results are carried out by Matlab (R2016b) with same hardware platform: Intel Core i7-6900k 3.2 GHz and 128 G memory. Please note that the running time of 2-D ESPRIT does not contain the estimation of model order. According to the table, 2-D ASSA performs higher efficiency than 2-D ESPRIT in processing the same data with the following differences.
(1) The size of the augmented Hankel matrix used by 2-D ASSA is approximate 1/(M/4) of the block-Hankel matrix used by 2-D ESPRIT when L = Q = N2 and P =M2. Moreover, the data size of the Hankel matrix could be further reduced by using the operation (22)-(24); (2) No extra pairing scheme is required by using multiple-range search strategy which can finish the pairing process once the pole estimation of the aspect dimension is confirmed. Moreover, calculation of the block-Hankel matrix for the aspect dimension in 2-D ESPRIT is simplified as calculating multiple small Hankel matrices with size M2 × (M − M2 + 1).

Conclusions
In this paper, we have presented a two-dimensional augmented state-space approach for sparse representation of 2-D wideband radar signatures collected on man-made targets. To do this, a two-step procedure, i.e., an augmented state-space approach followed by multiple-range search strategy, is proposed to estimate the complex amplitudes and poles in down-range and aspect dimensions. In general, there are mainly two contributions provided in this paper. First, we develop a computationally efficient approach by adopting several time-saving operations, whereas the pole-estimation accuracy is still at the same level; second, the proposed approach can apparently alleviate the pole-pairing problem by using the multiple-range search strategy.
Numerical as well as measured ISAR data are processed to validate the proposed approach. Experimental results demonstrate that 2-D ASSA is robust and accurate in pole-paring for different SNRs, and is applicable for sparse representation of 2-D wideband radar signatures collected on man-made targets with low computational cost.
Future works are considered to be as follows. On the one hand, physical meanings of the extracted parameters of scattering centers by 2-D ASSA might be further studied to the possible applications in automatic target recognition (ATR). On the other hand, the extension of the proposed approach to 3-D radar signatures obtained from different azimuth and elevation angles would also be an interesting study.
The coordinate form of the normalized eigenvalue vector is: where Q n denotes the n-th points. Based on information theory and principal factor analysis [28], the separable signals and noises in H often have the following hypothesis.
where K 1 is the number of separable signals which also denotes the model order. From (A5), it can be deduced that the point Q K 1 belonging to [Q 1 , Q 2 , . . . , Q n ] has the longest distance to line Q 1 D which has the following condition.
where i = 1, 2, . . . , K 1 , k → Q i Q K 1 +1 denotes the slope of line → Q i Q K 1 +1 , D is a point which satisfies the slope condition (A6).
Thus, Q K 1 can be confirmed by searching the longest distance from point Q i to line → Q 1 D.
where i = 1, 2, . . . , n; g (i) denotes the distance from point Q i to → Q 1 D.

Appendix B
From the observability matrix Ω given by (20), controllability matrix Γ given by (21) and the matrices Ω r f , Ω rl , Γ c f ,Γ cl in (27), (28), (29) and (30), it is not difficult to deduce that the augmented open-loop matrix P satisfies the following matrix equations.
Then the augmented open-loop matrix P can be obtained by least squares.
Moreover, here are two ways to compute the corresponding constant matrices S and P. The first way is for matrix P computed by (A14), the corresponding matrix S is defined as the first M rows of the observability matrix Ω. In this appendix, we provide two ways to estimate the parameter matrices (S, P, D). For consideration of the robustness to noise, (A23), (A14) (or (A15)) and (A19) are preferred. It is also worth mentioning that 2-D ASSA has two ways, which include the matrix form (S, P, D) and the pole form (a k , s k , p k ), to reconstruct the original signatures Y. The matrix form could provide more accurate reconstruction result than the pole form but with much more parameters.