ISAR Imaging for Maneuvering Targets with Complex Motion Based on Generalized Radon-Fourier Transform and Gradient-Based Descent under Low SNR

Abstract

The existing inverse synthetic aperture radar (ISAR) imaging algorithms for ship targets with complex three-dimensional (3D) rotational motion are not applicable because of continuous change of image projection plane (IPP), especially under low signal-to-noise-ratio (SNR) condition. To overcome this obstacle, an efficient approach based on generalized Radon Fourier transform (GRFT) and gradient-based descent optimal is proposed in this paper. First, the geometry and signal model for nonstationary IPP of ship targets with complex 3-D rotational motion is established. Furthermore, the two-dimensional (2D) spatial-variant phase errors caused by complex 3-D rotational motion which can seriously blur the imaging performance are derived. Second, to improve the computational efficiency for 2-D spatial-variant phase errors compensation, the coarse motion parameters of ship targets are estimated via the GRFT method. In addition, using the gradient-based descent optimal method, the global optimum solution is iteratively estimated. Meanwhile, to solve the local extremum for cost surface obtained via conventional image entropy, the image entropy combined with subarray averaging is applied to accelerate the global optimal convergence. The main contributions of the proposed method are: (1) the geometry and signal model for ship targets with a complex 3-D rotational motion under nonstationary IPP are established; (2) the image entropy conjunct with subarray averaging operation is proposed to accelerate the global optimal convergence; (3) the proposed method can ensure the imaging accuracy even with high imaging efficiency thanks to the sole optimal solution generated by using the subarray averaging and image entropy. Several experiments using simulated and electromagnetic data are performed to validate the effectiveness of the proposed approach.

1. Introduction

Inverse synthetic aperture radar (ISAR) is an applicable technique to obtain high-resolution ISAR images for targets because the structure, size, and shape can be reconstructed using the echoes reflected from it under all-day and all-weather condition. Therefore, it can be widely utilized in military and civilian fields [1,2,3].

Compared with cooperative targets, e.g., on-orbit satellite and airplane, non-cooperative targets such as ship targets have complex three-dimensional (3D) rotational motions, e.g., roll, pitch and yaw, and translational motions. Typically, the translational motions can be accurately compensated via a standard compensation algorithm [4]. However, according to the analysis presented in [5,6], the components of 3D rotational motions are time-varying in amplitude and direction vectors. As a result, the image projection plane (IPP) of the targets with complex 3-D rotational motion presents nonstationary characteristics, which would violate the assumption that the IPP is fixed during coherent processing interval (CPI), and followed by the inapplicable of existing ISAR imaging approach for non-cooperative targets.

Several imaging algorithms for maneuvering targets with complex motion have been proposed in recent years [7,8,9,10,11,12,13,14,15,16,17,18,19,20,21]. They can be roughly divided into parametric-based and non-parametric methods. For the parametric-based method, the typical methods are modeling the signals of a specific range bin as quadratic phase signals or cubic phase signals [7,8,9,10,11,12,13]. By using parameters estimation methods, e.g., cubic phase function (CPF) [8,9], non-uniform sampled CPF [10], integrated generalized cubic phase function (IGCPF) [11], scaled Fourier transform (SCFT)-based algorithm [12], generalized decoupling technique (GDT)-based algorithm [13], et al., the coefficients of those phase signal can be accurately estimated. However, the operation for selecting a specific range bin to estimate the coefficients is computationally extensive and infeasible under low SNR conditions. For non-parametric-based methods, the typical methods are time-frequency distribution (TFD) [14] based or polynomial-phase transform (PPT) based [15], which include short-time Fourier transform (STFT) [16], continuous wavelet transform (CWT) [17], Wigner-Ville distribution (WVD) et al. Those methods can reduce the order of signals by using a nonlinear transform operation. Nevertheless, the cross-term interference will be generated while processing multicomponent chirps signals. Besides, the resolution is low, which would also affect the application in the real world. To improve the robustness and effectiveness for ISAR imaging of ship targets, the autofocus algorithm based on data-driven is proposed in [18]. The approach can be divided into two kinds, e.g., estimating and compensating phase errors from the image domain or from the signal domain. The essence of the first one is that the phase errors are modeled as three orders or higher orders polynomial. With the evaluation indicators in the image domain, the phase errors can be estimated and compensated via existing optimum approaches. In [19], the phase gradient autofocus (PGA) is presented. However, the number of iteration and the width of the data window are troublesome problems. The key of the second one is accurately extracting phase errors from the signal domain. In [20], the phase errors between consecutive two pulses are estimated. However, the accumulation of phase errors will inevitably occur while processing multiple pulses. Additionally, the ship target images can be reconstructed using at least two times GRFT [21], where the GRFT method is used for coarse and fine estimation for motion parameters, which is time-consuming because the estimation accuracy is determined by the search step of the motion parameters.

In addition, the methods based on optimum coherent processing interval (CPI) selection are developed for ISAR imaging of ship targets with complex 3-D motion. By coarsely reconstructing the images and extracting the features from which, the maneuvering severity of ship targets is determined and followed by the selection of optimum CPIs [22,23]. However, it suffers from serious efficiency problems. Further, the methods of analyzing Doppler frequency from prominent scatterers are proposed to estimate the characteristics of target motion [24,25]. However, these algorithms either need to detect strong scatterers or have high computational complexity, thus the application in practice has limitations.

To obtain well-focused ISAR images for ship targets with a complex 3-D rotational motion under low SNR, a ship ISAR imaging algorithm based on the GRFT and gradient-based descent optimal is proposed in this paper. Considering the nonstationary characteristic of IPP during CPI, the radar LOS is modeled as the function in terms of slow time, which can accurately describe the motion characteristic for ship targets with complex 3-D rotational motion. Meanwhile, the 2-D spatial-variant phase errors caused by the 3-D rotational motion are derived. Additionally, the GRFT-based is introduced to a rough estimate of the motion parameters. Furthermore, the accurate motion parameters are estimated using the gradient-based descent optimal method. Considering the local convergence of cost surface obtained using conventional image entropy, the approach of image entropy combined with subarray averaging operation is used to improve the convergence efficiency for the global optimal solution. Accordingly, the 2-D spatial-variant phase errors can be precisely estimated and followed by the well-focused ISAR images. Simulated data and electromagnetic data are utilized to verify the effectiveness of the proposed approach. Compared with the existing imaging algorithm for ship targets, the main contribution is as follows: (1) the signal model for ship target with nonstationary IPP is derived, which can accurately describe the motion characteristic of ship targets with complex 3-D rotational motion; (2) the GRFT combined with gradient-based optimal estimation is proposed to improve the processing efficiency for motion parameters estimation; (3) the image entropy based on subarray averaging operation is applied to accelerate the global convergence for the optimum solution.

The rest of this work is organized as follows. In Section 2, the geometric and signal model for ship target with complex 3-D rotational motion are introduced, and 2-D spatial-variant phase errors with nonstationary IPP are also provided. An efficient parameters estimation approach based on the GRFT method and gradient descent approach is proposed in Section 3, where GRFT is adopted to a rough estimate of the motion parameters, and the gradient-based optimal combined with subarray averaging operation is proposed to exactly estimate the motion parameters. At the same time, some considerations in practical application are presented in this part. The experimental results and corresponding analysis with simulated and electromagnetic data are described in Section 4, and some conclusions are summarized in Section 5.

2. ISAR Imaging for Ship Targets

In this section, the geometry model and three-dimensional (3D) rotational motion model are given in Figure 1, where the Cartesian coordinate $\left(X,Y,Z\right)$ is established in the target body, and origin $O$ is the rotation center, ${\eta }_p$ and ${\phi }_p$ denote the elevation angle and azimuth angle of the radar line-of-sight (LOS), respectively. ${\theta }_y$, ${\theta }_r$, and ${\theta }_p$, respectively, represent the angular motion yaw, roll, and pitch, which are rotating around $Z$, $X$, and $Y$ axes, respectively.

2.1. Signal Model for Ship Targets

Now, we suppose the linear-frequency-modulated (LFM) signals are transmitted in the radar system, and it can be written as

$$\mathrm{s}\left(t_r\right)=\mathrm{rect}\left(\frac{t_r}{T_p}\right)exp\left(j2\pi \left(f_{c}t+\frac{1}{2}{K}_r{t}_r^2\right)\right)$$

(1)

where $\mathrm{rect}\left(\mathrm{x}\right)=\{\begin{array}{c}1, \left|\mathrm{x}\right|\le \frac{1}{2}\\ 0, \left|\mathrm{x}\right|>\frac{1}{2}\end{array}$, $t_r,t_a,f_c,K_r,T_p$ denote the fast time, slow time, carrier frequency, frequency modulation rate, and pulse width, respectively, $t=t_r+m·T_p,m=0,1,2,\dots \left(M-1\right)$, stands for the full time, and $M$ is the total number of the received pulses. As shown in Figure 1, arbitrary scatterer P is located on the target body, whose coordinate is $\left(x_p,y_p,0\right)$.

The echoes of the scatterer P after demodulation are given by

$$\begin{array}{ll}\mathrm{s}\left(t_r,t_a\right)={\sigma }_p·& \mathrm{rect}\left(\frac{t_r-2R_p\left(t_a\right)/c}{T_p}\right)\mathrm{rect}\left(\frac{t_a}{T_a}\right)\\ & \times exp\left(j2\pi \left(-f_{c}2R_p\left(t_a\right)/c+\frac{1}{2}{K}_r{\left(t_r-2R_p\left(t_a\right)/c\right)}^2\right)\right)\end{array}$$

(2)

where ${\sigma }_p,c,T_a,R_p\left(t_a\right)$ denote the reflected coefficients, speed of light, coherent integration time, and instantaneous slant range from radar to scatterer P, respectively.

Conducting Fourier transform (FT) along with $t_r$ and range compression to (2), one obtains

$$S\left(f_r,t_a\right)={\sigma }_p\cdot \mathrm{w}\left(f_r\right)·\mathrm{rect}\left(\frac{t_a}{T_a}\right)exp\left(-j\frac{4\pi R_p\left(t_a\right)}{c}\left(f_r+f_c\right)\right)$$

(3)

where $\mathrm{w}\left(f_r\right)$ denotes the frequency window function.

Generally speaking, the instantaneous slant range $R_p\left(t_a\right)$ of the target can be decomposed into translational motions part $R_T\left(t_a\right)$ and rotational motions part $R_r\left(t_a\right)$, given by

$$R_p\left(t_a\right)=R_r\left(t_a\right)+R_T\left(t_a\right)$$

(4)

The translational motion part $R_T\left(t_a\right)$ should be compensated for because all of the scatterers in the target body share the same part that has no contribution to ISAR imaging. The rotational motions $R_r\left(t_a\right)$ can be calculated as [26]

$$R_r\left(t_a\right)={\left[rot\left(t_a\right)\cdot P\right]}^T\cdot {\mathit{i}}_{los}$$

(5)

where $P={\left[x_p,y_p,0\right]}^T$, $rot\left(t_a\right)$, and ${\mathit{i}}_{los}$, respectively, are the coordinate of the scatterer P, rotational matrix, and radar LOS, ${[·]}^T$ denotes the transposition. It should be noted that, from (5), the rotational motion is related to ${\mathit{i}}_{los}$ and $rot\left(t_a\right)$ during the CPI. Obviously, the structure of the targets is projected to the 2-D image plane, e.g., IPP. The range dimension is defined as the direction of LOS, while the cross-range dimension is defined as the cross-product of the radar LOS direction and the effective vector. Therefore, the definition of IPP is related to the radar LOS and effective rotational vector. In general, ${\mathit{i}}_{los}$ is a constant during the CPI if the motions of targets are moderate. However, when the targets are involved in complex 3-D rotational motion, ${\mathit{i}}_{los}$ are varied by the azimuth time, which will cause the change of IPP. Therefore, to accurately describe the motion characteristic of maneuvering targets, in this work, the ${\phi }_p$ and ${\eta }_p$ can be modeled as the function of slow time $t_a$, which form the unit vector for the direction of ${\mathit{i}}_{los}$, given by

$$i_{los}={\left[cos\left({\phi }_p\left(t_a\right)\right)cos\left({\eta }_p\left(t_a\right)\right),cos\left({\phi }_p\left(t_a\right)\right)sin\left({\eta }_p\left(t_a\right)\right),sin\left({\phi }_p\left(t_a\right)\right)\right]}^T$$

(6)

where

$${\phi }_p\left(t_a\right)={\phi }_0+\kappa t_a$$

(7)

$${\eta }_p\left(t_a\right)={\eta }_0+\gamma t_a$$

(8)

Further, based on second-order Taylor series expression, we have

$$cos\left(x\right)\approx 1-\frac{x}{2}$$

(9)

$$sin\left(x\right)\approx x$$

(10)

Thus, substituting (7)–(10) into (6),${\mathit{i}}_{los}$ is as ${\mathit{i}}_{los}={\left[i_1\left(t_a\right),i_2\left(t_a\right),i_3\left(t_a\right)\right]}^T$, where

$$i_1\left(t_a\right)=k_1+k_2\cdot t_a$$

(11)

$$i_2\left(t_a\right)=k_3+k_4\cdot t_a$$

(12)

$$i_3\left(t_a\right)=k_5+k_6\cdot t_a$$

(13)

where

$$k_1=1-\frac{1}{2}{\eta }_0-\frac{1}{2}{\phi }_0+\frac{1}{4}{\eta }_0{\phi }_0$$

(14)

$$k_2=-\frac{1}{2}\left(\gamma +\kappa -\frac{1}{2}{\phi }_0\gamma -\frac{1}{2}{\eta }_0\kappa \right)$$

(15)

$$k_3={\eta }_0-\frac{1}{2}{\eta }_0{\phi }_0$$

(16)

$$k_4=\gamma -\frac{1}{2}{\phi }_0\gamma -\frac{1}{2}{\eta }_0\kappa$$

(17)

$$k_5={\phi }_0$$

(18)

$$k_6=\kappa$$

(19)

In addition, the form of $rot\left(t_a\right)$ [27] in (5) can be written as

$$rot\left(t_a\right)=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]$$

(20)

where

$$\{\begin{array}{c}{a}_{11}=cos\left({\theta }_p\left(t_a\right)\right)cos\left({\theta }_y\left(t_a\right)\right)\\ a_{12}=-cos\left({\theta }_p\left(t_a\right)\right)sin\left({\theta }_y\left(t_a\right)\right)\\ a_{13}=sin\left({\theta }_p\left(t_a\right)\right)\\ a_{21}=sin\left({\theta }_r\left(t_a\right)\right)sin\left({\theta }_p\left(t_a\right)\right)cos\left({\theta }_y\left(t_a\right)\right)+cos\left({\theta }_r\left(t_a\right)\right)sin\left({\theta }_y\left(t_a\right)\right)\\ a_{22}=-sin\left({\theta }_r\left(t_a\right)\right)sin\left({\theta }_p\left(t_a\right)\right)sin\left({\theta }_y\left(t_a\right)\right)+cos\left({\theta }_r\left(t_a\right)\right)cos\left({\theta }_y\left(t_a\right)\right)\\ a_{23}=-sin\left({\theta }_r\left(t_a\right)\right)cos\left({\theta }_p\left(t_a\right)\right)\\ a_{31}=-cos\left({\theta }_r\left(t_a\right)\right)sin\left({\theta }_p\left(t_a\right)\right)cos\left({\theta }_y\left(t_a\right)\right)+sin\left({\theta }_r\left(t_a\right)\right)sin\left({\theta }_y\left(t_a\right)\right)\\ a_{32}=cos\left({\theta }_r\left(t_a\right)\right)sin\left({\theta }_p\left(t_a\right)\right)sin\left({\theta }_y\left(t_a\right)\right)+sin\left({\theta }_r\left(t_a\right)\right)cos\left({\theta }_y\left(t_a\right)\right)\\ a_{33}=cos\left({\theta }_r\left(t_a\right)\right)cos\left({\theta }_p\left(t_a\right)\right)\end{array}$$

(21)

where ${\theta }_r\left(t_a\right),{\theta }_y\left(t_a\right),{\theta }_p\left(t_a\right)$ can also be expressed as

$$\{\begin{array}{c}{\theta }_r\left(t_a\right)\approx {\omega }_r{t}_a+\frac{{\omega }_r^{\prime }}{2}\cdot t_a^2\\ {\theta }_p\left(t_a\right)\approx {\omega }_p{t}_a+\frac{{\omega }_p^{\prime }}{2}\cdot t_a^2\\ {\theta }_y\left(t_a\right)\approx {\omega }_y{t}_a+\frac{{\omega }_y^{\prime }}{2}\cdot t_a^2\end{array}$$

(22)

where ${\omega }_r,{\omega }_y,{\omega }_p,{\omega }_r^{\prime },{\omega }_r^{\prime },{\omega }_r^{\prime }$ represent the constant rotation velocity of roll, yaw, and pitch, the rotation acceleration of roll, yaw, and pitch, respectively.

Therefore, $rot\left(t_a\right)$ can be written as

$$rot\left(t_a\right)\approx \left[\begin{array}{ccc}1-\frac{\left({\omega }_y^2+{\omega }_p^2\right)}{2}{t}_a^2& -{\omega }_y{t}_a-\frac{{\omega }_y^{\text{'}}}{2}{t}_a^2& {\omega }_p{t}_a+\frac{{\omega }_p^{\text{'}}}{2}{t}_a^2\\ {\omega }_y{t}_a+\left({\omega }_r{\omega }_p+\frac{{\omega }_y^{\text{'}}}{2}\right)t_a^2& 1-\frac{\left({\omega }_y^2+{\omega }_r^2\right)}{2}{t}_a^2& -{\omega }_r{t}_a-\frac{{\omega }_r^{\text{'}}}{2}{t}_a^2\\ -{\omega }_p{t}_a+\left({\omega }_r{\omega }_y-\frac{{\omega }_p^{\text{'}}}{2}\right)t_a^2& {\omega }_r{t}_a+\left({\omega }_y{\omega }_p+\frac{{\omega }_r^{\text{'}}}{2}\right)t_a^2& 1-\frac{\left({\omega }_p^2+{\omega }_r^2\right)}{2}{t}_a^2\end{array}\right]$$

(23)

Therefore, $R_r\left(t_a\right)$ can be re-expressed as

$$R_r\left(t_a\right)=K_0+K_1\cdot t_a+K_2\cdot t_a^2$$

(24)

where

$$K_0=k_1{x}_p+k_3{y}_p$$

(25)

$$K_1=\left(k_2+k_3{\omega }_y-k_5{\omega }_p\right)x_p+\left(k_4-k_1{\omega }_y+k_5{\omega }_r\right)y_p$$

(26)

$$\begin{gathered} K_2=\frac{1}{2}\left(-k_1\left({\omega }_y^2+{\omega }_p^2\right)+k_3\left({\omega }_y^{\prime }+2{\omega }_r{\omega }_p\right)-k_5\left({\omega }_p^{\prime }-2{\omega }_r{\omega }_y\right)-2k_6{\omega }_p\right)x_p \\ +\frac{1}{2}\left(-k_1{\omega }_y^{\prime }-2k_2{\omega }_y-k_3\left({\omega }_y^2+{\omega }_r^2\right)+k_5\left({\omega }_r^{\prime }+2{\omega }_p{\omega }_y\right)+2k_6{\omega }_r\right)y_p \end{gathered}$$

(27)

According to the analysis above, $K_0$ and $K_1$ can be written as

$$\left[\begin{array}{c}{K}_0\\ K_1\end{array}\right]=\left[\begin{array}{cc}{k}_1& k_3\\ k_2+k_3{\omega }_y-k_5{\omega }_p& k_4-k_1{\omega }_y+k_5{\omega }_r\end{array}\right]\left[\begin{array}{c}{x}_p\\ y_p\end{array}\right]$$

(28)

Furthermore, $K_2$ can be expressed as the function of $K_0$ and $K_1$, given by

$$\begin{gathered} K_2={\left[\begin{array}{c}\frac{1}{2}\left(-k_1\left({\omega }_y^2+{\omega }_p^2\right)+k_3\left({\omega }_y^{\prime }+2{\omega }_r{\omega }_p\right)-k_5\left({\omega }_p^{\prime }-2{\omega }_r{\omega }_y\right)-2k_6{\omega }_p\right)\\ \frac{1}{2}\left(-k_1{\omega }_y^{\prime }-2k_2{\omega }_y-k_3\left({\omega }_y^2+{\omega }_r^2\right)+k_5\left({\omega }_r^{\prime }+2{\omega }_p{\omega }_y\right)+2k_6{\omega }_r\right)\end{array}\right]}^T \\ \times {\left[\begin{array}{cc}{k}_1& k_3\\ k_2+k_3{\omega }_y-k_5{\omega }_p& k_4-k_1{\omega }_y+k_5{\omega }_r\end{array}\right]}^{-1}\cdot \left[\begin{array}{c}{K}_0\\ K_1\end{array}\right]=\alpha K_0+\beta K_1 \end{gathered}$$

(29)

Thus,$R_r\left(t_a\right)$ can be re-written as

$$R_r\left(t_a\right)=K_0+K_1\cdot t_a+\left(\alpha K_0+\beta K_1\right)\cdot t_a^2$$

(30)

Therefore, the echoes of targets can be re-expressed as

$$S\left(f_r,t_a\right)={\sigma }_p\cdot \mathrm{w}\left(f_r\right)\mathrm{rect}\left(\frac{t_a}{T_a}\right)exp\left(-j\frac{4\pi }{c}\left(f_r+f_c\right)\left(K_0+K_1\cdot t_a+\left(\alpha K_0+\beta K_1\right)\cdot t_a^2\right)\right)$$

(31)

2.2. Signal Analysis for Targets with Complex Motion

Rewriting (31), and given by

$$\begin{gathered} S\left(f_r,t_a\right)={\sigma }_p\cdot \mathrm{w}\left(f_r\right)\mathrm{rect}\left(\frac{t_a}{T_a}\right)exp\left(-j\frac{4\pi f_r}{c}{K}_0\right)exp\left(-j\frac{4\pi K_1{t}_a}{\lambda }\right)exp\left(-j\frac{4\pi K_0}{\lambda }\right) \\ \times exp\left(-j\frac{4\pi \left(K_1\cdot t_a+\left(\alpha K_0+\beta K_1\right)\cdot t_a^2\right)·f_r}{c}\right)exp\left(-j\frac{4\pi \left(\left(\alpha K_0+\beta K_1\right)\cdot t_a^2\right)}{\lambda }\right) \end{gathered}$$

(32)

What is noteworthy is that, from (32), the first phase term is the range compression term. The second one is a linear phase term that is related to the Doppler frequency. The third one is a constant independent to ISAR imaging. The fourth one and the last one are, respectively, the range migration term and the 2-D spatial-variant phase error term [28]. It is quite obvious by now that the well-focused ISAR images can be obtained via compensating the range migration term and 2-D spatial-variant phase error term. The standard compensation algorithm can be applied to compensate for the range migration term, while the 2-D spatial-variant phase error terms that are different from scatterer to scatterer should be compensated for in every pixel because of the range and azimuth spatial-variant features, which increases the difficulty of phase error compensation for.

As described in [28], the analytical expression of the polluted signal consists of a focused ISAR imagery and azimuth phase history data, given by

$$S\left(f_r,t_a\right)=\tilde{S}\left(f_r,t_a\right)\cdot exp\left(-j\frac{4\pi \left(\alpha K_0+\beta K_1\right)\cdot t_a^2}{\lambda }\right)$$

(33)

where $\tilde{S}\left(f_r,t_a\right)$ is the focused ISAR imagery, and it is

$$\tilde{S}\left(f_r,t_a\right)=\mathrm{w}\left(f_r\right)\mathrm{rect}\left(\frac{t_a}{T_a}\right)exp\left(-j\frac{4\pi f_r}{c}{K}_0\right)exp\left(-j\frac{4\pi K_1{t}_a}{\lambda }\right)$$

(34)

Performing inverse Fourier transform (IFT) for (33) in terms of $f_r$, it can be writing a discrete form as

$$s\left(m,n\right)=\tilde{s}\left(m,n\right)\cdot exp\left(-j\frac{4\pi \left(\alpha K_0+\beta K_1\right)\cdot n^2}{\lambda }\right)$$

(35)

where $m=0,2,\cdot \cdot \cdot \left(M-1\right)$ denotes the range indices and $M$ denotes the number of samples in range dimension, $n=0,2,\cdot \cdot \cdot \left(N-1\right)$ denotes the azimuth index in synthetic aperture time and $N$ denotes the number of azimuth dimension, $s\left(m,n\right)$ and $\tilde{s}\left(m,n\right)$ are the discrete forms of $s\left(t_r,t_a\right)$ and $\tilde{s}\left(t_r,t_a\right)$, respectively.

It should be noted that unless the azimuth phase history data is perfect compensated for, or we cannot obtain a well-focused ISAR imagery, given by

$$\tilde{s}\left(m,n\right)=s\left(m,n\right)\cdot exp\left(-j\frac{4\pi \left(\alpha K_0+\beta K_1\right)\cdot n^2}{\lambda }\right)$$

(36)

Suppose the parameters $\left(\alpha ,\beta \right)$ are exactly estimated, then the well-focused ISAR images $g\left(m,n\right)$ can be obtained via an inverse discrete Fourier transform (IDFT) along with azimuth dimension, and it is

$$\begin{gathered} g\left(m,n\right)=\frac{1}{N}{\displaystyle \sum }_{n=0}^N\tilde{s}\left(m,n\right)\cdot exp\left(j2\pi \frac{mn}{N}\right) \\ =\frac{1}{N}{\displaystyle \sum }_{n=0}^{N-1}exp\left(j2\pi \frac{mn}{N}\right)·s\left(m,n\right)\cdot exp\left(j\frac{4\pi \left(\alpha K_0+\beta K_1\right)\cdot n^2}{\lambda }\right) \end{gathered}$$

(37)

Once the 2-D spatial-variant phase error terms are accurately compensated for, the well-focused ISAR images can be obtained. Thus, we suppose the 2-D spatial-variant phase error terms are perfectly compensated for, and conduct IFT and FT along with $f_r$ and $t_a$, respectively, to (33), the well-focused ISAR imagery can be expressed as

$$\begin{gathered} \tilde{g}\left(m,n\right)=IFF{T}_{t_a}\left\{{FFT}_{f_r}\left\{\tilde{S}\left(f_r,t_a\right)\right\}\right\} \\ ={\sigma }_p\cdot sinc\left[\frac{2B_r}{c}\left(m-K_0\right)\right]\cdot sinc\left[\frac{2T_a}{\lambda }\left(n-K_1\right)\right] \end{gathered}$$

(38)

where $IFF{T}_{t_a}\{\cdot \}$ and $FF{T}_{f_r}\{\cdot \}$ denote the IFT and FT operator to $t_a$ and $f_r$, respectively. Obviously, based on the resolution relationship in range dimension and cross-range dimension, $K_0$ and $K_1$ can also be expressed as

$$K_0=\frac{c}{2B_r}\cdot m$$

(39)

$$K_1=\frac{\lambda }{2T_a}\cdot n$$

(40)

Substituting (39) and (40) to (37), it can be re-written as

$$g\left(m,n\right)=\frac{1}{N}{\displaystyle \sum }_{n=0}^{N}exp\left(j2\pi \frac{mn}{N}\right)·s\left(m,n\right)\cdot exp\left(j\left(\frac{2\pi c}{\lambda B_r}\alpha \cdot m+\frac{2\pi }{T_a}\beta \cdot h\right)n^2\right)$$

(41)

It is noticeable that once the 2-D spatial-variant phase error terms are precisely compensated for, the well-focused ISAR images can be obtained. Furthermore, the 2-D spatial-variant phase errors are determined by two unknown parameters $\left(\alpha ,\beta \right)$. Therefore, many existing optimization algorithms such as gradient-based Newton algorithms or extended algorithms can be utilized to estimate those parameters. However, those methods need a large processing cost. In addition, the selection of initial values largely determines the imaging efficiency and accuracy of those algorithms. As a result, a fast parameters estimation method is still needed to improve computational efficiency.

2.3. Proposed Approach Description

In this section, to improve the ISAR imaging efficiency, the Generalized Radon-Fourier transform (GRFT) is firstly utilized to coarsely estimate the two unknown parameters for a suitable initial value to fine global optimal estimation. Following the coarse estimation operation, the fine search operation via gradient-based descent algorithm, e.g., Broyden–Fletcher–Goldfarb–Shanno (BFGS) method, is conducted. Additionally, the image entropy based on subarray averaging operation is generated to accelerate the global optimal convergence.

2.3.1. Coarse Parameters Estimation with GRFT

In general, the GRFT [29] can be defined as

$$G\left(\alpha ,\beta \right)={{\displaystyle \int }}_{-T_{ob}/2}^{T_{ob}/2}S\left(f_r,t_a,\alpha ,\beta \right)\cdot H\left(f_r,t_a,\alpha ,\beta \right)d{t}_a$$

(42)

where $T_{ob}$ is the observation time, $\alpha \in \left[{\alpha }_{min},{\alpha }_{max}\right]$, $\beta \in \left[{\beta }_{min},{\beta }_{max}\right]$ are the definition domains of parameters $\left(\alpha ,\beta \right)$, and ${\alpha }_{min},{\alpha }_{max}, {\beta }_{min},{\beta }_{max}$ are the maximum and minimum values of $\alpha ,\beta$, respectively. $G\left(\alpha ,\beta \right)$ denotes the coherent peak of GRFT. The analytical expression of $H\left(f_r,t_a,\alpha ,\beta \right)$ can be written as

$$H\left(f_r,t_a,\alpha ,\beta \right)=exp\left(j\frac{4\pi }{\lambda }\left(\alpha \cdot K_0+\beta \cdot K_1\right)\cdot t_a^2\right)$$

(43)

It should be noted that, after GRFT processing, the signal of (36) are now projected to the parameters domain, where the sole coherent peak value is obtained once the real value is around the estimated value $\left(\widehat{\alpha },\widehat{\beta }\right)$, given by

$$G\left(\widehat{\alpha },\widehat{\beta }\right)={\sigma }_p\cdot B_r\cdot T_p\cdot T_{ob}\cdot PRF$$

(44)

where $PRF$ is the pulse repetition frequency (PRF).

To further describe the principle of GRFT, the sketch map is shown in Figure 2. For the sake of simplicity, it is assumed that the component sets $H\left(f_r,t_a,\alpha ,\beta \right)$ that calculated using different $\alpha$ and $\beta$ are presented in Figure 2a, where $\alpha \in \left[{\alpha }_{min},{\alpha }_{max}\right]$ with step $△\alpha$, and $\beta \in \left[{\beta }_{min},{\beta }_{max}\right]$ with step $△\beta$. The mapping results $of G\left(\widehat{\alpha },\widehat{\beta }\right)$ with the GRFT method are provided in Figure 2b. Therefore, if the real values are around $\left(\widehat{\alpha },\widehat{\beta }\right)$, the coherent peak $G\left(\widehat{\alpha },\widehat{\beta }\right)$ that larger than others can be detected, shown in Figure 2b.

Therefore, according to the position and value of the maximal coherent peak, the true value is thought to be around it and the general range of $\left(\alpha ,\beta \right)$ can be estimated as

$$\left(\widehat{\alpha },\widehat{\beta }\right)=\arg max\left(G\left(\alpha ,\beta \right)\right)$$

(45)

where $\left(\widehat{\alpha },\widehat{\beta }\right)$ denotes the estimated parameters around the true one. Besides, the GRFT method can be repeatedly utilized to narrow down the range of the parameters to be estimated.

2.3.2. Fine Parameters Estimation with Gradient-Based Optimal

Image Entropy Combined with Subarray Averaging Operation

The entropy refers to the degree of chaos in the system, and the entropy of a well-focused ISAR image is littler than that of the unfocused one. Hence, the littler the image entropy is, the clearer the images will become. Therefore, for existing parameters estimation methods based on optimization technique, image entropy [30] is widely utilized as an image quality indicator to evaluate the ISAR imaging performance, which is an effective method to signify the definition of ISAR image. However, the cost surface obtained via conventional image entropy contains many local optimal and global optima solutions, which increase the difficulty for the search of the global optimum solution because the convergence of optimal solution cannot be guaranteed. To overcome the barrier above, many optimization methods, such as evolutionary computation (EC), simulated annealing (SA), ant colony optimization, are utilized to find the global optimum solution. However, those methods require a large processing time, which is inapplicable in practical applications [31]. Additionally, their solutions are extremely susceptible to the number of search spaces. Taking the difficulty above into consideration, in this work, the subarray averaging combined with image entropy [32,33] is proposed to eliminate the local optimal solution.

Subsequently, the image entropy combined with subarray averaging is defined as

$$E=-{\displaystyle \sum }_{l=0}^{L}h\left(l\right)\ln h\left(l\right)$$

(46)

$$h\left(l\right)=\frac{{{\displaystyle \sum }}_{p=1}^P{\left|g_p\left(l\right)\right|}^2}{{{\displaystyle \sum }}_{p=1}^P{{\displaystyle \sum }}_{l=0}^{L-1}{\left|g_p\left(l\right)\right|}^2}$$

(47)

where $g_p\left(l\right)$ is $pth$ subarray image in (47), $L$ and $P$ denote the length and total number of the subarray sub-image, respectively, $h\left(l\right)$ denotes the summation of squared sub-image for $P$ subarrays, normalized by total energy calculated in $P$-squared sub-image. Therefore, $h\left(l\right)$ can be applied as a new sub-image obtained via the subarray averaging technique proposed in this work. As shown in Figure 3, the sketch map of the subarray averaging is provided, where the ISAR image is divided into many overlapped subarrays sub-images with length $L$ along with azimuth dimension, the interval of consecutive subarrays is $d$, shown in Figure 3. Besides, the cost surface using conventional image entropy and subarray averaging combined with image entropy are provided in Figure 4a,b, respectively, where the simulated parameters are the same as provided in Section 3.1. A notable feature, from Figure 4, is that the global minimum solution and local minimum solution coexist in the cost surface calculated using conventional image entropy compared with that of the image entropy combined with subarray averaging, shown in Figure 4a,b, respectively. Thus, the global optimal solution cannot be guaranteed while the local extremum solutions existed in the cost surface. Therefore, the image entropy combined with subarray averaging is generated as the cost surface for optimal searching, which would significantly improve the ISAR imaging efficiency.

Parameters Estimation Based on Gradient Descent Method

In Section 3.1, the approximate solution of the parameters $\left(\alpha ,\beta \right)$ is estimated. As a result, an effective method should be adopted to further accurately estimate parameters $\left(\alpha ,\beta \right)$. Thanks to the smooth cost surface shown in Figure 4b, the global optimal solution can be quickly detected along the gradient descent direction of which. Therefore, for the existing method, gradient-based algorithms, e.g., Newton, are a valid approach to iteratively search the global optimal solutions. However, both the Newton method and the damped Newton method are time-consuming because the Hesse matrix should be calculated. Thanks to the Hesse matrix can be increasingly approximated by using the gradient information in each iteration, the quasi-Newton method [34] is effective for solving the unconstrained optimization problems. Therefore, in this work, the quasi-Newton method based on the BFGS algorithm [35,36] is adopted to search for the global optimal solutions.

Based on (46) and (47), the Formula for the partial derivative of $E\left(\alpha ,\beta \right)$ to $\alpha$ is derived as

$$\begin{array}{ll}\frac{\partial E\left(\alpha ,\beta \right)}{\partial \alpha }& =-{\displaystyle {\displaystyle \sum }_{l=0}^L}\left(\frac{\partial g\left(l,\alpha ,\beta \right)}{\partial \alpha }\ln g\left(l,\alpha ,\beta \right)+g\left(l,\alpha ,\beta \right)\frac{\partial \ln g\left(l,\alpha ,\beta \right)}{\partial \alpha }\right)\\ & =-{\displaystyle {\displaystyle \sum }_{l=0}^L}\left(\frac{\partial g\left(l,\alpha ,\beta \right)}{\partial \alpha }\ln g\left(l,\alpha ,\beta \right)+\frac{\partial g\left(l,\alpha ,\beta \right)}{\partial \alpha }\right)\\ & =-{\displaystyle {\displaystyle \sum }_{l=0}^L}\left(\frac{\partial g\left(l,\alpha ,\beta \right)}{\partial \alpha }\left(\ln g\left(l,\alpha ,\beta \right)+1\right)\right)\end{array}$$

(48)

where $\frac{\partial g\left(l,\alpha ,\beta \right)}{\partial \alpha }$ can be expressed as

$$\begin{array}{l}\frac{\partial g\left(l,\alpha ,\beta \right)}{\partial \alpha }\\ =\frac{\left({{\displaystyle \sum }}_{p=0}^{P-1}\frac{\partial {\left|h_p\left(l,\alpha ,\beta \right)\right|}^2}{\partial \alpha }\right)\left({{\displaystyle \sum }}_{p=0}^{P-1}{{\displaystyle \sum }}_{l=0}^{L-1}{\left|h_p\left(l\right)\right|}^2\right)-{{\displaystyle \sum }}_{p=0}^{P-1}{\left|h_p\left(l,\alpha ,\beta \right)\right|}^2{{\displaystyle \sum }}_{p=0}^{P-1}{{\displaystyle \sum }}_{l=0}^{L-1}\frac{\partial {\left|h_p\left(l\right)\right|}^2}{\partial \alpha }}{{\left({{\displaystyle \sum }}_{p=0}^{P-1}{{\displaystyle \sum }}_{l=0}^{L-1}{\left|h_p\left(l\right)\right|}^2\right)}^2}\end{array}$$

(49)

where $\frac{\partial {\left|h_p\left(l,\alpha ,\beta \right)\right|}^2}{\partial \alpha }$ can be written as

$$\frac{\partial {\left|h_p\left(l,\alpha ,\beta \right)\right|}^2}{\partial \alpha }=2\cdot Re\left\{h_p^{*}\left(l,\alpha ,\beta \right)\cdot \frac{\partial h_p\left(l,\alpha ,\beta \right)}{\partial \alpha }\right\}$$

(50)

where $Re\{·\}$ denotes the real parts. $\frac{\partial h_p\left(l,\alpha ,\beta \right)}{\partial \alpha }$ can be written as

$$\frac{\partial h_p\left(l,\alpha ,\beta \right)}{\partial \alpha }=\frac{1}{N}{\displaystyle \sum }_{n=0}^{N-1}exp\left(j2\pi \frac{mn}{N}\right)s\left(m,n\right)exp\left(j\frac{4\pi }{\lambda }\left({\zeta }^x{K}_0+{\zeta }^y{K}_1\right)\cdot n^2\right)\left(j\frac{4\pi }{\lambda }{K}_0\cdot n^2\right)$$

(51)

Similar to the derivation of $\alpha$, the partial derivative of $E\left(\alpha ,\beta \right)$ to $\beta$ can also be derived. The sole difference of $\frac{\partial h_p\left(l,\alpha ,\beta \right)}{\partial \beta }$ is as

$$\frac{\partial h_p\left(l,\alpha ,\beta \right)}{\partial \beta }=\frac{1}{N}{\displaystyle \sum }_{n=0}^{N-1}exp\left(j2\pi \frac{mn}{N}\right)s\left(m,n\right)exp\left(j\frac{4\pi }{\lambda }\left({\zeta }^x{K}_0+{\zeta }^y{K}_1\right)\cdot n^2\right)\left(j\frac{4\pi }{\lambda }{K}_1\cdot n^2\right)$$

(52)

Thus, the gradient $\nabla E\left(\alpha ,\beta \right)$ of entropy $E\left(\alpha ,\beta \right)$ to $\left(\alpha ,\beta \right)$ can be expressed as

$$\nabla E\left(\alpha ,\beta \right)=\left[\frac{\partial E\left(\alpha ,\beta \right)}{\partial \alpha },\frac{\partial E\left(\alpha ,\beta \right)}{\partial \beta }\right]$$

(53)

Furthermore, the detailed implementation procedure of BFGS is as

• Set the initial parameter $x^0$, initial matrix $H_0$, and the precision of error $\epsilon$ as $\left(\widehat{\alpha },\widehat{\beta }\right)$, unit matrix, and $1{\mathrm{e}}^{-5}\ge 0$, respectively.
• Calculate the gradient $\nabla E\left(x^0\right)$. If $\Vert \nabla E\left(x^0\right)\Vert <\epsilon$, then stop the calculation and the optimal parameter is $x^{*}=x^0$. Otherwise, conduct the next step.
• Set $p_0=-H_0\cdot \nabla E\left(x^0\right),k\leftarrow 0$, and conduct the next step.
• Perform a one-dimensional search to obtain $t_k$ such that $E\left(x^k+t_k\cdot p_k\right)=\min E\left(x^k+t\cdot p_k\right)$ is satisfied. Set $x^{k+1}=x^k+t\cdot p_k$ and conduct the next step.
• Calculate $\nabla E\left(x^{k+1}\right)$, if $\Vert \nabla E\left(x^{k+1}\right)\Vert \le \epsilon$, then stop the calculation, and set the optimal parameter as $x^{*}=x^{k+1}$. Otherwise, conduct the next step.
• If $k+1=n$, then $x^0=x^n$ and conduct Step 3. Otherwise, conduct the next step.
• Calculate

  $$\begin{array}{ll}{H}_{k+1}& =H_k+\frac{∆x_k·∆x_k^T}{∆x_k^T\cdot ∆g_k}\left(1+\frac{∆g_k^T\cdot H_k\cdot ∆g_k}{∆x_k^T\cdot ∆g_k}\right)\\ & -\frac{1}{∆x_k^T\cdot ∆g_k}\left(∆x_k\cdot ∆g_k^T\cdot H_k+H_k\cdot ∆g_k\cdot ∆x_k^T\right)\end{array}$$

  (54)

  $$p_{k+1}=-H_{k+1}\cdot \nabla E\left(x^{k+1}\right)$$

  (55)

  where $∆x_k=x^{k+1}-x^k,∆g_k=E\left(x^{k+1}\right)-E\left(x^k\right)$, ${(\cdot )}^T$ denotes the transposition of $(\cdot )$, and $k\leftarrow k+1$, and repeat to step 4.

Finally, the whole ISAR imaging procedure of our proposed method is as follows, and the flowchart of the proposed method is provided in Figure 5.

• Obtain the raw echoes, and conduct preprocessing part, e.g., range alignment, phase adjustment, and RCM correction.
• Coarsely search the range of true parameters via detecting the coherent peak with GRFT.
• Finely estimate the optimal parameters by using the gradient descent method.
• Finishing the compensation of 2D spatial-variant phase errors and obtain the well-focused ISAR image.

2.4. Some Considerations for the Proposed Method in Applications

2.4.1. Computational Complexity Analysis

In this section, the computational complexity of our proposed method is compared with that of the PGA [37] method, STFT. Without loss of generality, the operations with less time cost are omitted. Generally, we suppose $N{log}_{2}N$ floating-point operations (FLOPs) are required to compute $N$-point FT or IFT, and $N$ FLOPs are needed to compute $N$-point data center shifting operation, $MN$ FLOPs are required to calculate the image entropy for ISAR image with size $\left(M\times N\right)$, and $N$ FLOPs are needed for the STFT of $N$-point data, and $\left(MN\right)$ FLOPs are needed to calculate gradient for ISAR image with size $\left(M\times N\right)$. Let $N_r,N_a,P,Q$ represent the samples in range dimension and azimuth dimension, the search number for the search space $\left(\alpha ,\beta \right)$, respectively.

The realization procedure of the PGA algorithm consists of five steps [36], which are coarsely RD imaging, center shifting operation, windowing operation, phase gradient estimation, and iterative phase error correction. The computational complexity of coarsely RD imaging that is composed of range compression and azimuth compression is $O\left\{N_a{N}_r{log}_2\left(N_r\right)+N_a{N}_r+N_a{N}_r{log}_2\left(N_a\right)\right\}$. The computational cost of center shifting operation, windowing operation, phase gradient estimation, and iterative phase error correction are $O\left\{N_a{N}_r\right\}$,$O\left\{N_a{N}_r{log}_2\left(N_a\right)\right\}$,$O\left\{N_a{N}_r+N_a{N}_r{log}_2\left(N_a\right)\right\}$, respectively. Thus, the total computational complexity of the PGA algorithm is

$$C_{\mathrm{PGA}}=O\left\{N_a{N}_r{log}_2\left(N_r\right)+2N_a{N}_r+3N_a{N}_r{log}_2\left(N_a\right)\right\}$$

(56)

According to [16], the signal after pulse compression is processed by the STFT method. Hence, the procedure consists of range compression, and iteratively STFT processing along with azimuth dimensional in each range cell. In this paper, the length of the frequency smoothing window of STFT is $N_a/4$. Thus, the computational cost of time-frequency processing for $N_a$ points data is $N_a\cdot N_a/8+N_a{log}_2\left(N_a\right)$. Therefore, the total computational complexity of STFT is

$$C_{\mathrm{STFT}}=O\left\{N_a{N}_r{log}_2\left(N_a\right)+2N_a{N}_r+N_r\left( N_a{N}_a/8+N_a{log}_2\left(N_a\right)\right)\right\}$$

(57)

Based on the imaging procedure shown in Figure 5, the computational cost of range compression is $O\left(N_a{N}_r{log}_2\left(N_r\right)+N_a{N}_r\right)$. Suppose the number of subarray sub-image is $P$, and the size of the subarray image is $N_r\times L$, and the iterative time of the BFGS procedure is $T_I$. Considering the algorithm based on image entropy combined with subarray averaging operation, the computation of the gradient in (30) consumes $O\left(P{N}_{r}L\right)$. The considerable time-consuming procedure is the decision of the step size using GRFT in each iteration, which needs $\left(N_a{N}_r\right)$ FLOPs complex multiplications and $O\left( N_r\left(N_a\cdot {log}_2\left(N_a\right)\right)\right)$ FLOPs FT operations. Suppose the number of conducting GRFT method is $J$. Thus, the all computational load of the proposed method is

$$C_{prop}=O\left\{\begin{array}{c}{N}_a{N}_r{log}_2\left(N_r\right)+3N_a{N}_r+N_r\left(N_a{log}_2\left(N_a\right)\right)\\ +P{N}_{r}L{T}_I+J\left(N_a{N}_r{log}_2\left(N_a\right)+N_a{N}_r\right)\end{array}\right\}$$

(58)

According to the analysis above, the computational complexity of our proposed method is larger than that of the PGA. In addition, due to time-frequency processing in all range-cell, the computational complexity of the STFT method is the largest. Though the computational complexity of PGA is superior to that of our proposed method, the imaging quality of the PGA’s is poor, which have limitation in a real application. Considering the tradeoff between imaging quality and imaging efficiency, the proposed method has superiority in contrast to that of the PGA and STFT methods.

2.4.2. Doppler Frequency Spectrum Analysis

In this section, the time-varying Doppler frequency spectrum is analyzed. To illustrate the Doppler spectrum, suppose targets have yaw motion only, e.g.,${\theta }_r={\theta }_p=0$. Thus, the rotation matrix becomes

$$rot\left(t_a\right)=\left[\begin{array}{ccc}cos\left({\theta }_y\right)& -sin\left({\theta }_y\right)& 0\\ sin\left({\theta }_y\right)& cos\left({\theta }_y\right)& 0\\ 0& 0& 1\end{array}\right]$$

(59)

Substituting (59) into (5), then the rotational motion part ${ℛ}_r\left(t_a\right)$ can be written as

$$\begin{array}{ll}{ℛ}_r\left(t_a\right)=& i_1\left(x_p\left(1-\frac{1}{2}\left({\omega }_y{t}_a+\frac{1}{2}{\omega }_y^{\prime }{t}_a^2\right)\right)-y_p\left({\omega }_y{t}_a+\frac{1}{2}{\omega }_y^{\prime }{t}_a^2\right)\right)\\ & +i_2\left(x_p\left({\omega }_y{t}_a+\frac{1}{2}{\omega }_y^{\prime }{t}_a^2\right)+y_p\left(1-\frac{1}{2}\left({\omega }_y{t}_a+\frac{1}{2}{\omega }_y^{\prime }{t}_a^2\right)\right)\right)\end{array}$$

(60)

Thus, the phase of the returned signal from targets can be re-written as

$$\mathsf{\Phi }\left(t_a\right)=\frac{4\pi \left(f_c+f_r\right){ℛ}_r\left(t_a\right)}{c}$$

(61)

Conducting the derivation of $\mathsf{\Phi }\left(t_a\right)$ in terms of $t_a$, the Doppler frequency spectrum of the targets can be derived as

$$\begin{array}{ll}{f}_d\left(t_a\right)& =-\frac{1}{2\pi }\frac{d\left\{\mathsf{\Phi }\left(t_a\right)\right\}}{d{t}_a}\\ & =\frac{2\left(f_c+f_r\right)}{c}\frac{d\left\{ℛ\left(t_a\right)\right\}}{d{t}_a}\end{array}$$

(62)

where

$$\begin{gathered} \frac{d\left\{\mathsf{\Phi }\left(t_a\right)\right\}}{d{t}_a}=i_1^{\prime }\left\{x_p\left(1-\frac{1}{2}\left({\omega }_y{t}_a+\frac{1}{2}{\omega }_y^{\prime }{t}_a^2\right)\right)-y_p\left({\omega }_y{t}_a+\frac{1}{2}{\omega }_y^{\prime }{t}_a^2\right)\right\} \\ +i_1\left(-\frac{x_p}{2}\left({\omega }_y+{\omega }_y^{\prime }{t}_a\right)-y_p\left({\omega }_y+{\omega }_y^{\prime }{t}_a\right)\right) \\ +i_2^{\prime }\left\{x_p\left(1-\frac{1}{2}\left({\omega }_y{t}_a+\frac{1}{2}{\omega }_y^{\prime }{t}_a^2\right)\right)-y_p\left({\omega }_y{t}_a+\frac{1}{2}{\omega }_y^{\prime }{t}_a^2\right)\right\} \\ +i_2\left(x_p\left({\omega }_y+{\omega }_y^{\prime }{t}_a\right)-\frac{y_p}{2}\left({\omega }_y+{\omega }_y^{\prime }{t}_a\right)\right) \end{gathered}$$

(63)

where $i_j^{\prime },j=1,2$ denote the derivation of $i_j,j=1,2$.

Similarly, if the radar LOS is a constant during CPI, then the Doppler frequency $f_d^{con}\left(t_a\right)$ becomes

$$\begin{array}{ll}{f}_d^{con}\left(t_a\right)& =-\frac{1}{2\pi }\frac{d\left\{{\mathsf{\Phi }}^{con}\left(t_a\right)\right\}}{d{t}_a}\\ & =\frac{2\left(f_c+f_r\right)}{c}\frac{d\left\{{ℛ}^{con}\left(t_a\right)\right\}}{d{t}_a}\end{array}$$

(64)

where

$$\begin{gathered} {ℛ}^{con}\left(t_a\right)=l_1\left(x_p\left(1-\frac{1}{2}\left({\omega }_y{t}_a+\frac{1}{2}{\omega }_y^{\prime }{t}_a^2\right)\right)-y_p\left({\omega }_y{t}_a+\frac{1}{2}{\omega }_y^{\prime }{t}_a^2\right)\right) \\ +l_2\left(x_p\left({\omega }_y{t}_a+\frac{1}{2}{\omega }_y^{\prime }{t}_a^2\right)+y_p\left(1-\frac{1}{2}\left({\omega }_y{t}_a+\frac{1}{2}{\omega }_y^{\prime }{t}_a^2\right)\right)\right) \end{gathered}$$

(65)

where $l_1=cos\left({\phi }_p\right)cos\left({\eta }_p\right),l_2=cos\left({\phi }_p\right)sin\left({\eta }_p\right)$.

Accordingly, the Doppler frequency spectrum for targets with roll or pitch only has the same expression form. Additionally, comparing with the assumption that the radar LOS is a constant, from (62)–(65), the Doppler frequency has a higher-order phase term which is determined by the higher-order coefficients of the radar LOS that would affect the Doppler frequency spectrum. Therefore, the Doppler frequency spectrum calculated using time-varying radar LOS can accurately present the non-stationary IPP.

2.4.3. Sampling Rate and PRF

According to (61), the Doppler frequency $f_d\left(t_a\right)$ of targets can be written as

$$\begin{array}{ll}{f}_d\left(t_a\right)& =-\frac{1}{2\pi }\frac{d\left\{\mathsf{\Phi }\left(t_a\right)\right\}}{d{t}_a}\\ & =\frac{2}{\lambda }\left(K_1+2K_2\cdot t_a\right)\end{array}$$

(66)

Based on the Nyquist sampling theorem, the pulse repetition frequency (PRF) must satisfy $PRF\ge 2f_{dmax}$, where $f_{dmax}$ is the maximum of the Doppler frequency, and it is

$$\begin{gathered} PRF\ge 2f_{dmax} \\ =max\left\{\frac{4}{\lambda }\left(K_1+2K_2\cdot T_a\right)\right\} \end{gathered}$$

(67)

where $max\{\cdot \}$ denotes the maximum value related to the maximum size of the targets.

2.4.4. Phase Error Analysis

For the instantaneous slant range $R_r\left(t_a\right)$, given by

$$R_r\left(t_a\right)={\left[rot\left(t_a\right)\cdot \mathit{P}\right]}^T\cdot i_{los}$$

(68)

The approximation of second-order Taylor series expression will cause the phase error $\epsilon$, and it can be written as

$$\epsilon =\frac{4\pi }{\lambda }\left|R_r^1-R_r^2\right|$$

(69)

where $R_r^1$ and $R_r^2$, respectively, denote the instantaneous slant range before and after Taylor series expression, and it can be written as

$$\begin{gathered} R_r^1=x_{p}sin\left({\eta }_p\right)\left(sin\left({\theta }_r\right)sin\left({\theta }_y\right)-cos\left({\theta }_r\right)cos\left({\theta }_y\right)sin\left({\theta }_p\right)\right) \\ +y_{p}sin\left({\theta }_p\right)\left(cos\left({\theta }_y\right)sin\left({\theta }_r\right)+cos\left({\theta }_r\right)sin\left({\theta }_p\right)sin\left({\theta }_y\right)\right) \\ +x_{p}cos\left({\eta }_p\right)sin\left({\phi }_p\right)\left(cos\left({\theta }_r\right)sin\left({\theta }_y\right)+cos\left({\theta }_y\right)sin\left({\theta }_p\right)sin\left({\theta }_r\right)\right) \\ +y_{p}cos\left({\eta }_p\right)sin\left({\phi }_p\right)\left(sin\left({\theta }_r\right)cos\left({\theta }_y\right)-sin\left({\theta }_p\right)sin\left({\theta }_r\right)sin\left({\theta }_y\right)\right) \\ +x_{p}cos\left({\eta }_p\right)cos\left({\phi }_p\right)cos\left({\theta }_p\right)cos\left({\theta }_y\right)-y_{p}cos\left({\eta }_p\right)cos\left({\phi }_p\right)cos\left({\theta }_p\right)sin\left({\theta }_y\right) \end{gathered}$$

(70)

$$\begin{gathered} R_r^2=k_1{x}_p+k_3{y}_p+\left(\left(k_2+k_3{\omega }_y-k_5{\omega }_p\right)x_p+\left(k_4-k_1{\omega }_y+k_5{\omega }_r\right)y_p\right)·t_a \\ +\left(\begin{array}{c}\frac{1}{2}\left(-k_1\left({\omega }_y^2+{\omega }_p^2\right)+k_3\left({\omega }_y^{\prime }+2{\omega }_r{\omega }_p\right)-k_5\left({\omega }_p^{\prime }-2{\omega }_r{\omega }_y\right)-2k_6{\omega }_p\right)x_p\\ ++\frac{1}{2}\left(-k_1{\omega }_y^{\prime }-2k_2{\omega }_y-k_3\left({\omega }_y^2+{\omega }_r^2\right)+k_5\left({\omega }_r^{\prime }+2{\omega }_p{\omega }_y\right)+2k_6{\omega }_r\right)y_p\end{array}\right)\cdot t_a^2 \end{gathered}$$

(71)

Substituting (70) and (71) into (69), the phase error can be written as

$$\epsilon =\frac{4\pi }{\lambda }\left|\begin{array}{c}{x}_{p}sin\left({\eta }_p\right)\left(sin\left({\theta }_r\right)sin\left({\theta }_y\right)-cos\left({\theta }_r\right)cos\left({\theta }_y\right)sin\left({\theta }_p\right)\right)\\ +y_{p}sin\left({\theta }_p\right)\left(cos\left({\theta }_y\right)sin\left({\theta }_r\right)+cos\left({\theta }_r\right)sin\left({\theta }_p\right)sin\left({\theta }_y\right)\right)\\ +x_{p}cos\left({\eta }_p\right)sin\left({\phi }_p\right)\left(cos\left({\theta }_r\right)sin\left({\theta }_y\right)+cos\left({\theta }_y\right)sin\left({\theta }_p\right)sin\left({\theta }_r\right)\right)\\ +y_{p}cos\left({\eta }_p\right)sin\left({\phi }_p\right)\left(sin\left({\theta }_r\right)cos\left({\theta }_y\right)-sin\left({\theta }_p\right)sin\left({\theta }_r\right)sin\left({\theta }_y\right)\right)\\ +x_{p}cos\left({\eta }_p\right)cos\left({\phi }_p\right)cos\left({\theta }_p\right)cos\left({\theta }_y\right)-y_{p}cos\left({\eta }_p\right)cos\left({\phi }_p\right)cos\left({\theta }_p\right)sin\left({\theta }_y\right)\\ -k_1{x}_p-k_3{y}_p-\left(\left(k_2+k_3{\omega }_y-k_5{\omega }_p\right)x_p+\left(k_4-k_1{\omega }_y+k_5{\omega }_r\right)y_p\right)·t_a\\ -\left(\begin{array}{c}\frac{1}{2}\left(-k_1\left({\omega }_y^2+{\omega }_p^2\right)+k_3\left({\omega }_y^{\text{'}}+2{\omega }_r{\omega }_p\right)-k_5\left({\omega }_p^{\text{'}}-2{\omega }_r{\omega }_y\right)-2k_6{\omega }_p\right)x_p\\ ++\frac{1}{2}\left(-k_1{\omega }_y^{\text{'}}-2k_2{\omega }_y-k_3\left({\omega }_y^2+{\omega }_r^2\right)+k_5\left({\omega }_r^{\text{'}}+2{\omega }_p{\omega }_y\right)+2k_6{\omega }_r\right)y_p\end{array}\right)\cdot t_a^2\end{array}\right|$$

(72)

In general, the phase errors caused by the motion can be neglected if it is confined within $\frac{\pi }{4}$, given by

$$\epsilon \le \frac{\pi }{4}$$

(73)

To further present the phase error caused by the operation of approximation, the simulated result based on (72) is shown in Figure 6, where the simulated parameters are the same as Section 3.1. It is worthy to be noted that the phase errors of approximation are confined with $\frac{\pi }{4}$, shown in Figure 6.

3. Experimental Results and Analysis

In this section, several experiments and corresponding analyses using simulated data and electromagnetic data are presented to verify the availability and robustness of our proposed approach.

3.1. Simulation Results

In this part, some experimental results based on simulated echoes are presented to validate the effectiveness of our proposed method. The simplified ship model is shown in Figure 7, where the model consists of 73 scatterers located on the surface of the target. Suppose the linear frequency modulation (LFM) signal is transmitted, whose carrier frequency, bandwidth, pulse width, and pulse repetition frequency are 5 GHz, 500 MHz, 4us, and 1000 Hz, respectively. For the sake of simplicity, the amplitudes of echoes are all ones. The target motion parameters and radar LOS parameters are given in

Table 1. Motion parameters of the targets and radar LOS parameters.

Target motion parameters (rad/s)
${\omega }_y$	${\omega }_p$		${\omega }_r$	${\omega }_r^{\prime }$	${\omega }_y^{\prime }$		${\omega }_p^{\prime }$
−0.121	0.1		−0.01	3.5 × 10−2	5.5 × 10−2		−0.025
LOS motion parameter (rad/s)
${\phi }_0$		${\eta }_0$		$\kappa$		$\gamma$
1.0472		0.7854		0.35		0.5

. The number of pulses utilized for ISAR imaging is 640, and each pulse contains 3000 samples. All experiments are conducted on the platform of MATLAB 2014a on an Intel(R) Core(TM) i5-8400 CPU @ 2.8 GHz, 2808 MHz, RAM 8 GB, and Microsoft Window 10 operating system.

Based on the realization procedure of our proposed method, the GRFT is utilized to coarsely estimate the range of true parameters. In general, we can roughly calculate the range of $\alpha$ and $\beta$ according to the extreme motion condition of the ship targets. By doing so, in this work, the coarse result with the GRFT method is shown in Figure 8. Noticeably, the sole coherent peak is obtained from Figure 8, which means that the true parameters are around it. Thus, the coarse parameters can be obtained by detecting the coordinates of the sole coherent peak. Additionally, the fine parameters are estimated by using the BFGS method, and the convergence process of fine parameters estimation operation is provided in Figure 9, where the convergence of image entropy in terms of iteration number with BFGS method is compared with that of the proposed method. The initial value using the BFGS method is set as (1, 1), and the cost function is obtained with conventional image entropy. The search range of coarse estimation using the GRFT approach is (−1, 1) with step 0.25, and the cost function is calculated using the image entropy combined with subarray averaging. It is worthy to be pointed out that the convergence of iteration number with the proposed method and BFGS, respectively, are 4-time and 12-time. Our proposed method converges faster than that of the BFGS method. Thus, our proposed method can improve the convergent efficiency for global convergence.

And then, we are focusing on the imaging performance of PGA, STFT, and our proposed method. The ISAR imaging results using PGA, STFT, and our proposed method are provided in Figure 10a–c, respectively. In addition, the image entropy of different imaging methods is provided in

Table 2. Image entropy of corresponding imaging method.

Approach	Image Entropy
PGA algorithm	9.3565
STFT method	11.3896
Proposed method	8.0824

, where the image entropy of the imaging result with the proposed approach is littler than that of the PGA and STFT methods. To further demonstrate the focal performance, a scatterer marked with a circular is partially enlarged. Notably, the imaging performance of our proposed method is superior to that of the PGA and STFT methods.

3.2. Imaging Result under Different SNRs

In this section, different zero-mean complex Gaussian noises are added to the range-compressed echoes to verify the anti-noise performance of our proposed method, and the signal-to-noise ratio (SNR) is defined as

$$\mathrm{SNR}=10{log}_{10}\left(\frac{Mean power of signal}{Power of noise}\right)$$

(74)

The imaging results under different SNRs from 0 to −14 dB are provided in Figure 11, where Figure 11a–c, Figure 11d–f, Figure 11g–i are the imaging results via PGA, STFT, and our proposed method under SNR = 0 dB, SNR = −7 dB, SNR = −14 dB, respectively. It should be observed that, from Figure 11, the anti-noise performance of our proposed method is better than that of the RD, PGA method, especially at low SNR conditions.

Additionally, the conventional image entropy [30] is utilized as a quality indicator for imaging performance, and the image entropy $E$ is defined as

$$E=-{\displaystyle \sum }_{m=0}^{M-1}{\displaystyle \sum }_{n=0}^{N-1}\frac{{\left|g\left(m,n\right)\right|}^2}{G}log\frac{{\left|g\left(m,n\right)\right|}^2}{G}$$

(75)

where $G={\displaystyle \sum }_{m=0}^M{\displaystyle \sum }_{n=0}^N{\left|g\left(m,n\right)\right|}^2$, and $g\left(m,n\right)$ is the ISAR image. 50-time Monte Carlo experiments are performed to obtain the average value of image entropy, provided in Figure 12, where the SNRs are from −14 dB to 14 dB. It is seen from Figure 12, that the image entropy of our imaging results is lower than that of the PGA and STFT method, which means that the anti-noise performance of our proposed method is more obvious as the decrease of SNRs, as demonstrated in Figure 12.

3.3. Electromagnetic Data

To verify the performance of the proposed method in a more realistic environment, in this section, the radar cross-section (RCS) data are produced using an electromagnetic (EM) scattering prediction technique, which is an effective and economical way to obtain the radar echoes from ship target due to the practical difficulties in measuring real-world radar data. The well-known physical optical (PO) [38] technique which is one of the widely adopted techniques for high-frequency EM computation is utilized to generate the RCS data for ship targets. As shown in Figure 13, the scene contains a destroyer, where the entity and 3-D computer-aided design (CAD) is presented in Figure 13a,b, respectively. In addition, the model of original RCS data is provided in Figure 14a. For comparison, the imaging results using PGA, STFT, and our proposed method are provided in Figure 14b–d, respectively. Obviously, although three algorithms can roughly reconstruct the outline of the destroyer, the result of our proposed method is closer to the real outline in Figure 14a, especially for the head of the target which is labeled using a white box. Meanwhile, to further quantitative analysis for the imaging results, the conventional image entropy in (47) is utilized to evaluate the imaging results, provided in

Table 3. Image entropy of corresponding imaging method.

Approach	Image Entropy
PGA algorithm	8.7717
STFT method	10.1647
Proposed method	8.3391

. Notably, from Table 3, the image entropy of our proposed method is obviously littler than that of the PGA and STFT approach, which is consistent with the imaging results presented in Figure 14b–d.

4. Discussion

Inverse synthetic aperture radar (ISAR) plays an important role in target detection and recognition thanks to the all-weather, all-day, high-resolution. During imaging for ship targets with moderate motion characteristics, the IPP is modeled as a constant, which is accurate. However, if the ship targets are combined with a highly maneuvering motion, the assumption that the IPP is stationary during CPI is invalid and followed by the disabling of the existing imaging approach. Therefore, to accurately present the motion characteristic, in this work, the nonstationary IPP is introduced by using the time-varying instantaneous slant range where the radar LOS is modeled as a function about slow time. In addition, the 2-D spatial-variant phase errors which can severely blur ISAR images are derived. By exploring the relationship that the coordinate of scatterer can be described via the resolution in range dimension and azimuth dimension, the 2-D spatial-variant phase errors are estimated using two parameters that corresponding to the velocity of targets. Furthermore, the GRFT method and gradient-based optimal are proposed to coarsely and finely estimate those two parameters, respectively. Besides, the approach of image entropy combined with subarray averaging operation is generated to accelerate the global optimal convergence. In conclusion, thanks to the coarse estimation operation of the parameters via the GRFT method, the global convergence can be finished quickly, which can effectively improve the imaging efficiency. Meanwhile, the subarray averaging operation can eliminate the local optimal, which can not only improve the imaging efficiency but also ensures the accuracy of parameters estimation. The imaging results of our proposed method are compared with that of the PGA method and STFT method. The experimental results using simulations and electromagnetic data demonstrate that the proposed method has a tradeoff between imaging efficiency and imaging quality.

5. Conclusions

To reconstruct the ISAR images for ship targets with nonstationary IPP under complex 3-D rotation motion, an efficient imaging approach based on GRFT and gradient-based optimal is proposed in this work. First, the geometry and signal model for ship targets with nonstationary IPP is established. The 2-D spatial-variant phase errors caused by complex 3-D rotation motion are derived because they can seriously blur the ISAR images. Second, for the accuracy and effectiveness of compensating the 2-D spatial-variant phase errors, the GRFT in conjunction with gradient-based optimal is utilized to coarsely and finely estimate the motion parameters, respectively. Considering the local convergence of cost surface obtained with conventional image entropy, the image entropy combined with subarray averaging is introduced to improve the convergence efficiency for the global optimal solution. Finally, the 2-D spatial-variant phase errors can be precisely estimated followed by the well-focused ISAR images. Simulated data and EM data are adopted to verify the effectiveness of our proposed approach. In conclusion, the proposed method can take a tradeoff between imaging performance and computational efficiency under a low noisy environment.