1 Introduction

ISAR imaging system obtains the high range resolution by transmitting wideband signals and the azimuth resolution by the coherent accumulation of the target echo to achieve a Range-Doppler image. Due to the fact that the actual system may have a limited number of echoes or some of the echoes are subject to strong interference, the resolution of the traditional FFT-based Range-Doppler algorithm may not be able to meet the requirements. When the coherent accumulation angle is limited, it is necessary to use the super-resolution method for ISAR imaging^[1-2].

At present, most common ISAR super-resolution imaging methods are mainly classified into the nonparametric spectral estimation^[3-5] and the Compressive Sensing (CS)^[6-8]. The nonparametric spectrum estimation has a wide range of applications because no parameter needs to be set, while its resolution improvement is limited and the calculation is cumbersome. The CS algorithm exploits the sparseness of the echo at the scattering point, which makes a breakthrough to the Nyquist sampling theorem and can obtain a higher resolution. However, the use of the algorithm needs to predict the sparsity level of the signal and may lose part of the information.

In this paper, a nonparametric sparse spectrum estimation algorithm for ISAR imaging based on SPICE^[9-11] is proposed. The algorithm uses the sparseness of the signal without setting any parameters, which can effectively improve the resolution of ISAR imaging. In order to reduce the amount of calculation, this paper also introduces a fast implementation of the SPICE algorithm based on the G-S factorization^[12-13]. The ISAR imaging of simulated and measured data can illustrate the effectiveness of the SPICE algorithm.

The structure of this paper is as follows. The ISAR imaging principle is introduced in Section 2. Section 3 presents the principle of SPICE and the fast implementation of the SPICE algorithm is described in Section 4. The imaging results are discussed in Section 5 with the conclusions drawn in Section 6.

2 Signal Model

The general situation of ISAR imaging is that the radar does not move while the target moves. To simplify the model, as shown in Fig. 1, we assumed that the translational component of the target motion has been compensated, leaving only the rotational component.

It is assumed that during the observation process, the target turns a small angle against the reference point. At a slow time, the signal in each range cell can be expressed as

$ \begin{array}{*{20}{c}} {y\left( {{t_m}} \right) = \sum\limits_{i = 1}^K {{\sigma _i}} \exp \{ \frac{{ - {\rm{j}}4\pi }}{\lambda }\left[ {{x_i}\sin \left( {\omega {t_m}} \right) + } \right.}\\ {{y_i}\left( {1 - \cos \left( {\omega {t_m}} \right)} \right)]\} } \end{array} $

(1)

where λ denotes the wavelength, ω is the angular velocity of rotation, σ_i denotes the magnitude of the scattering point, and (x_i, y_i) denotes the position of the scattering point. If t_m is small, then sin(ωt_m)≈ω t_m and cos(ω t_m)≈1. Therefore, the signal within each range cell can be represented as single frequency signals as follows:

$ y\left( {{t_m}} \right) = \sum\limits_{i = 1}^K {{\sigma _i}} \exp \left( {\frac{{ - {\rm{j}}4\pi {x_i}\omega {t_m}}}{\lambda }} \right) $

(2)

In practice, there are only a few effective scattering points in a range cell, so the SPICE algorithm can be used to improve the resolution.

3 SPICE Algorithm

Let K denote the number of sampling points in the frequency domain. If y =[y(t₁), …, y(t_N)]^T denotes a time domain signal of length N, a_k=[exp(j ω_k t₁), …, exp(j ω_k t_N)]^T denotes a Fourier vector at ω_k, then we can rewrite formula (2) in the sparse form as follows:

$ \boldsymbol{y}=\sum\limits_{k=1}^{K} s_{k} \boldsymbol{a}_{k}+{\varepsilon}=\left[\boldsymbol{a}_{1}, \cdots, \boldsymbol{a}_{K}, \boldsymbol{I}\right]\left[\begin{array}{l}{\boldsymbol{s}} \\ {\boldsymbol{\varepsilon}}\end{array}\right]=\boldsymbol{A} \boldsymbol{\beta} $

(3)

where

$ \begin{array}{l}{\boldsymbol{s}=\left[s_{1}, \cdots, s_{k}\right]^{\mathrm{T}}, \boldsymbol{\beta}=\left[\boldsymbol{s}^{\mathrm{T}}, \boldsymbol{\varepsilon}^{\mathrm{T}}\right]^{\mathrm{T}}} \\ {\boldsymbol{A}=\left[\boldsymbol{a}_{1}, \cdots, \boldsymbol{a}_{K}, \boldsymbol{I}\right]=\left[\boldsymbol{a}_{1}, \cdots, \boldsymbol{a}_{K+N}\right]}\end{array} $

(4)

where ε denotes noise term, {s_k}_k=1^K denotes the corresponding amplitude at ω_k. Since y is a sparse signal, most of the elements in {s_k}_k=1^K are zero. If the noise covariance matrix is expressed as diag{σ₁, σ₂, …, σ_N}, the covariance matrix of the signal can be expressed as

$ \boldsymbol{R}= E\left(\boldsymbol{y} \boldsymbol{y}^{\mathrm{H}}\right)=\sum\limits_{k=1}^{K}\left|s_{k}\right|^{2} \boldsymbol{a}_{k} \boldsymbol{a}_{k}^{\mathrm{H}}+\\ \operatorname{diag}\left\{\sigma_{1}, \sigma_{2}, \right., \cdots, \sigma_{N} \} \triangleq \boldsymbol{A} \boldsymbol{P} \boldsymbol{A}^{\mathrm{H}} $

(5)

where H represents the conjugate transpose and

$ \begin{aligned} \boldsymbol{P}=& \operatorname{diag}\left\{\left|s_{1}\right|^{2}, \cdots, \left|s_{K}\right|^{2}, \sigma_{1}, \cdots, \sigma_{N}\right\}=\\ & \operatorname{diag}\left\{p_{1}, \cdots, p_{N+K}\right\} \end{aligned} $

(6)

Using the covariance fitting criterion, the following optimization problem is obtained

$ f=\left\|\boldsymbol{R}^{-1 / 2}\left(\boldsymbol{y} \boldsymbol{y}^{\mathrm{H}}-\boldsymbol{R}\right)\right\|^{2} $

(7)

where ‖·‖ represents the Frobenius norm of the matrix, and R^-1/2 is the square root of the Hermitian positive definite matrix R^-1. After a calculation^[4], it can be deduced that the minimization of f is the same as that of the following problem

$ \mathop {\min }\limits_p {\mathit{\boldsymbol{y}}^{\rm{H}}}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{y}} + \sum\limits_{k = 1}^{K + N} {w_k^2} {p_k};{w_k} = \frac{{\left\| {{\mathit{\boldsymbol{a}}_k}} \right\|}}{{\left\| \mathit{\boldsymbol{y}} \right\|}} $

(8)

For the convenience of calculation, Eq. (8) is changed to the following constraint minimization problem

$ \begin{aligned} \min\limits_{p, \beta} \boldsymbol{\beta}^{\mathrm{H}} \boldsymbol{P}^{-1} \boldsymbol{\beta} &+\sum\limits_{k=1}^{K+N} w_{k}^{2} p_{k} \\ \text { s.t } \boldsymbol{A} \boldsymbol{\beta} &=\boldsymbol{y} \end{aligned} $

(9)

During the iteration, we can fix P to calculate the update of β, and the results can be obtained from Ref. [4] as follows:

$ \boldsymbol{\beta}=\boldsymbol{P} \boldsymbol{A}^{\mathrm{H}} \boldsymbol{R}^{-1} \boldsymbol{y} $

(10)

After completing the update of β, P needs to be calculated by using formula (9). Since the scattering points of the target are independent from each other, the minimum update of each frequency point can be calculated separately. So Eq.(9) can be reduced to the following form

$ \min\limits_{p \geqslant 0} \frac{|\beta|^{2}}{p}+w^{2} p $

(11)

The optimal solution is

$ p_{k}=\frac{\left|\beta_{k}\right|^{2}}{w_{k}}, k=1, 2, \cdots, K+N $

(12)

It can be seen from the above analysis that, through iterative operations, the values of {p_k}_k=1^K+N and β are updated alternately so that the algorithm runs stably. According to Ref. [9], for any initial value {p_k≥0}, the results of the iteration will converge to the optimal solution. However, the initialization of p_k=|a_k^* y |²/‖a_k‖⁴(k=1, 2, …, K+N) by the single-frequency least-squares (SFLS) method can accelerate the convergence speed.

4 Fast Implementation of SPICE

The iteration of SPICE needs the calculation of the inverse of the covariance matrix R^-1 in Eq.(10), and the computational cost of computing R^-1 is O(N³). Since the running time of the SPICE algorithm is dominated by calculating the inverse of the covariance matrix, a method is introduced in this section for the fast calculation of R^-1 in the SPICE algorithm. R is rewritten as

$ \boldsymbol{R}_{N}=\sum\limits_{k=1}^{K}\left|s_{k}\right|^{2} \boldsymbol{a}_{k} \boldsymbol{a}_{k}^{*}+\operatorname{diag}\left\{\sigma_{1}, \sigma_{2}, \cdots, \sigma_{N}\right\}=\\ \left[\begin{array}{cccc}{r_{0}+\sigma_{1}} & {r_{1}} & {\cdots} & {r_{N-1}} \\ {r_{1}^{*}} & {r_{0}+\sigma_{2}} & {\cdots} & {r_{N-2}} \\ {\vdots} & {\vdots} & {\ddots} & {\vdots} \\ {r_{N-1}^{*}} & {r_{N-2}^{*}} & {\cdots} & {r_{0}+\sigma_{N}}\end{array}\right] $

(13)

Suppose that the noise variance is constant throughout the slow time, i.e., σ₁=…=σ_N=σ. R_N becomes a Toeplitz and Hermitian matrix, so we can use the method in Refs. [6-7] to perform G-S factorization to reduce the amount of computation for inversion. First, R_N is expressed as

$ \boldsymbol{R}_{N}=\left[\begin{array}{cc}{r_{0}+\sigma} & {\boldsymbol{r}_{N-1}^{\mathrm{H}}} \\ {\boldsymbol{r}_{N-1}} & {\boldsymbol{R}_{N-1}}\end{array}\right]=\left[\begin{array}{cc}{\boldsymbol{R}_{N-1}} & {\hat{\boldsymbol{r}}_{N-1}^{*}} \\ {\hat{\boldsymbol{r}}_{N-1}^{\mathrm{T}}} & {r_{0}+\boldsymbol{\sigma}}\end{array}\right] $

(14)

where r_N-1=[r₁^*, r₂^*, …, r_N-1^*]^T and $ \hat{\boldsymbol{r}}_{N-1}$ denote the reverse order of r_N-1. Then calculate

$ \boldsymbol{w}_{N-1}=-\boldsymbol{R}_{M-1}^{-1} \boldsymbol{r}_{M-1} $

(15)

$ \alpha_{N-1}=r_{0}+\sigma-\boldsymbol{r}_{N-1}^{\mathrm{H}} \boldsymbol{R}_{N-1}^{-1} \boldsymbol{r}_{N-1} $

(16)

Let L_N(t, Z)=[t, Zt, …, Z^N-1t]. Z_N denotes an M × M strictly lower-triangular matrix with ones on the first sub-diagonal and zeros everywhere else. So we can use G-S factorization to write R^-1 as

$ \boldsymbol{R}_{N}^{-1}=L_{N}\left(\boldsymbol{t}_{1}, \boldsymbol{Z}_{N}\right) L_{N}^{\mathrm{H}}\left(\boldsymbol{t}_{1}, \boldsymbol{Z}_{N}\right)-L_{N}\left(\boldsymbol{t}_{2}, \boldsymbol{Z}_{N}\right) L_{W}^{\mathrm{H}}\left(\boldsymbol{t}_{2}, \boldsymbol{Z}_{N}\right) $

(17)

where

$ \begin{aligned} \boldsymbol{t}_{1} &=\frac{1}{\sqrt{\alpha_{M-1}}}\left[1, \boldsymbol{w}_{M-1}\right]^{\mathrm{T}} \\ \boldsymbol{t}_{2} &=\frac{1}{\sqrt{\alpha_{M-1}}}\left[0, \hat{\boldsymbol{w}}_{M-1}^{*}\right]^{\mathrm{T}} \end{aligned} $

(18)

From the above analysis, it can be seen that only w_N-1 and α_N-1 are needed to get R^-1, and the calculation of w_N-1 and α_N-1 by the Levinson-Durbin algorithm^{[12, 14]} requires only N² operations. In Eq.(17), both L_N(t₁, Z_N) and L_N(t₂, Z_N) are triangular Toeplitz matrices, so the FFT can be used to calculate R^-1y in Eq.(10)^[14]. Compared with the direct calculation of R^-1, the proposed method greatly reduces the amount of computation for each iteration.

The fast implementation ofthe SPICE algorithm can be summarized as follows:

Initialization:

Calculate {p_k}_k=1^K+N using SFLS and calculate w_k in Eq.(8).

The (i+1)th iteration:

1) Construct the covariance matrix R in Eq.(5) using FFT with O(Klog₂K) flops^[13];

2) Calculate w_N-1 and α_N-1 using the Levinson-Durbin (L-D)-type algorithm with O(N²) flops;

3) Calculate R^-1y using FFT with O(N²log₂N) flops, calculate β by Eq.(10), and then update {p_k}_k=1^K+N by Eq.(12) with O(N²K) flops;

4) Stop the iterative operation according to the predetermined threshold.

5 Experiments and Discussion

To verify the proposed method, we simulated a multi-component sinusoidal signal. The amplitude and frequency of the four components are (0.5, 1.46 Hz), (1.0, 2.34 Hz), (1.5, 3.26 Hz), and (1.5, 3.46 Hz). SPICE algorithm was used for spectral estimation of this signal. The number of sampling points is 128. Fig. 2 shows the relationship between the energy difference and the number of iterations for adjacent iterations of the SPICE algorithm.

In practical applications, we can assume that the iteration converges when the energy difference is less than one percent of the total energy and the number of iterations is usually less than ten. Fig. 3 (a) is the result of spectrum estimation by FFT, which indicates that if FFT is used directly, the resolution cannot meet the requirement due to the small number of sampling points. Fig. 3 (b) is the result of SPICE, which shows that the SPICE algorithm had lower sidelobes and is more accurate for the amplitude estimates than the FFT.

Fig. 4 shows the flowchart for ISAR imaging using the SPICE algorithm. Motion compensation was performed on the high resolution range image, and the azimuth compression of the signal in each range cell was then performed using the improved SPICE algorithm proposed in this paper. In this experiment, the FFT, Iterative Adaptive Approach (IAA)^{[3, 15]}, and the SPICE algorithm were compared.

The simulated MIG-25 data was used for the ISAR imaging experiment. The stepped-frequency radar works at 9 GHz and the pulse repetition frequency is 15 kHz. A total of 512 pulses, 64 samples were collected from each pulse. Figs. 5 (a)-(c) show the results of MIG-25 imaging by FFT, IAA, and the fast implementation of SPICE algorithm, respectively. It can be seen that both IAA and SPICE algorithm can improve the resolution of the image, while the SPICE algorithm works better than the IAA.

The entropy of the three images was calculated to compare their resolutions, whose results were 7.859 9, 7.272 5, and 6.549 6, respectively, which indicates that the ISAR image obtained by the SPICE algorithm has the best focus performance.

In order to verify the effectiveness of the SPICE in practice, the measured data from the Yak-42 plane was used. The radar transmits an LFM signal at a carrier frequency of 5 520 MHz with a bandwidth of 400 MHz. The pulse repetition frequency is 400 Hz. Here we only use 128 consecutive pulses. Fig. 6 (a) shows the result of Yak-42 imaging by FFT, which was blurry due to the short accumulation time. Fig. 6 (b) illustrates the result of IAA, which indicates that the IAA algorithm can improve the resolution of ISAR imaging. Fig. 6 (c) is the imaging result of the fast implementation of SPICE algorithm, which shows that the image quality was further improved.

The entropy of the three images was calculated, in which the image using FFT was 7.327 3, using IAA was 7.042 7, and using the fast implementation of SPICE algorithm was 6.613 1, among which the image resolution using the SPICE algorithm was the highest.

6 Conclusions

In this paper, a novel high-resolution imaging algorithm based on SPICE is proposed to solve the problem of short accumulation time of ISAR imaging. The basic principles of ISAR imaging and SPICE algorithm were briefly described, and a fast implementation of the SPICE was introduced. The amount of calculation was reduced by G-S factorization and FFT. Experiments on simulated MIG-25 data and measured Yak-42 data show that the proposed algorithm can significantly improve the resolution of ISAR imaging.