Inverse Synthetic Aperture Radar Imaging of Targets with Complex Motion based on Optimized Non-Uniform Rotation Transform

: Focusing on the inverse synthetic aperture radar (ISAR) imaging of targets with complex motion, this paper proposes a modiﬁed version of the Fourier transform, called non-uniform rotation transform, to achieve cross-range compression. After translational motion compensation, the target’s complex motion is converted into non-uniform rotation. We deﬁne two variables—relative angular acceleration (RAA) and relative angular jerk (RAJ)—to describe the rotational nonuniformity. With the estimated RAA and RAJ, rotational nonuniformity compensation is carried out in the non-uniform rotation transform matrix, after which, a focused ISAR image can be obtained. Moreover, considering the possible deviation of RAA and RAJ, we design an optimization scheme to obtain the optimal RAA and RAJ according to the optimal quality of the ISAR image. Consequently, the ISAR imaging of targets with complex motion is converted into a parameter optimization problem in which a stable and clear ISAR image is guaranteed. Compared to precedent imaging methods, the new method achieves better imaging results with a reasonable computational cost. Experimental results verify the effectiveness and advantages of the proposed algorithm.


Introduction
Inverse synthetic aperture radar (ISAR) imaging is a key approach for the observation, surveillance, and recognition of space, aerial, and marine targets, with advantages of long working distance, high resolution, and all-day working capability [1][2][3][4]. The Fourier transform-based range-Doppler algorithm works well for imaging of targets with smooth motion, which has been widely used in ISAR imaging systems. In practice, however, the target's motion is often complex, such as fluctuation of ships and high maneuvering of airplanes, which leads to the blurring of ISAR images obtained through the range-Doppler algorithm. It is believed that the complex motion of the target results in quadratic and even cubic phase terms in cross-range signals; thus, the conventional range-Doppler algorithm fails. To solve this problem, several approaches have been proposed, including slow time interpolation [5], joint time-frequency representation [5][6][7][8][9][10], parameter estimation of signal models [11][12][13][14][15], and super resolution methods [16][17][18][19][20][21].
The slow time interpolation method aims to obtain the uniform rotation angle through interpolation, after which conventional range-Doppler algorithm can be carried out. Given that translational motion has been compensated, the target's complex motion can be regarded as non-uniform rotation. The idea of transforming non-uniform rotation into uniform rotation through suitable for imaging of a real target. Wang et al. [24,25], Li et al. [26], and Liu et al. [27] indicated that the non-uniform rotation of the target leads to multiple LFM or PPS components in cross-range signals, and used the slow time resampling approach for rotational nonuniformity compensation. However, the coefficient estimation of the PPS is not very accurate due to noise and mutual interference, which adversely affect the following slow time resampling and final ISAR imaging. In fact, slow time resampling is a type of interpolation that is costly in computation and inaccurate.
In this paper, we explore new approaches for ISAR imaging of targets with complex motion. Our studies focus on the estimation and compensation of rotational nonuniformity. Inspired by the parameter estimation methods of PPS, we further compensate the high-order phase terms caused by non-uniform rotation in the Fourier transform matrix. The Fourier transform containing high-order phase term compensation is called non-uniform rotation transform. Owing to the rigid body property of most targets, the high-order phase terms of PPS can be divided into several groups according to the scatterers' azimuth positions, and each group is compensated in the corresponding row of the Fourier transform matrix. The construction of the compensation phase terms depends on two parameters-relative angular acceleration (RAA) and relative angular jerk (RAJ)-which are determined by the motion status of the target and are identical for all scatterers. The rough values of RAA and RAJ are estimated through the coefficient estimation of PPS, and the accurate values are obtained through parameter optimization according to optimal imaging results. Moreover, the optimal ISAR image is unique for the selected pulses, thus the parameter optimization is a convex problem that can be solved using efficient methods such as gradient descent [28]. With non-uniform rotation transform, the ISAR imaging of targets with complex motion is equivalent to a convex optimization problem that guarantees an optimal imaging result without blurring. Besides clear and stable imaging results, ISAR imaging based on non-uniform rotation does not require heavy computation due to the rough coefficient estimation of PPS and the efficient optimization method.
The remainder of this paper is organized as follows. Section 2 describes the ISAR imaging model for targets with rigid body, and discusses the influence of the target's complex motion. Section 3 analyzes the properties of PPS in the cross-range direction and rotational nonuniformity estimation. Section 4 elaborates our new imaging method based on the optimized non-uniform rotation transform. Section 5 validates the proposed method through experiments, and Section 6 presents the conclusions. Figure 1 illustrates the ISAR imaging model of targets with rigid body. After translational motion compensation, the target's motion is transformed into rotation around an axis. We establish the target coordinate system o-xyz with the origin located at the target's centroid, where the z-axis is assumed to be the target's equivalent rotation axis and the y-axis is parallel to the line of sight (LOS) of the radar. While imaging, the target's rotation angle relative to the LOS is denoted as θ(t m ), where t m represents slow time. The distance between radar and target's centroid is denoted as R 0 . A scatterer in the target, denoted as q, is selected, whose coordinates in the target coordinate system are q x q , y q , z q . The distance between the radar and scatterer q can be expressed as

ISAR Imaging Model for Targets with Rigid Body
(1) The Taylor expressions of sin θ(t m ) and cos θ(t m ) are where o θ k (t m ) represents the sum of high-order terms of θ(t m ), whose minimum order is k + 1. Since θ(t m ) is small during the coherent processing time, the quadratic and higher-order terms in (2) and (3) can be ignored. Then, R q (t m ) is approximated as For targets with smooth motion, the uniform rotation model is used to describe θ(t m ) after translational motion compensation. The equivalent angular velocity is denoted as ω , and R q (t m ) becomes In general, an LFM signal is adopted as the transmitted signal, which can be expressed as where T p is the pulse width, f c is the carrier frequency, and K is the frequency modulation rate. The bandwidth of the signal is KT p . The rectangular window function rect(·) satisfies The echo is received through wideband direct sampling, and then, the matched filtering containing high-speed compensation [29] is adopted for pulse compression. After that, the obtained range profile is expressed as where t is fast time, Q is the total number of scatterers, A i is the backward scattering coefficient of the ith scatterers, and c is light velocity. Function sin c(·) satisfies sin c(x) = sin(πx) πx . Assuming that translational motion compensation has been accomplished, the phase at the peak of the sinc function corresponding to scatterer q is Clearly, ϕ q (t m ) contains two parts: constant phase component and linear phase component. The Fourier transform in the cross-range direction can be used to separate the scatterers in the Doppler frequency domain. The Doppler frequency of scatterer q is f m = − 2 f c c ω x q , where x q represents its azimuth position. We denote the range profile sequence after translational motion compensation as G(m, n), where m represents the pulse number related to t m and n represents the range cell related to t . There are M pulses and N range cells in total. We denote the discrete ISAR image as I(m, n). The Fourier transform matrix is denoted as P. The process of the range-Doppler algorithm can be expressed as For targets with complex motion, however, a non-uniform rotation model is used to describe θ(t m ). In general, the angular acceleration, denoted as α , and the angular jerk, denoted as β , are considered. Herein, the target's rotation angle is expressed as Therefore, ϕ q (t m ) becomes [25][26][27] Equation (12) indicates that ϕ q (t m ) satisfies the PPS model. The target's complex motion results in quadratic and cubic phase terms in (12), and the corresponding Doppler frequency Consequently, Doppler spreads arise and the ISAR image blurs. In the following section, we analyze and compensate the high-order phase terms of the PPS through non-uniform rotation transform.

Properties of the PPS and Rotational Nonuniformity Estimation
In ISAR imaging, the properties of PPS in the cross-range direction need to be addressed. According to (12), the linear, quadratic, and cubic coefficients of a PPS are p depend on the target's motion status, and are identical for all scatterers in the target. We name p and p as RAA and RAJ, respectively. Then, ϕ q (t m ) can be expressed as In (13), angular velocity ω represents the uniform rotation of the target, and RAA and RAJ represent the non-uniform rotation components. Herein, we select RAA and RAJ as the parameters indicating rotational nonuniformity. To obtain RAA and RAJ, we need to estimate coefficients of PPS.
When the target's maneuvering is not severe, the PPS model degrades into an LFM signal. In this case, only RAA is considered based on the estimation of linear and quadratic coefficients.
The linear and quadratic coefficients of the LFM signal are also referred to as centroid frequency and modulation rate, respectively. According to previous studies, fractional Fourier transform (FrFT) [30], Radon-Wigner transform [31], Lv's distribution [32], and centroid frequency-chirp rate distribution (CFCRD) [33] are effective approaches for parameter estimation of LFM signals. These methods transform a discrete signal into a two-dimensional parameter domain, and each LFM signal component becomes a two-dimensional peak function. The peak coordinates represent the centroid frequency and modulation rate, respectively. The FrFT algorithm is cross-term free for parameter estimation of multiple LFM signals, but costly for computation, especially for obtaining precise results. CFCRD suffers from cross-term interference under multiple LFM signal components, and yet is efficient in computation cost with fine precision. In this paper, CFCRD is preferred for parameter estimation of multiple LFM signals, and then, a rough value of RAA can be obtained.
When the target's maneuvering is severe, both RAA and RAJ need to be considered. In this case, coefficients of the CPF signal are estimated. According to previous studies, several algorithms for quadratic and cubic coefficient estimation can be retrieved, such as higher-order ambiguity function-integrated cubic phase function [34], scaled Fourier transform [35], keystone time-chirp rate distribution (KTCRD) [12], and chirp rate-quadratic chirp rate distribution [13]. Similarly, these algorithms transform a discrete signal into a two-dimensional parameter domain whose two dimensions represent quadratic and cubic coefficients, respectively. Considering the trade-off between calculation complexity and estimation precision, the KTCRD method is preferred for quadratic and cubic coefficient estimation of the CPF signal. For a rough estimation, the CFCRD method is still used for linear coefficient estimation of the CPF signal.
Owing to the possible residual error of translational motion compensation, noise interference, cross-term interference, etc., the estimated RAA and RAJ based on individual scatterers are inaccurate. To improve the accuracy, we adopt the averaging scheme. First, several range cells containing strong scatterers are selected for RAA and RAJ estimation. Second, in each range cell, the coefficients of several predominant PPS components are estimated. Then, rough values of RAA and RAJ are obtained. Third, based on the estimated RAAs and RAJs, average RAA and RAJ values are calculated. Through averaging, the estimation errors obtained from individual scatterers are alleviated.

New ISAR Imaging Approach based on Optimized Non-Uniform Rotation Transform
In this section, we discuss the implementation of the optimized non-uniform rotation transform and the corresponding ISAR imaging procedures.
With RAA and RAJ, high-order phase terms of the PPS can be constructed and then compensated uniformly in the Fourier transform matrix. Primarily, we concentrate on the compensation of the quadratic phase term in (13), which contains four parameters: x q , ω, p , t m . Among them, p and t m are known and x q ω indicates the Doppler frequency of scatter q. Through Fourier transform, x q ω is projected into cross-range cells in the discrete Doppler domain. Here, we introduce a new variable m to represent the discrete value of x q ω . Each m value corresponds to a row vector in the Fourier transform matrix. For a coherent processing interval containing M pulses, the Fourier transform matrix and m values are shown as follows: (14) where W = exp −j 2π M Note that m is proportional to the phase gradient of the row vector in the Fourier transform matrix. With m , a row vector in matrix P can be expressed as v row m, m = W mm , (15) where m = 0, 1, 2, . . . , M − 1 represents the discrete slow time t m . In fact, v row (m, m ) is the discrete version of the linear phase term in (13). Similarly, the discrete versions of quadratic and cubic phase terms in (13) can be expressed as where m − M/2 indicates that t m = 0 is located at the center of the pulse sequence. PRF is short for the pulse repeat frequency. In (16) and (17), PRF is induced by the discretization of t m . Following are the features of (16) and (17): (1) For all scatterers in the target, the values of m, p , and p are identical.
(2) For scatterers located in the same Doppler cell, their quadratic and cubic phase terms are approximately the same. Consequently, these phase terms can be compensated uniformly using a row vector. According to (16) and (17), two compensation matrices can be constructed with the variation of m and m , which are denoted as P q and P c , respectively. The size of P q and P c is the same as that of P. The Fourier transform matrix containing high order phase term compensation can be expressed as where represent Hadamard multiplication and P is the non-uniform rotation transform matrix. The imaging formula for targets with complex motion can be expressed as Given that RAA and RAJ are accurate, an optimal imaging result without blurring can be obtained through (19). However, the estimated RAA and RAJ based on parameter estimation of the PPS are not accurate enough, which degrades the imaging result. To solve this problem, we adopt the parameter searching scheme to obtain optimal RAA and RAJ and an optimal imaging result. First, we set the searching ranges of RAA and RAJ according to the roughly estimated values. Second, we set the searching step lengths according to the accuracy demands. Third, we carry out non-uniform rotation transform based on different RAA and RAJ values, and evaluate the corresponding imaging results using entropy or contrast function (Appendix A). A small entropy and a large contrast indicate good ISAR image quality. According to the imaging result with least entropy or largest contrast, optimal RAA and RAJ and an optimal ISAR image are obtained.
The parameter searching process mentioned above is a brute force solution. Though effective, it is costly in computation. Fortunately, optimal RAA and RAJ corresponding to the optimal imaging result are unique, and thus, the parameter searching process is a convex optimization problem that can be solved through a fast method such as gradient descent. We denote the entropy function of an ISAR image under different RAA and RAJ values as H(p , p ). Partial derivatives H(p ,p ) ∂p and H(p ,p ) ∂p are calculated and then used to update the values of p and p . When the absolute value of the entropy gradient is smaller than the preset tiny threshold, the entropy is regarded to reach the bottom with optimal RAA and RAJ and an optimal ISAR image.
In the new ISAR imaging approach proposed in this paper, the non-uniform rotation transform instead of conventional Fourier transform is used for cross-range pulse compression. High-order phase terms caused by the target's complex motion are compensated uniformly in the transform matrix according to (14)- (19). Consequently, the blurring of the ISAR image is eliminated and an optimal imaging result is obtained. Note that the non-uniform rotation transform is a generalized version of conventional Fourier transform. When both RAA and RAJ are zero, the non-uniform rotation transform degrades into conventional Fourier transform. Therefore, our new approach is feasible for ISAR imaging of targets with both smooth and complex motion. The procedure of the new ISAR imaging approach is summarized in Figure 2.

Experimental Results
In this section, the new imaging approach is validated on the echo pulses of a target with different complex motion. The imaging result of the RID algorithm is considered as comparison.
In the simulation, the parameter settings of the radar system are illustrated in Table 1. The target shown in Figure 3a is a plane model containing 42 scatterers. The backward scattering coefficients of all scatterers are set as 1. Let us assume that the radar is located at the origin of the radar coordinate system O-XYZ, where the O-XY plane is parallel to the horizontal plane. At the starting of the coherent processing interval, the target's coordinates in the radar coordinate system are (3000, 3000, 7000 m) and its velocity is v = (225, 300, 0)m/s. We set the acceleration of the target to zero. The ISAR image of the target obtained through range-Doppler algorithm is shown in Figure 3b. Obviously, the ISAR image is well focused and it is used as reference for the following experiments.

ISAR Imaging of Target with Constant Acceleration
We set the acceleration of the target as a = (90, 120, 0)m/s 2 , while the initial velocity and start position of target stay unchanged. The ISAR image obtained through the range-Doppler algorithm is shown in Figure 4a. Compared to Figure 3b, blurring occurs in the ISAR image and the image quality is poor, as indicated by the larger entropy and smaller contrast. The cross-range signal in range cell 485 is extracted and its time-frequency representation is illustrated in Figure 4b. The smoothed pseudo Wigner-Ville distribution [36] is adopted to obtain the time-frequency representation, where it is identified that two scatterers are located in this range cell and the Doppler frequencies of both are time-varying. The slope of the Doppler frequency line is induced by the target's acceleration. Moreover, the slope of the Doppler frequency line is related to its intercept: larger the intercept, larger is the slope. The relation between the slope and intercept stems from the rigid property of the target body. The linear Doppler frequency in Figure 4b indicates that the scatterers result in LFM signals, and only RAA needs to be considered. Signals in two range cells are selected, and their parameter distribution through the CFCRD method is shown in Figure 5. The calculation time for the two signals is 3.03 s. Four peaks can be identified in Figure 5a, and three peaks can be identified in Figure 5b. The estimated centroid frequency and modulation rate are illustrated in Table 2. The rough RAA of each scatterer is calculated and the average RAA is found to be 0.1262.  Next, the gradient descent method is used to obtain the optimal RAA and imaging result. The initial value of RAA is set as 0.1262. Through the iteration process, RAA is updated using p ← p − ∆·H(p )/∂p , where H(p )/∂p is the entropy gradient of the ISAR image and the step length ∆ is set as 0.001. When the entropy gradient is smaller than , the iteration stops. The RAA and ISAR image in the last iteration are output as the final results. The process and output of the gradient descent process are shown in Figure 6, which shows the entropy to descend gradually. After 51 iterations, the termination condition is met with a calculation time of 4.21 s. The output optimal RAA value is 0.1225. The final ISAR image based on the optimized non-uniform rotation transform is shown in Figure 6b. Compared to that in Figure 4a, the blurring is eliminated effectively and a high-quality ISAR image is obtained, as indicated by the smaller entropy and larger contrast.

ISAR Imaging of Target with Curve Flight Path
In this section, the complex motion of the target is set as the curve flight path, as shown in Figure 7a. Non-uniform rotation transform is performed to compensate for the adverse influence of the curve flight path. For simplicity, the attitude of the airplane model is kept constant during imaging. After pulse compression and translation motion compensation, the ISAR image based on the range-Doppler algorithm is obtained as shown in Figure 7b. It can be seen that the curve flight path leads to blurring of the ISAR image. The cross-range signals in range cells 485 and 508 are extracted and their time-frequency representations are illustrated in Figure 8. Two and four scatterers can be identified in these two range cells, respectively. Similar to Figure 4b, the Doppler frequencies in Figure 8 are time-varying. Moreover, the nonlinear feature of the Doppler curve in Figure 8 indicates that the changing rates of the Doppler frequency are not constant. Therefore, the cross-range signals satisfy the CPF model, and both RAA and RAJ need to be estimated.  The cross-range signal in range cell 485 is selected to estimate RAA and RAJ. To obtain the coefficients of the CPF signal, CFCRD and KTCRD are carried out as shown in Figure 9. According to Figure 9a, the centroid frequency and modulation rate of the two scatterers are (10, 2.745) and (24.12, 6.667), respectively. Thus, the rough RAA value is 0.27. According to Figure 9b, the modulation rate and quadratic modulation rate of the two scatterers are (2.5, 2) and (8, 10), respectively. Then, we obtain p 3 /p 2 ≈ 1. In other words, RAJ/RAA ≈ 1, so a rough RAJ value is also 0.27. To demonstrate the effect of RAA and RAJ compensation, two schemes are carried out. In the first scheme, only RAA compensation is considered in the gradient descent process, and the results are illustrated in Figure 10. The minimum ISAR image entropy is 5.8337 after 71 iterations. The calculation time is 5.84 s and the output RAA is 0.1443. Compared to Figure 7b, the ISAR image in Figure 10b shows a significant improvement. However, slight blurring can still be identified, as indicated by a red dash circle in Figure 10b. In the second scheme, both RAA and RAJ are considered and a two-dimensional gradient descent process is carried out. In each iteration, RAA is updated using p ← p − ∆· H(p ,p ) ∂p and RAJ is updated using p ← p − ∆· H(p ,p ) ∂p , where the step length ∆ is set as 0.002. Between two adjacent iterations, if the entropy change is smaller than 0.0001, the iteration process stops. The initial values of both RAA and RAJ are set as 0.27. The termination condition is met after 61 iterations, as shown in Figure 11a, and the calculation time is 16.43 s. The outputs of optimal RAA and RAJ are 0.1766 and 0.2075, respectively. The entropy of the output ISAR image shown in Figure 11b is smaller than that in Figure 10b. A larger contrast of Figure 11b also indicates that a better image quality is obtained.  We can conclude that for the target's complex motion on such a curve flight path, the cross-range signal satisfies the CPF signal model. Using a non-uniform rotation transform containing only RAA compensation, the blurring of the ISAR image is alleviated, but slight blurring still exists. When both RAA and RAJ are considered in the non-uniform rotation transform, an optimal ISAR image is obtained. Herein, we conclude that the proposed optimized non-uniform rotation transform is effective for ISAR imaging of targets with complex motion.

ISAR Imaging of Boeing 727 with Complex Motion
We also test the non-uniform rotation transform on the simulated data of Boeing B727 airplane provided by the U.S. Naval Research Laboratory [5]. The data contain 256 successive pulses, as shown in Figure 12a. The carrier frequency is 9 GHz, the bandwidth is 150 MHz, and PRF is 20 KHz. The total coherent processing time is 0.0128 s. Translational motion compensation has been accomplished. Because of the maneuvering of the target, the ISAR image obtained through the range-Doppler algorithm is blurred as shown in Figure 12b. To overcome the blurring, the optimized non-uniform rotation transform is carried out and its imaging result is compared with that of the RID algorithm. Primarily, rotational nonuniformity needs to be estimated. The time-frequency distribution of the cross-range signal in range cell 31 is obtained through the smoothed pseudo Wigner-Ville distribution, as shown in Figure 13a. The linear curves in the figure indicate that the signal satisfies the LFM model, thus only RAA needs to be estimated. Due to the short duration of the cross-range signal, CFCRD is not suitable for coefficient estimation in this case. Alternatively, FrFT [37] is used to estimate the signal coefficient as shown in Figure 13b. The centroid frequency and modulation rate of the two scatterers shown in Figure 13b are (−745, −24,546) and (1375.7, 98,304), respectively. Then, the rough RAA value is 26.2. Next, the gradient descent process is carried out to obtain the optimal RAA and imaging result. The step length in iteration is set as 1 and the termination condition is set such that the entropy gradient is smaller than 1 × 10 −6 . After 352 iterations, the termination condition is met and the calculation time is found to be 7.22 s; the output optimal RAA value is 32.24. The entropy change during the iteration and the output optimal ISAR image are shown in Figure 14. Compared to Figure 12b, the blurring in Figure 14b is eliminated significantly. The change in entropy and contrast also verifies the effectiveness of the proposed algorithm.  In contrast, an RID algorithm is used to obtain the instantaneous ISAR image of Boeing 727. In each range cell, the cross-range signal is processed with smoothed pseudo Wigner-Ville distribution to obtain the two-dimensional time-frequency distribution. At each time instant, the instantaneous Doppler frequencies of the cross-range signals in 64 range cells combine an RID ISAR image. There are 256 RID ISAR images in total, of which two representative images, 8th and 30th, are illustrated in Figure 15. Compared to the image shown in Figure 12b, blurring is eliminated significantly in the images shown in Figure 15. However, the defects are also obvious. First, Doppler spread of individual scatterers still exists, which can be identified in Figure 15a. Second, weak scatterers become almost invisible, which leads to the loss of target shape information. It is difficult to interpret an airplane from Figure 15b. Third, the RID image is unstable since images at different time instants are different. It is found that entropy and contrast are not proper norms for selecting the best one from 256 RID images. This leads to confusion and inconvenience in practice. The ISAR image observed in Figure 14b is still better than the RID images observed in Figure 15. Though entropy and contrast of the latter indicate better image quality, the former provides clearer and more stable target information. This confirms the advantages of non-uniform rotation transform over the RID algorithm. Table 3 compares three algorithms. In summary, the proposed method outperforms range-Doppler and RID algorithms in terms of the image quality. Moreover, the imaging result of the proposed method is unique and more stable compared to that of RID.

Conclusions
In this paper, a new imaging approach based on an optimized non-uniform rotation transform is proposed for ISAR imaging of targets with complex motion. Two parameters, RAA and RAJ, are defined to describe the equivalent rotational nonuniformity of the target. Based on RAA and RAJ, the quadratic and cubic phase terms resulting from the target's complex motion can be compensated uniformly in the Fourier transform matrix. A rough estimation of RAA and RAJ is realized through coefficient estimation of the PPS. The optimal RAA and RAJ are obtained through gradient descent optimization, after which an optimal ISAR image is obtained. In this way, ISAR imaging of the target with complex motion is converted into the parameter optimization problem. Through the optimized nonuniformity rotation transform, the blurring in the ISAR image of the target with complex motion is eliminated effectively. Compared to precedent imaging algorithms such as RID algorithm, the imaging result of the proposed method is clearer and more stable. Moreover, the additional computation cost of the proposed method is reasonable. Imaging results of targets with different complex motion and simulated data of Boeing 727 verify the effectiveness and advantages of the proposed algorithm. Our future work focuses on the possible adverse influence and countermeasures of translational motion compensation error on non-uniform rotation transform.  where E(·) represents the average operation.