A Novel 2-D Geometry Reconstruction Approach for Space Debris via Interpolation-Free Operation under Low SNR Conditions

Two unsolved key issues in inverse synthetic aperture radar (ISAR) imaging for non-cooperative rapidly spinning targets including high computational complexity and poor imaging performance in the case of low signal-to-noise ratio (SNR) are addressed in this work. In the strip-map imaging mode of SAR, it is well known that azimuth spatial invariant characteristics exist, and inspired by this, we propose a fast ISAR imaging approach for spinning targets. Our approach involves two steps. First, a precise analytic expression in the range-Doppler (RD) domain is produced using the principle of stationary phase (POSP). Second, a novel interpolation kernel function is designed to remove two-dimensional (2-D) spatial-variant phase errors, and the corresponding fast implementation steps that only require Fourier transform and multiplications are also presented. Finally, a well-focused ISAR image is obtained by compensating the azimuth high-order terms. Compared with current imaging methods, our approach avoids multi-dimensional search and interpolation operations and exploits the 2-D coherent integrated gain; the proposed method is of low computational cost and robustness in the low SNR condition. The effectiveness of the proposed approach is confirmed by numerically simulated experiments.


Introduction
The rapid developments of the aerospace industry raised the issue of space debris, which pose a serious threat to spacecraft or satellites in orbit, with potentially disastrous consequences [1][2][3]. In turn, detecting, identifying, and cleaning space debris is currently an important research topic. Among the different methods to monitor space debris, inverse synthetic aperture radar (SAR) (ISAR) represents a powerful tool able to produce high-resolution images of spinning targets, and it has been widely adopted to obtain space debris images [4][5][6]. However, the trajectories of space debris are complex, and besides translational motion, they are usually characterized by rapid spinning along the main rotating axis. This makes range-Doppler (RD) [7,8] and range instantaneous Doppler (RID) algorithms for ISAR imaging [9,10] unusable and, in turn, it makes it challenging to produce well-focused high-resolution ISAR imaging for the spinning target. One of the problems is that fast rotating targets cause great range cell migration (RCM) and non-stationary Doppler frequency modulation (DFM). Moreover, the RCM and DFM exhibit the two-dimensional (2-D) spatial-variant features [11][12][13][14][15][16][17][18], which are correlated to fast and slow times dynamics, and with the positions of scattering centers. If the 2-D spatial-variant dependent phases are not properly corrected, the final RD-based images would be blurred. Another issue, arising from the spinning motion of space debris and the long imaging distance, is the strong noise, which makes the signal-to-noise ratio (SNR) of the echo signal very low and prevents a precise estimation of the motion parameters [11,17]. The overall effect is that the performance of the conventional ISAR imaging for spinning target is seriously degraded.
To overcome the above mentioned issues and produce images of rapidly spinning targets, several ISAR approaches with narrowband and wideband transmitted signals have been suggested and developed [11][12][13][14][15][16][17][18][19][20][21]. In [11] a single-range Doppler interferometry approach (SRDI) for spinning targets was developed building on the properties of the Doppler spectrum in one rotation period. SRDI detects the spinning targets by transforming the 1-D Doppler signal into a 2-D time-frequency signal, from which the energy of spinning target may be accumulated into a sharp peak. However, its performance is degraded from cross-terms, interferences, and resolution loss [11,12]. Moreover, its complexity is high, due to the need for a multi-dimensional search. To overcome the shortcomings of SRDI, a single-range matching filtering method (SRMF) was developed [13,14] to produce well-focused an ISAR image. To this aim, a match function was designed with different rotating radii and initial phases to compensate the transverse echo signal cell-by-cell. The clean technique was also exploited to migrate the high sidelobes effects. In [15], the authors developed a single-range image fusion method to produce 2-D images in the Cartesian coordinates. The polar format algorithm (PFA) [16], the modified complex valued back-projection algorithm (BPA) [17], and the tomography algorithm [18] were also used to obtain the images of rotating targets. However, all those techniques have high computational complexities due to the involved integrals, and still lack fast implementation. In [19], a real-valued generalized Radon transform (GRT) approach was presented to achieve imaging of spinning targets, which employs the sinusoidal envelop of the echo signal for spinning targets to obtain an ISAR image via non-coherent accumulation operation. As such, it cannot obtain satisfactory imaging performance in low SNR conditions. In [20], an extended Hough transform (HT) was introduced to realize ISAR imaging of rotating targets. However, the algorithm is time-consuming, due to the multi-dimensional parameter search. Upon using the sparse property of the echo signal and the matching pursuit (MP) algorithm, a micro-motion imaging approach was proposed in [21], which, however, has a high computational cost and fails when Doppler ambiguity is present. In [22], well-focused ISAR images of rapidly spinning targets were generated using the segmental pseudo Keystone transform (KT) (SPKT). However, this technique is also time-consuming, and it suffers from accuracy loss due to the need for an initial phase estimation stage and to the interpolation operation used for each scattering point of the spinning target. Overall, despite all the aforementioned algorithms having their merits, a globally effective ISAR imaging method for non-cooperative spinning targets is still lacking, especially for those situations where the SNRs are low.
In this paper, exploiting the azimuth spatial invariance of strip-map SAR imaging mode [23,24], an effective imaging approach for rapidly spinning targets is suggested, which works also with a low SNR environment. To this aim, we first formulate the RCM and time-varying DFM caused by the target's spinning motion in terms of range and azimuth, i.e., a 2-D structure. Then, using the principle of stationary phase (POSP) [25,26], a precise analytic expression in RD domain is obtained, where the scatters located in one range cell with different azimuth positions are assigned to the same range migration trajectory to effectively eliminate the azimuth spatial variance of the initial phase. Third, a novel interpolation kernel function is designed to remove the range variance of the RCM and DFM, and an interpolation-free mapping algorithm is suggested to further reduce the computational burden. Finally, by compensating the azimuth high-order terms, a high resolution ISAR image for rapidly spinning targets is obtained. Compared to currently available imaging methods, our technique can be efficiently implemented, since it only requires Fourier transform (FT) without the need of any interpolation or multi-dimensional search operations. Furthermore, thanks to the use of a 2-D coherent integrated gain, a satisfactory imaging performance in low SNR condition may also be obtained. We also performed numerically simulated experiments to verify the correctness and the effectiveness of the proposed algorithm.
The paper is organized as follows. In Section 2, the geometry and signal model for rapidly spinning targets are provided. In Section 3, a novel efficient ISAR imaging algorithm in the RD domain using the SAR technique is described. The fast implementation of the method and some remarks about its practical applications are presented in Section 4. In Section 5, results from simulated experiments are given together with their analysis. In Section 6, we discuss the computational complexity of the method, whereas Section 7 closes the paper with some concluding remarks.

Imaging Model of Spinning Target and Problem Formation
The geometric configuration of ISAR imaging for spinning space debris is depicted in Figure 1, where the origin O denotes the position of the center of the rotating target. During the imaging period, the motion of the target may be decomposed into translation and rotation. The translational motion is parallel to the line-of-sight (LOS) and does not contribute to the cross-range resolution. On the other hand, rotation changes the Doppler frequency and affects the azimuth imaging. In Figure 1, we denote by P an arbitrary scattering point that rotates around the origin O with constant rotation rate ω. The spinning angular rate ω is the sum of the self-spinning angular velocity and a constant term due to the time-variant LOS angle [18][19][20][21][22]. The distance between the spinning target center and the radar is denoted by R 0 R 0 r p , where ϕ p is the initial phase of the spinning scatterer P.
Remote Sens. 2020, 05, x FOR PEER REVIEW 3 of 18 without the need of any interpolation or multi-dimensional search operations. Furthermore, thanks to the use of a 2-D coherent integrated gain, a satisfactory imaging performance in low SNR condition may also be obtained. We also performed numerically simulated experiments to verify the correctness and the effectiveness of the proposed algorithm. The paper is organized as follows. In Section 2, the geometry and signal model for rapidly spinning targets are provided. In Section 3, a novel efficient ISAR imaging algorithm in the RD domain using the SAR technique is described. The fast implementation of the method and some remarks about its practical applications are presented in Section 4. In Section 5, results from simulated experiments are given together with their analysis. In Section 6, we discuss the computational complexity of the method, whereas Section 7 closes the paper with some concluding remarks.

Imaging Model of Spinning Target and Problem Formation
The geometric configuration of ISAR imaging for spinning space debris is depicted in Figure 1, where the origin denotes the position of the center of the rotating target. During the imaging period, the motion of the target may be decomposed into translation and rotation. The translational motion is parallel to the line-of-sight (LOS) and does not contribute to the cross-range resolution. On the other hand, rotation changes the Doppler frequency and affects the azimuth imaging. In Figure  1, we denote by P an arbitrary scattering point that rotates around the origin with constant rotation rate . The spinning angular rate is the sum of the self-spinning angular velocity and a constant term due to the time-variant LOS angle [18][19][20][21][22]. The distance between the spinning target center and the radar is denoted by ≫ , where is the initial phase of the spinning scatterer . During the rotation, shadowing of the scatterers will cause the loss of the echo signal of the spinning target. The non-uniform rotating rate, or residual translational movement, will result in defocused ISAR images. However, the current paper is focused on designing an efficient imaging algorithm based on the presented imaging model. Here we assume that translational motion compensation (TMC) has been implemented by the available TMC algorithm [27,28] and focus on the effects of target rotation. According to Figure 1, the instantaneous slant range ; between the radar phase center and the scatterer in is given by where represents the azimuth slow-time variable. Assuming that the radar transmits a linear frequency-modulated (LFM) signal, the received baseband echo signal of the point target after demodulation is During the rotation, shadowing of the scatterers will cause the loss of the echo signal of the spinning target. The non-uniform rotating rate, or residual translational movement, will result in defocused ISAR images. However, the current paper is focused on designing an efficient imaging algorithm based on the presented imaging model. Here we assume that translational motion compensation (TMC) has been implemented by the available TMC algorithm [27,28] and focus on the effects of target rotation. According to Figure 1, the instantaneous slant range r t a ; r p between the radar phase center and the scatterer in P is given by r t a ; r p = R 0 + r p sin ωt a + ϕ p (1) where t a represents the azimuth slow-time variable. Assuming that the radar transmits a linear frequency-modulated (LFM) signal, the received baseband echo signal of the point target after demodulation is where σ p is the reflection coefficient of the scatterer in P, t r is the range fast-time variable, and ω r (·) and ω a (·) represent the envelope in the range and azimuth direction, respectively. c, f c , and γ denote the speed of light, the carrier frequency of radar signal, and the chirp rate, respectively. Performing FT with respect to the range fast-time variable t r and removing the range frequency-modulated term, the received echo signal becomes where λ denotes the carrier wavelength, and f r is the range frequency variable with respect to t r . In Equation (3), the phase consists of three terms. The first is the linear part corresponding to the range direction position, the second term is the RCM, and the third one is due to the azimuth time-varying DFM. As a matter of fact, the RCM and DFM terms contain a 2-D spatial dependence, which is due to rotation of the target and makes it difficult to obtain an effective imaging scheme for spinning targets. Indeed, unless the second and third terms are satisfactorily compensated, the final ISAR image will be blurred. Notice that the 2-D spatially dependent phase terms are different for each scatterer, and therefore, the compensation should be performed pixel by pixel. In order to face this problem while utilizing the sinusoidal envelope of the spinning targets in the range slow-time domain, existing GRT methods correct the phase via multi-dimensional searching operation along the rotating radius r p and initial phase ϕ p for all scatterers in the spinning target [19]. Performance is, however, poor due to the non-coherent accumulation. In the SPKT method, the initial phase estimation and the interpolation operation are required for each scatterer in the spinning target [22], thus leading to a great computational burden and a low imaging quality in the environment, inducing a low SNR. Here, by exploiting the inherent azimuth invariance of SAR, a fast ISAR imaging method in the case of spinning targets is developed to produce a well-focused ISAR image for spinning targets.
The instantaneous Doppler frequency f r is calculated by taking the derivative of the phase term in (3) with respect to t a as From the Nyquist sampling theorem, the pulse repetition frequency (PRF) must satisfy PRF ≥ 2 f dmax . Thus, the maximum of Doppler frequency is f dmax = 2ωr pmax /λ, where r pmax is the longest rotation radius of the targets. The Doppler bandwidth of scatterers is then

Proposed Imaging Approach Description
In this section, we focus on the 2-D spatial-variant RCM and DFM fast correction under low SNR conditions. A fast imaging approach for non-cooperative spinning targets is developed. Before going into the details of the imaging algorithm, let us summarize the principles of azimuth invariant SAR strip-map imaging mode. A schematic representation is given in Figure 2; panel 2(a) shows the geometric model of multiple targets located in different range and azimuth cells at a given time, whereas panels 2(b) and 2(c) show the situation after range compression. For targets located in the same range cell at different azimuth times, the migration trajectories are overlapped in the RD domain, and this is the well-known azimuth invariance of SAR strip-map imaging mode [23][24][25][26].
SAR strip-map imaging mode. A schematic representation is given in Figure 2; panel 2(a) shows the geometric model of multiple targets located in different range and azimuth cells at a given time, whereas panels 2(b) and 2(c) show the situation after range compression. For targets located in the same range cell at different azimuth times, the migration trajectories are overlapped in the RD domain, and this is the well-known azimuth invariance of SAR strip-map imaging mode [23][24][25][26].
In our proposed method, in the Doppler domain and based on the above-mentioned invariance, echo signals of the same distance are only processed once, thus reducing the computational complexity.

Analytical Expression Derivation in RD Domain
Scatterers with the same rotating radius have overlapped trajectories in the range Doppler domain, and thus the RCM correction located in same range gate with different azimuth time instant can be realized just once, thus improving imaging efficiency. In order to obtain an expression in the range Doppler domain after range compression, let us perform FT of the slow-time variable in Equation (2), thus arriving at where is the transmitted LFM signal bandwidth, and denotes the range envelope. The integral in (6) cannot be evaluated analytically due to the sin function term in the slant range formula ; , and to obtain an analytic expression we resort to POSP as in the following. The phase expression in (6) is According to the POSP, the value in (7) is calculated by the stationary point, that is, By solving the Equation (8), one obtains the stationary point * = 1 arccos − 2 − , In our proposed method, in the Doppler domain and based on the above-mentioned invariance, echo signals of the same distance are only processed once, thus reducing the computational complexity.

Analytical Expression Derivation in RD Domain
Scatterers with the same rotating radius have overlapped trajectories in the range Doppler domain, and thus the RCM correction located in same range gate with different azimuth time instant can be realized just once, thus improving imaging efficiency. In order to obtain an expression in the range Doppler domain after range compression, let us perform FT of the slow-time variable t a in Equation (2), thus arriving at where B r is the transmitted LFM signal bandwidth, and sinc denotes the range envelope. The integral in (6) cannot be evaluated analytically due to the sin function term in the slant range formula r t a ; r p , and to obtain an analytic expression we resort to POSP as in the following. The phase expression in (6) is According to the POSP, the value in (7) is calculated by the stationary point, that is, By solving the Equation (8), one obtains the stationary point Remote Sens. 2020, 12, 2059 6 of 17 which is then substituted in Equation (1), leading to the instantaneous slant distance expression in the RD domain Equation (10) reveals that r( f a ) includes a trigonometric function sin, which is generated by the fast rotation of the targets. Finally, the expressions of the echo signal in the range-Doppler in Equation (11) is obtained, where the second exponential term contains the initial phase ϕ p , and the first and the third exponentials represent high order phase terms.

The Proposed ISAR Imaging Method
In Equation (11), in the high order phase term or the sin c envelop, the initial phase ϕ p is absent. The range sin c envelop is, in turn, only related to the rotating radius, showing that scatterers in the same range cell and different azimuth instantaneous times have the same trajectories in the RD domain. Let us now describe in some detail the steps needed to derive the proposed ISAR imaging method.
In Equation (11), the range envelope term sin c[·] depends on the radial range r( f a ) and the azimuth frequency f a , which cause the RCM. Because of the RCM, the energy of the same target is scattered in different range cells in the original (r − f a ) plane. In order to efficiently remove the RCM, we introduce a novel interpolation kernel function [29,30], given by After applying the interpolation, the signal in the RD domain is given in (13), where RCM has been removed. The phenomenon is illustrated in Figure 3. The original discrete signal in the (r − f a ) plane is shown in Figure 3a, where the sample points of the Doppler frequency are equally spaced. The signal after traditional KT in the (r − f a ) plane is illustrated in Figure 3b, where Doppler frequency intervals are varied with range sampling points, but the intervals between range sampling points are identical. The signal of the proposed method in the (r − f a ) plane is demonstrated in Figure 3c, where different Doppler frequencies have different range sampling intervals.
Upon comparison of the original RD domain expression in (11) to the RD domain expression after transformation in (13), we see that along the sinc trajectory, the energies of the scatterers are transformed into one range cell , with the sinc • term depending only on the rotating radius . The azimuth frequency has been removed from the range envelope, which means that the RCM has been properly corrected (i.e., compensated).
Let us now address the azimuth processing procedure. The compensation function , is designed to compensate the high order phase terms in (13) as Multiplying (13) by (14), the signal is given by such that the energy of a single scatterer is still in the same range cell . Finally, by applying an azimuth inverse fast Fourier transform (IFFT), the ISAR image expression becomes From (16), one obtains the coordinates , of the scatterer in and, in turn, an accurate imaging of the target.

Fast Implementation.
In general, the novel interpolation operation in the RD domain proposed in Section 3 may increase the computational complexity and also introduce loss of accuracy. On the other hand, the data for non-sampling points may be obtained through a zero-padding in another domain, where the interpolation operation is replaced by FFTs and phase multiplication [31][32][33]. As a consequence, Upon comparison of the original RD domain expression in (11) to the RD domain expression after transformation in (13), we see that along the sin c trajectory, the energies of the scatterers are transformed into one range cell r p , with the sin c[·] term depending only on the rotating radius r p .
The azimuth frequency f a has been removed from the range envelope, which means that the RCM has been properly corrected (i.e., compensated).
Let us now address the azimuth processing procedure. The compensation function H(r, f a ) is designed to compensate the high order phase terms in (13) as Multiplying (13) by (14), the signal is given by such that the energy of a single scatterer P is still in the same range cell r p . Finally, by applying an azimuth inverse fast Fourier transform (IFFT), the ISAR image expression becomes S(r, t a ) = σ p sin c 2πγT p c r − r p sin c πB a ω ωt a + ϕ p From (16), one obtains the coordinates r p, ϕ p of the scatterer in P and, in turn, an accurate imaging of the target.

Fast Implementation
In general, the novel interpolation operation in the RD domain proposed in Section 3 may increase the computational complexity and also introduce loss of accuracy. On the other hand, the data for non-sampling points may be obtained through a zero-padding in another domain, where the interpolation operation is replaced by FFTs and phase multiplication [31][32][33]. As a consequence, the drawback mentioned above may be eliminated. In fact, the interpolation function may be reformulated as ∆r( f a ) = r p ·α f a r p (17) where the nonlinear shifting coefficient α f a r p is given by In order to obtain a fast implementation of the interpolation-free mapping method, the shift length ∆r( f a ) for the azimuth bin is decomposed into two terms, given by where n i is an integer, δ i ∈ (0, 1), and ∆r is the range cell size. Then, the interpolation operation may be accomplished through the following steps: (1) According to (7), one extracts the range cell data for each azimuth bin in the RD domain, i.e., S r = S r, f i a (20) and calculates n i and δ i according to (10) and (12); (2) Zero-padding S i r to S r of length 4 to achieve and indexing circularly in [1, · · · , N], where N is original data length; (3) FFT of S r to obtain S f r ; (4) Phase multiplication where f r is the new range frequency; (5) IFFT of S f r to complete the sub-bin shift of S r ; where F is the interpolation mapped data; (7) Iterating step (1) to step (6) for all bins to achieve interpolation-free mapping of the full image.
The interpolation-free mapping method replaces the RD domain interpolation operation through zero-padding in the wavenumber domain, which reduces the computational complexity and ISAR image errors. The flowchart of the algorithm and its fast implementation is shown in Figure 4.

Some Remarks about the Practical Application of the Proposed Method
(1) Range and Azimuth Resolution Analysis: The range resolution of ISAR imaging is determined by the transmitted signal bandwidth = 2 (25) whereas the azimuth resolution depends on the azimuth Doppler frequency resolution where ∆ is Doppler bandwidth, i.e., the inverse of the coherent imaging accumulation time ∆ = 1 ⁄ . The azimuth resolution may be thus rewritten as where ∆ is the rotation angle of the spinning target relative to the radar platform during the coherent accumulation time. In summary, the range resolution of ISAR is determined by the radar transmitting signal bandwidth, whereas the azimuth resolution is obtained by the target's rotation angle relative to the radar's line of sight during the coherent imaging accumulation.
(2) Effects of Estimation of on ISAR Images: As mentioned before, the proposed method is based on the knowledge of the angular velocity . In practice, however, the estimation of the angular velocity may not completely accurate.

Some Remarks about the Practical Application of the Proposed Method
(1) Range and Azimuth Resolution Analysis: The range resolution ρ r of ISAR imaging is determined by the transmitted signal bandwidth B r whereas the azimuth resolution ρ a depends on the azimuth Doppler frequency resolution where ∆ f d is Doppler bandwidth, i.e., the inverse of the coherent imaging accumulation time ∆ f d = 1/T. The azimuth resolution may be thus rewritten as where ∆θ is the rotation angle of the spinning target relative to the radar platform during the coherent accumulation time. In summary, the range resolution of ISAR is determined by the radar transmitting signal bandwidth, whereas the azimuth resolution is obtained by the target's rotation angle relative to the radar's line of sight during the coherent imaging accumulation.
(2) Effects of Estimation of ω on ISAR Images: As mentioned before, the proposed method is based on the knowledge of the angular velocity ω. In practice, however, the estimation of the angular velocity may not completely accurate. In order to assess the effects of this uncertainty on the final ISAR images, let us remember that the imaging results are obtained by coherent accumulation through the FT, which means the final phase error ∆ψ should be restricted in ∆ψ < π 2 . Assuming that the uncertainty of ω is δω, the phase error of the spinning target during the coherent time ∆t(0 < ∆t < π/ω) may be expressed by where r max is maximum rotating radius for the spinning target. Assuming δω∆t/2 1, we obtain sin(δω∆t/2) ≈ δω∆t/2 and sin ω∆t + ϕ p + δω∆t/2 ≤ 1. The constrained condition in (28) can be written as δω ω π < λ 2r max (29) Equation (29) provides a sufficient and unnecessary condition for obtaining well-focused ISAR images for the spinning target and should be taken into account in selecting the system parameters.
(3) Removing the effects of the sidelobes: It may happen that the sidelobes of a strong target hide the mainlobe of a weak target. Therefore, clean technology [34] is employed to suppress the sidelobes. The process of clean method proceeds as follows: (a) Estimate the initial phase ϕ p and rotation radius r p of all scatterers. The foregoing imaging algorithm may have already estimated these parameters. (b) According to the estimated initial phase ϕ p and rotation radius r p , a point spread function (PSF) is designed. Assuming that the coordinate of the maximum peak is r pξ , ϕ pξ , the PSF is given by (c) Using a minimum norm criterion, the scattering coefficient at r pξ , ϕ pξ is estimated by whereσ is the estimated value of the scattering coefficient. Performing the derivative of I(σ) with respect toσ leads to where X * (r, ψ) represents the conjugate function of X(r, ψ). When ∂I/∂σ = 0, the reflection coefficient of scatterer in coordinate r pξ , ϕ pξ is obtained. Under this condition, the estimated value of the scattering coefficient isσ Using this technique, the resulting r pξ , ϕ pξ ,σ are accurate estimates of the coordinates and the reflection coefficient. Subtract the contribution of scatterer r pξ , ϕ pξ from the original signal, thus arriving at s clean (r, ψ) = s −σX (35) (e) Use s clean (r, ψ) as the new input to step (a), and repeat until convergence.

Simulation Results and Discussion
In this section, results from numerically simulated experiments are provided to confirm the effectiveness and the robustness of the proposed ISAR imaging algorithm for rapidly spinning targets. The simulation settings are given in Table 1.

Imaging Results for a Single Scatterer Target
Results for the case of a single scatterer are shown in Figure 4. The rotation radius of the target was 3 m, and the initial phase was 0 rad. The distance between the rotating center and the antenna phase center was set to 20 km. Figure 5a shows the scatterer distribution in Cartesian coordinates. Figure 5b provides the 2-D echo signal data. The profile after range compressed for echo signal is shown in Figure 5c, where a trigonometric function may be clearly recognized. The envelop of signal in the RD domain is shown in Figure 5d. Obviously, the energy of a single spinning target in the RD domain spanned multiple range cells. We now separated the energies in the RD domain into ψ ∈ [0, π] and ψ ∈ [−π, 0] parts and performed an interpolation-free mapping operation according to (12) for each part. The results are depicted in Figures 5e and 5f, respectively. It can be seen that the main energy of a single target was concentrated into one range cell, which shows that the RCM was well corrected using the proposed algorithm. The ISAR image after high order phase compensation and combining the ψ ∈ [0, π] and ψ ∈ [−π, 0] parts of energy are presented in Figure 5g. The imaging result after clean post-processing is shown in Figure 5h. Finally, the interpolation-free mapping is shown in Figure 5i.

Imaging Results for Surface Targets
The simulation results for a set of targets on a surface are depicted in Figure 6. The surface was composed of eight scatterers, whose rotation radii were 2 , 2 , 2 , 2 , 3 , 3 , 3 , 3 , respectively, and the initial phases were 0, 2 ⁄ , − 2 ⁄ , , 4 ⁄ , − 4 ⁄ , 3 4, −3 4 ⁄ ⁄ , respectively. The reflection coefficients were set to one, and the distance between the antenna phase center and the rotating center was set to 20 . Figure 6a illustrates the scatterer distribution under the Cartesian coordinates. The profile after compressing range for echo signal is shown in Figure 6b, which resembles a set of trigonometric functions with all the trajectories for the scattering points intersecting with each other. In the GRT algorithm, the energy of the target is accumulated along the profile, as shown in Figure 6b. Thus, searching different initial phases and rotating radii was needed here. The ISAR image after high order phase compensation is shown in Figure 6c, and after clean processing in Figure 6e. Meanwhile,

Imaging Results for Surface Targets
The simulation results for a set of targets on a surface are depicted in Figure 6. The surface was composed of eight scatterers, whose rotation radii were 2m, 2m, 2m, 2m, 3m, 3m, 3m, 3m, respectively, and the initial phases were 0, π/2, −π/2, π, π/4, −π/4, 3π/4, −3π/4, respectively. The reflection coefficients were set to one, and the distance between the antenna phase center and the rotating center was set to 20 km. Figure 6a illustrates the scatterer distribution under the Cartesian coordinates. The profile after compressing range for echo signal is shown in Figure 6b, which resembles a set of trigonometric functions with all the trajectories for the scattering points intersecting with each other. In the GRT algorithm, the energy of the target is accumulated along the profile, as shown in Figure 6b. Thus, searching different initial phases and rotating radii was needed here. The ISAR image after high order phase compensation is shown in Figure 6c, and after clean processing in Figure 6e. Meanwhile, to demonstrate the effectiveness of the proposed algorithm, the comparisons with the GRT imaging result and GRT-clean imaging result are presented in Figures 6d and 6f, respectively.
As it is apparent from the figures, both of the images obtained after clean processing may be used to reconstruct the structure of the surface target. However, the proposed method is computationally efficient since the 2-D search is not needed here. to demonstrate the effectiveness of the proposed algorithm, the comparisons with the GRT imaging result and GRT-clean imaging result are presented in Figure 6d and Figure 6f, respectively. As it is apparent from the figures, both of the images obtained after clean processing may be used to reconstruct the structure of the surface target. However, the proposed method is computationally efficient since the 2-D search is not needed here.

Imaging Results under Different SNRs Condition
To further show the robustness of different methods, Gaussian white noise was added. Figure 7 shows the imaging results obtained by using the GRT method of [19], the SPKT method of [22], and our proposed method for two values of the SNR: 0 dB and -35 dB, respectively. All the GRT, SPKT, and the proposed method have coherent integration in the range direction.

Imaging Results under Different SNRs Condition
To further show the robustness of different methods, Gaussian white noise was added. Figure 7 shows the imaging results obtained by using the GRT method of [19], the SPKT method of [22], and our proposed method for two values of the SNR: 0 dB and -35 dB, respectively. All the GRT, SPKT, and the proposed method have coherent integration in the range direction. However, thanks to the coherent integration in the azimuth direction, our method and the SPKT showed a higher SNR gain than the GRT method due to the incoherent integration in the azimuth direction. The imaging results under different SNRs were consistent with the output SNR gain analysis, which verifies the correctness and effectiveness of the proposed method. However, compared with the SPKT method [22], our proposed method has great advantages in computational cost, since it avoids the initial phase estimation and interpolation operation used for each scatterer of the spinning target.

Computational Complexity
The running times of GRT, SPKT, and the proposed algorithm are plotted versus imaging scene sizes. The results are illustrated in Figure 8, where MATLAB running on an Intel quad-core processor, CPU clocked frequency at 3.2 GHz, memory 8 GB, and Windows 10 were used. As it is apparent from Figure 8, our method was faster than GRT and SPKT methods. Moreover, as the size increased, the difference became more striking. However, thanks to the coherent integration in the azimuth direction, our method and the SPKT showed a higher SNR gain than the GRT method due to the incoherent integration in the azimuth direction. The imaging results under different SNRs were consistent with the output SNR gain analysis, which verifies the correctness and effectiveness of the proposed method. However, compared with the SPKT method [22], our proposed method has great advantages in computational cost, since it avoids the initial phase estimation and interpolation operation used for each scatterer of the spinning target.

Computational Complexity
The running times of GRT, SPKT, and the proposed algorithm are plotted versus imaging scene sizes. The results are illustrated in Figure 8, where MATLAB running on an Intel quad-core processor, CPU clocked frequency at 3.2 GHz, memory 8 GB, and Windows 10 were used. As it is apparent from Figure 8, our method was faster than GRT and SPKT methods. Moreover, as the size increased, the difference became more striking.

Discussion
In this section, we compare our method with the existing GRT [19] and SPKT [22] methods in terms of their computational complexity and denote by and the number of samples in the range and the azimuth dimension, respectively. In GRT, the computational complexity of range compression consisting of -times -point FT and complex multiplications is log + . Let and denote, respectively, the search numbers for rotating radius and initial phase. The computational complexity of energy accumulation is . Therefore, the overall

Discussion
In this section, we compare our method with the existing GRT [19] and SPKT [22] methods in terms of their computational complexity and denote by N r and N a the number of samples in the range and the azimuth dimension, respectively. In GRT, the computational complexity of range compression consisting of N a -times N r -point FT and complex multiplications is O N a N r log 2 (N r ) + N a N r . Let Q and M denote, respectively, the search numbers for rotating radius and initial phase. The computational complexity of energy accumulation is O(MQN a ). Therefore, the overall computation complexity of GRT is C GRT = O N a N r log 2 (N r ) + N a N r + MQN a The SPKT method is based on the segmental pseudo Keystone transform, which is designed to realize migration through resolution cell correction. The whole data is divided into P segments (P ≥ 1, and N a /P is integer) and the number of scatterers is H. By ignoring the low computational cost of a few steps, the computational cost of the SPKT method is approximately given by In our proposed algorithm, implementation mainly includes range compression, N r -times azimuth FFT, N a N r point interpolation operation, complex multiplication, and azimuth IFFT. Among them, the range compression step consists of a range FFT with complexity O N a N r log 2 (N r ) , a complex multiplication with complexity O(N a N r ), and a range IFFT with complexity O N a N r log 2 (N r ) , then followed by N r -times azimuth FFT with a computational complexity of O N r N a log 2 (N a ) . The computational complexity of the interpolation operation is 2(2N ker − 1)N a N r , where N ker is the length of the interpolation kernel. This step could be time consuming, and therefore, the interpolation-free mapping method is considered here. The interpolation-free mapping method involves zero-padding, FFT, phase multiplication, and IFFT. Since only a few zeros are added to the data, the total computational cost of this step is O 2N a N r log 2 (N r ) + N a N r . After this, the high order phase is compensated through a complex multiplication with complexity O(N a N r ). Finally, N r -times azimuth IFFT with complexity O N r N a log 2 (N a ) is conducted to obtain the ISAR image. Therefore, the overall computation complexity of the proposed method is C proposed = O 4N a N r log 2 (N r ) + 2N r N a log 2 (N a ) + 3N a N r (38) The above analysis confirms that, in general, our proposed algorithm is convenient in terms of computational complexity. In addition, compared to GRT and SPKT methods, the proposed algorithm does not need to estimate the initial phase.

Conclusions
Rapidly spinning space debris generates echo signals with great RCM during a short coherent processing interval. Upon exploiting the inherent azimuth invariance of SAR, we suggested and discussed in detail an efficient ISAR imaging algorithm. At first, a precise analytic expression in the RD domain was developed based on the POSP principle. Then, a novel interpolation kernel function was designed to remove 2-D space dependent phase errors. At the same time, an efficient interpolation-free mapping algorithm was exploited to reduce the computational complexity, and the energy spanning multiple range cells were transformed into one range cell. Finally, a well-focused ISAR image was produced by compensating the azimuth high-order modulating terms. Simulation results prove that the proposed ISAR imaging algorithm is robust, effective, and time efficient compared to GRT and SPKT methods.