1 Introduction

The ISAR imaging system transmits wideband signals to improve the range resolution, and improves the azimuth resolution by the coherent accumulation of the target echo. In practice, there is various perturbation factors that make the data received in the azimuth direction missing^[1]. However, the instability of the radar system may cause radar echo randomly lost. At this point, the use of traditional FFT methods will lead to the target defocus and other issues, resulting in unidentifiable imaging results.

At present, the main method for ISAR imaging of missing data is the Compression Sensing (CS)^[2-4], such as an iterative or non-iterative Orthogonal Matching Pursuit(OMP)^[5] algorithm. As long as the signal satisfies a certain sparsity level, the time-domain can be sampled at a sampling frequency lower than the Nyquist theorem, so that the signal can be approximately reconstructed with high efficiency. However, the algorithm needs to predict the sparsity level of the signal, and when the percentage of missing is large, the recovery performance will be severely degraded.

Since the core issue is using the missing data to estimate the spectrum of the signal in each range unit, the spectral estimation algorithms can be used for ISAR imaging. In the past few years, modern spectral estimation algorithms have been widely used for ISAR super-resolution imaging^[6-10]. As a branch of modern spectral estimation method, adaptive spectral estimation has been greatly developed recently. Among them are the Capon^[11] and the Amplitude and Phase Estimation of a Sinusoid (APES)^[12-13], which both perform spectrum estimation by adaptively constructing a filter bank.

Recently, the non-parametric weighted least-squares Iterative Adaptive Approach (IAA)^[14-15] has further improved the accuracy of the adaptive spectral estimation method. IAA constructs the signal covariance matrix using the spectral estimation results of the previous iteration and takes its inverse matrix into the weighted least-squares optimization problem. Because this method is suitable for non-uniform sampling of data, IAA is applied to recover missing data^[16]. The method is not affected by the signal parameters, the signal can be recovered in the case of a large proportion of missing data, accurate imaging with small rotational angle^[10]. In this paper, we use IAA for missing data ISAR imaging. Using IAA for the azimuth compression of missing data can effectively solve the problem of spectrum leakage. The imaging results of simulated and measured data illustrate the effectiveness of IAA.

Structure of this paper: In Section 2, the ISAR imaging principle and the random missing data model is presented; Section 3 presents IAA; In Section 4, the ISAR missing-data imaging results based on IAA are presented, and through the comparison with OMP illustrates the superiority of IAA; Section 5 is the conclusion of the paper.

2 ISAR Imaging Principle and ISAR Random Missing Data Model

In order to improve the range resolution, the radar usually emits wideband signals and the most widely used is the linear frequency modulated (LFM) signal. We set the pulse emitting time as t_m=mT, m=0, 1, 2, ... is the slow time and the fast time is $\hat{t} $. The transmitted signal can be expressed as:

$ s\left( {\hat t,{t_m}} \right) = rec{t_{{T_p}}}\left( {\hat t} \right)\exp \left[ {{\rm{j}}2{\rm{ \mathsf{ π} }}\left( {{f_c}t + \frac{1}{2}k{{\hat t}^2}} \right)} \right] $

(1)

where k is the frequency modulation slope, f_c is center frequency, rect_Tp(u) is a rectangular window of width T_p. Let the distance between the scattering points and the radar is R_i, the echo is:

$ \begin{array}{*{20}{c}} {{s_r}\left( {\hat t,{t_m}} \right) = \sum\limits_{i = 1}^K {{A_i}} \cdot rec{t_{{T_p}}}\left( {\hat t} \right)\exp \left[ {{\rm{j}}2{\rm{ \mathsf{ π} }}\left( {\hat t - } \right.} \right.}\\ {\left. {\left. {\frac{{2{R_i}\left( {{t_m}} \right)}}{c}} \right) + \frac{1}{2}k{{\left( {\hat t - \frac{{2{R_i}\left( {{t_m}} \right)}}{c}} \right)}^2}} \right]} \end{array} $

(2)

The signal after pulse compression by STRETCH processing in the fast time is:

$ \begin{array}{*{20}{c}} {{S_r}\left( {\hat f,{t_m}} \right) = \sum\limits_{i = 1}^K {{A_i}} \cdot {\rm{sinc}}\left[ {{T_p}\left( {\hat f - \frac{{2k\mathit{\Delta }{R_i}\left( {{t_m}} \right)}}{c}} \right)} \right] \cdot }\\ {\exp \left[ {\frac{{ - {\rm{j}}4{\rm{ \mathsf{ π} }}\mathit{\Delta }{R_i}\left( {{t_m}} \right)}}{\lambda }} \right]} \end{array} $

(3)

where λ is the wavelength of the transmitted signal, ΔR_i(t_m) is the distance between each scattering point and the reference point. After the motion compensation, the turntable imaging geometry is shown in Fig. 1.

In Fig. 1, ω is the angular velocity of rotation, the relative distance can be expressed as

$ \mathit{\Delta }{R_i}\left( {{t_m}} \right) = {x_i}\sin \left( {\omega {t_m}} \right) + {y_i}\cos \left( {\omega {t_m}} \right) $

If t_m is very small, then we can obtain

$ \sin \left( {\omega {t_m}} \right) \approx \omega {t_m},\cos \left( {\omega {t_m}} \right) \approx 1 $

Therefore, at a slow time, the signal within each range cell can be represented as the sum of the single frequency signals:

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

(4)

If the signal s(t_m) in each range cell is converted to the Doppler domain, the distribution of the scattering points in the azimuth can be obtained. In practice, there is various perturbation factors that make the data received in the azimuth direction appear missing. When part of the echo is missing, the signal at slow time can be expressed as^[2]:

$ {s_g}\left( {{m_g}} \right) = \sum\limits_{i = 1}^K {{\sigma _i}\exp \left( {\frac{{ - {\rm{j}}4{\rm{ \mathsf{ π} }}{x_i}\omega {m_g}}}{\lambda }} \right)} $

(5)

where {m_g}_g=1^G are the receipt times of valid data. At this point, if we still use the fast Fourier transform to analyze the Doppler spectrum of the signal, there will be serious spectrum leakage.

3 Weighted Least-Squares Iterative Adaptive Approach

As aforementioned analysis, the missing data imaging problem can be translated into spectral estimation of non-uniformly sampled signal s_g. Let K denote the number of sampling points in the frequency domain. Usually, K is more than the number of the points in the time domain. Let:

$ {\mathit{\boldsymbol{A}}_g} = \left[ {{\mathit{\boldsymbol{a}}_g}\left( {{\omega _1}} \right), \cdots ,{\mathit{\boldsymbol{a}}_g}\left( {{\omega _K}} \right)} \right] $

(6)

$ {\mathit{\boldsymbol{a}}_g}\left( {{\omega _k}} \right) = \left[ {\exp \left( {{\rm{j}}{\omega _k}{m_1}} \right), \cdots ,\exp \left( {{\rm{j}}{\omega _k}{m_G}} \right)} \right] $

(7)

where {a_g(ω_k)}_k=1^K is the Fourier vectors corresponding to s_g at {ω_k}_k=1^K. The valid data s_g can be modeled as:

$ {\mathit{\boldsymbol{s}}_g} = {\mathit{\boldsymbol{A}}_g}\mathit{\boldsymbol{\alpha }} + e $

(8)

where α= [α(ω₁), ..., α(ω_K)]^T is the vector of complex spectral amplitudes at frequency {ω_k}, e is Gaussian complex noise, (·)^T denotes the transpose. Assuming that the signal is not correlated with noise, the noise and interference covariance matrix for the valid data can be defined as:

$ \mathit{\boldsymbol{Q}} = {\mathit{\boldsymbol{R}}_g} - {\left| {\alpha \left( {{\omega _k}} \right)} \right|^2}{\mathit{\boldsymbol{a}}_g}\left( {{\omega _k}} \right)\mathit{\boldsymbol{a}}_g^ * \left( {{\omega _k}} \right) $

(9)

$ \mathit{\boldsymbol{R}} = \sum\limits_{k = 1}^K {{{\left| {\alpha \left( {{\omega _k}} \right)} \right|}^2}{\mathit{\boldsymbol{a}}_g}\left( {{\omega _k}} \right)\mathit{\boldsymbol{a}}_g^ * \left( {{\omega _k}} \right)} = \mathit{\boldsymbol{AP}}{\mathit{\boldsymbol{A}}^ * } $

(10)

where (·)^* denotes conjugate transpose, R is the covariance matrix of the valid data, P is a K × K diagonal matrix of diagonal elements {|α(ω_k)|²}_k=1^K of signal power at frequency {ω_k}. The weighted least square (WLS) estimation is used to estimate the complex amplitudes {α(ω_k)}_k=1^K, and the following optimization problem is obtained:

$ \mathop {\min }\limits_{\alpha \left( {{\omega _k}} \right)} {\left( {{\mathit{\boldsymbol{s}}_g} - \alpha \left( {{\omega _k}} \right){\mathit{\boldsymbol{a}}_g}\left( {{\omega _k}} \right)} \right)^ * }{\mathit{\boldsymbol{Q}}^{ - 1}}\left( {{\mathit{\boldsymbol{s}}_g} - \alpha \left( {{\omega _k}} \right){\mathit{\boldsymbol{a}}_g}\left( {{\omega _k}} \right)} \right) $

(11)

$ \hat \alpha \left( {{\omega _k}} \right) = \frac{{{\mathit{\boldsymbol{a}}^ * }\left( {{\omega _k}} \right){\mathit{\boldsymbol{Q}}^{ - 1}}{\mathit{\boldsymbol{s}}_g}}}{{{\mathit{\boldsymbol{a}}^ * }\left( {{\omega _k}} \right){\mathit{\boldsymbol{Q}}^{ - 1}}\mathit{\boldsymbol{a}}\left( {{\omega _k}} \right)}} $

(12)

The detailed implementation of it can be found in Ref. [16]. Using matrix inversion lemma, we can simplify Eq. (12) to:

$ \hat \alpha \left( {{\omega _k}} \right) = \frac{{{\mathit{\boldsymbol{a}}^ * }\left( {{\omega _k}} \right){\mathit{\boldsymbol{R}}^{ - 1}}{\mathit{\boldsymbol{s}}_g}}}{{{\mathit{\boldsymbol{a}}^ * }\left( {{\omega _k}} \right){\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{a}}\left( {{\omega _k}} \right)}} $

(13)

Using the iterative method, first let P be the unit matrix as the initial value, and then calculate R by Eq.(10) and $ \left\{ \hat{\alpha }\left( {{\omega }_{k}} \right) \right\}_{k=1}^{K}$ by Eq.(13) until it converges. Usually, it is converged after 15 iterations. The method can be summarized through the following steps:

Step 1  Find the receipt times {m_g}_g=1^G of valid data S_g, construct the Fourier matrix A in the form of Eq.(6).

Step 2  Set P as the unit matrix, calculate the initial value of the covariance matrix R by Eq.(10).

Step 3  Calculate $ \left\{ \hat{\alpha }\left( {{\omega }_{k}} \right) \right\} $ by Eq.(13).

Step 4  Calculate R by Eq. (10), return to the Step 3, until the iteration converges.

The computational cost of computing a one-dimensional signal using the proposed algorithm is N(2KG²+KG+G³)^[17]. Where N is the number of iterations, G is the valid data length. The running time of the IAA algorithm is dominated by calculating the inverse of the covariance matrix. When valid data is large in length, iterations take a long time.

We simulate a multi-component complex sinusoidal signal, and use it to verify the missing data iterative adaptive method. The amplitude and frequency of these four complex sinusoidal signals are (0.5 Hz, 0.25 Hz), (1.0 Hz, 0.39 Hz), (1.5 Hz, 0.67 Hz), (1.5 Hz, 0.77 Hz). The number of sampling points is 100, and the complex Gaussian white noise with the mean value of 0 and variance 0.01 is added to the signal. To verify the ability of the IAA method to estimate line spectra with missing data, 50% of the data is randomly lost. Fig. 2(a) is the result of missing data spectrum estimated by FFT, it can be seen that the directly use of FFT will produce serious spectrum leakage. Fig. 2(b) is the result of IAA, it can be seen that the IAA algorithm can accurately estimate the spectrum of signal in the case of a large proportion of missing.

4 ISAR Missing-Data Imaging Based on Iterative Adaptive Approach 4.1 Simulated Data Imaging

To illustrate the superiority of the IAA algorithm in ISAR imaging, we used simulated MIG-25 data provided by the Naval Research Laboratory-ABD (NRL) for the missing data imaging experiments. The step frequency radar works at 9 GHz with a bandwidth of 512 MHz and the Pulse repetition frequency is 15 kHz. A total of 512 pulses, 64 samples are collected from each pulse. This experiment only takes the first 64 echoes, and randomly lost part of the echo. The proportions of missing data are 25%, 50% and 75%. Fig. 3 (a), (d) and (g) is the results of missing data imaging by FFT with 75%, 50% and 25% valid data, as expected, the distortion of the image becomes more and more serious as the valid data decreases, so we cannot obtain the information of the target; Fig. 3(b), (e) and (h) show the results of missing data imaging by IAA with 75%, 50% and 25% valid data, it can be seen that the IAA method can be used to get a clear image even in the absence of a large amount of echo, the image is enough to discern the shape of the target. Fig. 3(c), (f) and (i) are the results of missing data imaging by OMP with 75%, 50% and 25% valid data, it shows that OMP can also use missing data imaging, but the reduction effect is slightly worse than IAA.

4.2 Measured Data Imaging

To prove the validity of IAA method in practice, we used the real data from the Yak-42 plane to perform the missing data imaging experiment. The radar transmits a LFM signal at a carrier frequency of 5520 MHz with a bandwidth of 400 MHz. The pulse repetition frequency is 400 Hz and the pulse width is 25.6 μs. As with the previous data, we randomly discard some of the echoes. The proportions of missing data are 25%, 50% and 75%. Fig. 4(a), (d) and (g) are the results of missing data imaging by FFT with 75%, 50% and 25% valid data; Fig. 4(b), (e) and (h) show the results of missing data imaging by IAA with 75%, 50% and 25% valid data; Fig. 4(c), (f) and (i) is the results of missing data imaging by OMP with 75%, 50% and 25% valid data. It can be seen that both IAA and OMP can perform missing data imaging. Due to the sparseness of OMP algorithm, the resolution of the results is higher than that of IAA. However, the IAA has a clearer graphical outline, so the recovery capability is better.

4.3 Imaging Effect Analysis

The results of three methods of missing data imaging were quantitatively analyzed by correlation coefficients. The correlation coefficient of 2d image is defined as:

$ r = \frac{{\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {\left( {{a_{i,j}} - \bar a} \right)\left( {{b_{i,j}} - \bar b} \right)} } }}{{\sqrt {\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {{{\left( {{a_{i,j}} - \bar a} \right)}^2}} } } \cdot \sqrt {\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {{{\left( {{b_{i,j}} - \bar b} \right)}^2}} } } }} $

(14)

where

$ \bar a = \frac{1}{{NM}}\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {{a_{i,j}}} } ,\bar b = \frac{1}{{NM}}\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {{b_{i,j}}} } $

where a_{i, j} and b_{i, j} are the grayscale of each point in the corresponding image. The correlation coefficient indicates the degree of correlation between two images. The larger the correlation coefficient is, the closer the two images are. Conversely, the smaller the correlation coefficient is, the greater the difference between the two images is. Therefore, in the case of using the same algorithm, the correlation coefficient between the missing data image and the corresponding complete data image indicates the ability of the algorithm for the missing data imaging. Here we use the simulated MIG-25 data in Section 4.1 to increase the data loss rate from 0 to 75%. For each data loss rate, the sampling loss time is randomly generated. 100 Monte Carlo experiments were performed to calculate the average correlation coefficient. The experimental results is shown in Fig. 5.

The three curves in Fig. 5 show that the correlation coefficients between the original image and the restored image change with the missing data percentage when using FFT, IAA and OMP. When there is no missing data, the correlation coefficients of all three algorithms are 1, indicating that the two images are exactly the same. As the percentage of missing data increases, the correlation coefficient of the three algorithms starts to decline, indicating that images show different degrees of distortion. Overall, both IAA and OMP are better than FFT. When missing data is less than 40%, IAA's recovery is slightly better with OMP. When missing data is more than 40%, OMP's resilience declines rapidly because of the increased signal uncertainty over the transform domain as the missing data increases, and the signal will not satisfy the sparsity level. However, IAA's performance is relatively stable. Similar results are obtained for the measured data in Section 4.2.

5 Conclusions

ISAR data loss may occur due to unstable radar systems and azimuth interference. The use of traditional FFT processing methods can lead to azimuth spectrum leakage. In this paper, an ISAR missing data imaging method based on IAA is proposed. The algorithm applies the weighted least squares estimation to obtain the spectrum of the non-uniformly sampled signal. The results of the simulated MIG-25 and the Yak-42 show that the algorithm can obtain a clear target image under the condition of random missing data.