An Inverse Synthetic Aperture Ladar Imaging Algorithm of Maneuvering Target Based on Integral Cubic Phase Function-Fractional Fourier Transform

Yakun Lv 1,* ID , Yanhong Wu 1, Hongyan Wang 2, Lei Qiu 1, Jiawei Jiang 1 and Yang Sun 3 ID 1 Department of Electronic and Optical Engineering, Space Engineering University, Beijing 101416, China; mail2wyh@163.com (Y.W.); qiulei2013@nudt.edu.cn (L.Q.); herojjw@163.com (J.J.) 2 School of Space Information, Space Engineering University, Beijing 101416, China; yhgnaw@163.com 3 Science and Technology on Complex Electronic System Simulation Laboratory, Beijing 101416, China; fireflypd@buaa.edu.cn * Correspondence: qiul06@mails.tsinghua.edu.cn; Tel.: +86-010-6636-5117


Introduction
With the development of radar detection technology, high-precision target imaging has become an important aspect of the detection task.Inverse synthetic aperture ladar (ISAL) combines coherent laser technology and inverse synthetic aperture technology, overcoming the limitations of the actual aperture and diffraction.ISAL also overcomes the shortfalls of traditional microwave imaging radars that cannot provide enough range resolution for remote target and small target imaging and solves the problem experienced by traditional laser imaging radar, which cannot perform the high-resolution imaging of a moving target [1].ISAL is the only optical means by which centimeter-level resolution can be obtained at a range of thousands of kilometers [2].Therefore, ISAL imaging can fulfil the requirement for high precision imaging and quasi real-time imaging for target surveillance.
ISAL imaging is similar to the traditional inverse synthetic aperture radar (ISAR) imaging principle but due to the use of laser as a radiation source, ISAL has an ultra-high carrier frequency, ultra-large bandwidth and extremely short wavelength.Compared with ISAR, ISAL has higher resolution, smaller imaging angle and shorter imaging time [3].Research on ISAL has mainly focused on principal analysis and algorithm simulation.Some close-range field tests have been reported [4][5][6][7][8].

ISAL Signal Echo Model
The three-dimensional (3D) imaging geometry of the maneuvering target is shown in Figure 1a.Where the coordinate origin O is the target turntable center, point P(x p , y p , z p ) is any scattering point on the target and r p is the scattering point P position vector starting from O. ω is the rotational angular velocity vector of the target three-dimensional motion.The ISAL imaging projection plane Γ is determined by the vector ω and the radar line-of-sight direction (LOS) unit vector R, ω can be decomposed into the radial rotational component, ω r along the LOS and the rotational component, ω e perpendicular to the LOS.ω r cannot cause the radial movement of the scattering point, that is, it will not cause the phase change of the echo and ω e will cause the scattering point to move radially, resulting in Doppler frequency variation, which can achieve high-resolution ISAL imaging of the target, ω e is called effective rotation component.The three-dimensional motion velocity of the target can be decomposed into a component v in the Γ plane and a component perpendicular to the Γ plane and the vertical component does not affect the imaging of the target, so this component can be ignored.For the parallel component v, it can be decomposed into the radial component v r along the LOS and the component v e perpendicular to the LOS.v r causes the Doppler shift of the target echo which cause phase change, while v e does not generate Doppler shift.
After the above analysis, the effective component in the three-dimensional (3D) imaging geometry can be projected onto the imaging plane Γ to obtain a two-dimensional (2D) turntable imaging geometry as shown in Figure 1b [18].In two-dimensional (2D) imaging geometry, only the relative motion between radar and target is considered.Where O is the reference point and XOY is the rectangular coordinate system fixed on the target, the Y axis is the direction of the radar LOS.The target moves along the Y axis with the speed of v r and rotates around the O point at the angular velocity of ω e .Suppose that at time t (t is full-time and satisfies the equation t = t k + t m , where t k is the range fast time and t m is the azimuth slow time, m = 1, 2, . . ., M), the range between the target geometry center and the radar is R 0 (t).The rotation angle of the target relative to the radar is θ(t).Then, the range R P (t) between any point P(x, y) in the target and the radar is: Considering the inertia of the conventional target motion and the short imaging time of ISAL, which is less than 1 s, maneuvering target motion can only approximate to the second-order motion component.In other words, R 0 (t) and θ(t) can be approximated as the second-order function of t 2 : (2) where R 0 is the initial range, v 0 is the initial radial velocity, a is the radial acceleration, θ 0 is the initial rotation angle, ω is the rotation angular velocity and Ω is the rotation angular acceleration.
Electronics 2018, 7, x FOR PEER REVIEW 3 of 20 relative motion between radar and target is considered.Where O is the reference point and XOY is the rectangular coordinate system fixed on the target, the Y axis is the direction of the radar LOS.The target moves along the Y axis with the speed of r v and rotates around the O point at the angular velocity of e  .Suppose that at time t ( t is full-time and satisfies the equation  km t t t , where k t is the range fast time and m t is the azimuth slow time, 1, 2,..., mM  ), the range between the target geometry center and the radar is 0 () Rt.The rotation angle of the target relative to the radar is () t  .
Then, the range ()

P
Rt between any point ( , ) P x y in the target and the radar is: Considering the inertia of the conventional target motion and the short imaging time of ISAL, which is less than 1 s, maneuvering target motion can only approximate to the second-order motion component.In other words, 0 () Rt and () t  can be approximated as the second-order function of 2 t : where 0 R is the initial range, 0 v is the initial radial velocity, a is the radial acceleration, 0  is the initial rotation angle,  is the rotation angular velocity and Ω is the rotation angular acceleration.Because wavelength of ISAL is on the micron scale, it is sensitive to the motion of the target, so the effect of the fast time k t on the radial range cannot be ignored [19].When the pulse duration is short, the influence of k t on rotation components can be neglected.So, Equations ( 2) and (3) can be resolved as: Because wavelength of ISAL is on the micron scale, it is sensitive to the motion of the target, so the effect of the fast time t k on the radial range cannot be ignored [19].When the pulse duration is short, the influence of t k on rotation components can be neglected.So, Equations ( 2) and (3) can be resolved as: where R 0 (t m ) and v(t m ) represent the radial range and velocity varies with the azimuth time t m , respectively.According to the above equation, Equation (1) can be rewritten as: where R d (t m ) = x sin θ(t m ) + y cos θ(t m ).
ISAL usually uses ultra-bandwidth Linear Frequency Modulated (LFM) signals to achieve high range resolution, the expression of the transmitted signal can be written as: where T p is the width of the pulse, f c is the carrier frequency and k is the chirp-rate.Considering c v r + at, we ignore the effect of the target speed of the irradiating and receiving radar signals.Suppose the radar receives the echo signal of point P after time delay τ = 2R p (t)/c, the radar receiving signal is: To reduce the data rate, ISAL often uses optical heterodyne coherent detection [20] to handle the echo signals.Suppose the reference delay of the coherent pulses is τ re f = 2R re f (t m )/c, the reference signal is: Therefore, the output signal after optical heterodyne detection is: We substitute Equations ( 4) and (6) into Equation (11) to obtain a polynomial function about time t k : where , which is only related to the scattering point on the target in azimuth time.In practice, considering the impact of the target velocity on τ and τ re f , the envelope in the Equation (12) will cause a contraction in time.But the impact does not affect the analysis of the range dispersion, which is ignored here [1].P 0 is the phase term related only to the azimuth time t m , in which 2 f c ∆R mP /c is a necessary term for azimuth compression and has no influence on range compression.P 1 is the first-order phase term.The first item in P 1 is mainly affected by the high carrier signal frequency, which produces the signal in the pulse Doppler.For the uniformly accelerated moving target, the Doppler coupling time shift with azimuth time change is produced, resulting in a range move.The second item contains the range information 2k∆R mP /c, which is key to attaining range compression.P 2 is the chirp-rate phase term, mainly influenced by the ultra-high carrier frequency and large bandwidth.It is the root cause of the division and broadening of the peaks of the range.The range dispersion effect occurs if the conventional DFT is used for the compression processing of the range direction.From the expression P 2 , the chirp-rate term of all scattering points in the single pulse echo is the same fixed value, whereas for different pulse echoes, the chirp-rate rate varies with slow time t m .Therefore, processing the pulse echo sequence one at a time is necessary.P 3 and P 4 are the high-order phase terms.Because in a pulse period c av(t m ), T p 2 c/2k r T p and T p 3 a 2 c 2 /2k r T p , the influence of the P 3 and P 4 on the intra-pulse Doppler spectrum broadening can be ignored.
According to the above analysis, in the imaging time, the second-order polynomial approximation can appropriately reflect the motion state of the maneuvering target and meet the imaging needs.The effect of the third-and fourth-order terms can be ignored.The range echo signals after heterodyne detection can be approximated to multicomponent LFM signals with the same frequency modulation slope:

Range Imaging Based on ICPF-FRFT
FRFT is a kind of generalized Fourier transform that better focuses LFM signals [21].The FRFT of the signal s(t) is defined as: where p is the order of the FRFT, which can be any real number and α is the rotation angle; α = pπ/2.When α = π/2, FRFT becomes a traditional Fourier transform.K α (t, u) is the transformation operator, the expression is: FRFT is equivalent to projecting the signal on the frequency axis after the counterclockwise rotation α of the signal in the time-frequency plane.When the u axis of the FRFT is rotated to the time-frequency ridge of the signal, the amplitude of the signal projection to the fractional frequency u axis is maximized and the rotation angle at this time is called the best angle α k of rotation.Therefore, the projection of FRFT at the best angle of rotation can be used for range imaging and the imaging principle is shown in Figure 2.
FRFT is equivalent to projecting the signal on the frequency axis after the counterclockwise rotation  of the signal in the time-frequency plane.When the u axis of the FRFT is rotated to the time-frequency ridge of the signal, the amplitude of the signal projection to the fractional frequency u axis is maximized and the rotation angle at this time is called the best angle k  of rotation.
Therefore, the projection of FRFT at the best angle of rotation can be used for range imaging and the imaging principle is shown in Figure 2. FRFT requires the peak search method in the two-dimensional (2D) plane ( , ) u  to obtain the optimal rotation angle.Therefore, the effect of range image compression depends on the value of k  and its precision is easily influenced by the resolution of the search angle.The computation requirement is considerable in the high precision search.So, accuracy and computation are difficult to achieve.Paper and colleagues [22][23][24][25] proposed that ICPF can quickly estimate the chirp rate of LFM signals.The method only requires a 1D search and has good anti-noise performance and high estimation accuracy without being affected by subjective factors such as search resolution.Therefore, in this paper, ICPF was firstly used to estimate the modulation frequency of the optical heterodyne output signal and then the optimal rotation angle and the order of the FRFT were calculated.Finally, the range compression imaging was completed by FRFT at the optimal rotation angle.The ICPF definition of signal () xt is as follows: From the definition, ICPF is a kind of transformation that detects the chirp-rate of the signal, which can concentrate the signal energy on the chirp-rate of the signal, in line with the energy distribution of the linear frequency modulation signal.Since ICPF needs to calculate the 2  of the signal, using FFT for fast calculations is not possible.Therefore, we used the non-uniform fast Fourier transform (NUFFT) [26,27] to overcome the rigorous data sampling requirements of the FFT and to improve the algorithm's calculation speed.The NUFFT is defined as: FRFT requires the peak search method in the two-dimensional (2D) plane (α, u) to obtain the optimal rotation angle.Therefore, the effect of range image compression depends on the value of α k and its precision is easily influenced by the resolution of the search angle.The computation requirement is considerable in the high precision search.So, accuracy and computation are difficult to achieve.Paper and colleagues [22][23][24][25] proposed that ICPF can quickly estimate the chirp rate of LFM signals.
The method only requires a 1D search and has good anti-noise performance and high estimation accuracy without being affected by subjective factors such as search resolution.Therefore, in this paper, ICPF was firstly used to estimate the modulation frequency of the optical heterodyne output signal and then the optimal rotation angle and the order of the FRFT were calculated.Finally, the range compression imaging was completed by FRFT at the optimal rotation angle.The ICPF definition of signal x(t) is as follows: From the definition, ICPF is a kind of transformation that detects the chirp-rate of the signal, which can concentrate the signal energy on the chirp-rate of the signal, in line with the energy distribution of the linear frequency modulation signal.Since ICPF needs to calculate the τ 2 of the signal, using FFT for fast calculations is not possible.Therefore, we used the non-uniform fast Fourier transform (NUFFT) [26,27] to overcome the rigorous data sampling requirements of the FFT and to improve the algorithm's calculation speed.The NUFFT is defined as: where z l is non-uniform sampling time and x l is the corresponding non-uniform sampling position.
Here, interpolation time domain non-uniform sampling data z l is replaced by an interpolation index term to achieve fast non-uniform Fourier transform.
, where, K is the length of interpolation kernel function.According to P. O'Shea [27], the exponential function can be expanded as shown in Equation ( 23): where |ξ| ≤ π/c, c is the oversampling factor.Suppose . Substituting Equation (23) into Equation ( 22) yields a uniform frequency output: where The specific NUFFT implementation process is shown in Figure 3. First, the intermediate parameters µ l , ϕ k , φlm are calculated from the input non-uniform sampling data z l and the corresponding position x l .Then, the intermediate variable u j is calculated according to Equations (25).Finally, the corresponding frequency output value ẑk is calculated using Equation ( 24) with fast FFT.
is the oversampling factor.Suppose . Substituting Equation (23) into Equation ( 22) yields a uniform frequency output: where The specific NUFFT implementation process is shown in Figure 3.According to the above NUFFT principle, the non-uniform Fourier transform in the proposed ICPF is quickly implemented with NUFFT to reduce the computational complexity of the algorithm, so Equation ( 26) can be written as: Assuming N is the sampling point of a single pulse echo and M is the number of DFT search points, the complexity of the ICPF-based DFT direct calculation is O(N 2 M), whereas the complexity of the nonuniform Fourier transform calculation method is O(Nlog2(N)) [28].Assuming that K is the number of scanning points in the FRFT transform domain  , which is determined by the step size and range of  , the complexity of the discrete FRFT is O(NKlog2(N)).If we want an accurate k  , K is usually large.
The transformation and 2D search need to be coordinated [21,24], so the proposed NUFFT-based ICPF-FRFT algorithm does not need to perform any parameter search and has high anti-noise performance.These features enable the implementation of the ISAL imaging algorithm in real time.
As a result, the ICPF transform of the optical heterodyne output signal results in spikes only at its chirp-rate slope.The chirp-rate at the peak is k  : According to the above NUFFT principle, the non-uniform Fourier transform in the proposed ICPF is quickly implemented with NUFFT to reduce the computational complexity of the algorithm, so Equation ( 26) can be written as: where NUFFT τ 2 indicates the NUFFT operation on the variable τ 2 and FFT t indicates FFT operation on the variable t.
Assuming N is the sampling point of a single pulse echo and M is the number of DFT search points, the complexity of the ICPF-based DFT direct calculation is O(N 2 M), whereas the complexity of the non-uniform Fourier transform calculation method is O(Nlog 2 (N)) [28].Assuming that K is the number of scanning points in the FRFT transform domain α, which is determined by the step size and range of α, the complexity of the discrete FRFT is O(NKlog 2 (N)).If we want an accurate α k , K is usually large.The transformation and 2D search need to be coordinated [21,24], so the proposed NUFFT-based ICPF-FRFT algorithm does not need to perform any parameter search and has high anti-noise performance.These features enable the implementation of the ISAL imaging algorithm in real time.
As a result, the ICPF transform of the optical heterodyne output signal results in spikes only at its chirp-rate slope.The chirp-rate at the peak is µ k : Calculate the rotation angle α k and FRFT order p k corresponding to the tuning frequency, which are the best FRFT rotation angle and order, respectively.When using discrete FRFT calculations, the signal parameters need to be dimension normalized [29].The relationship between the chirp-rate and the rotation angle is provided in Equation ( 28) and the FRFT order is provided in Equation ( 29): The optical heterodyne signal in Equation ( 12) is subjected to p k order FRFT: ) where A(u) = 1 + j2P 2 exp(−j2πP 2 u 2 ).From the result, the peak value of the signal is obtained at u = P 1 sin α k .That is, the range compression of the echo signal is achieved by one FRFT and the phase exp(j2πP 0 ) of the azimuth compression is retained.

Feature of Azimuth Echo Signal
In order to facilitate analysis, we assumed that the radar echo has completed range compression and motion compensation, so the echo signal can be converted into a turntable model with centroid as the reference point and the azimuth echo signal at the point can be expressed as: where σ is the amplitude of the signal after the motion compensation.When the target is maneuvering, θ(t) can be expanded into a function of time t according to Taylor [30] due to the inertia of space targets.For a space target with certain inertia, the ISAL imaging time is shorter and the cumulative rotation angle required by the imaging is smaller, so the motion of the target and radar can approximate to the second-order component, meaning it approximates the uniform acceleration motion.
where θ 0 is the initial rotation angle, ω is the rotation angular velocity and Ω is the rotation angular acceleration.
As the ISAL wavelength is in the order of µm, to achieve the imaging resolution of mm magnitudes, the required rotation accumulation angle is in the order of mrad, so the following small angle approximation conditions are satisfied: According to Equations ( 29) and ( 31), the P point azimuth echo can be approximated to a linear frequency modulated signal. where Electronics 2018, 7, 148 9 of 20 In practice, multiple scattering points with different intensities are distributed in the same range cell, so the azimuth echo becomes a multicomponent LFM signal with a different LFM rate: where K is the number of scattered points and φ 0i , k ai and f ai satisfy Equations ( 35)-(37), respectively.

Azimuth Compression Based on ICPF-FRFT
ICPF can effectively suppress the cross and pseudo peaks caused by the interference of multicomponent signals, so the ICPF-FRFT can be used for imaging azimuth signals but the strong signal components affect the detection of the weak signal components.Therefore, we combined the CLEAN technique with ICPF-FRFT to estimate the strong to weak signals.The frequency modulation slope of the signal was calculated and then FRFT was used to image the signal components of different frequency modulation slopes.The imaging procedure for the azimuth compression based on ICPF-FRFT is shown in Figure 4.
In practice, multiple scattering points with different intensities are distributed in the same range cell, so the azimuth echo becomes a multicomponent LFM signal with a different LFM rate: where K is the number of scattered points and 0i  , ai k and ai f satisfy Equations ( 35)-(37), respectively.

Azimuth Compression Based on ICPF-FRFT
ICPF can effectively suppress the cross and pseudo peaks caused by the interference of multicomponent signals, so the ICPF-FRFT can be used for imaging azimuth signals but the strong signal components affect the detection of the weak signal components.Therefore, we combined the CLEAN technique with ICPF-FRFT to estimate the strong to weak signals.The frequency modulation slope of the signal was calculated and then FRFT was used to image the signal components of different frequency modulation slopes.The imaging procedure for the azimuth compression based on ICPF-FRFT is shown in Figure 4.

Calculate the optimum rotation angle
Range  The concrete steps are as follows: Step 1: Calculate the ICPF of a range cell of the echo signal and estimate the frequency modulation rate of the strongest signal component: where ŝ(φ 0i , f ai , k ai , t m ) is the ith signal component.
Step 2: Calculate the best rotation angle α i and the corresponding order p ki , according to Equations ( 28) and ( 29), respectively.Then, calculate FRFT of the range cell signal in the range [α i − ∆, α i + ∆] where ∆ is the calculation error of α i .Search the peak to obtain the corresponding position u i : Ŝα,l (u) (41) Step 3: Separate the peak point by using the CLEAN technique to construct a narrowband filter W(u i ) centered on u i .Filter the strongest component and the peak value Ŝi αi (u) is considered the azimuth focusing image of the ith component.
Step 4: Transform the rest of the signal to the time domain using FRFT with a rotation angle of −α i .
Step 5: Repeat the above steps until all the scattered points in the current range cell are separated.This separation can be judged by when the residual signal component energy E of the ith range cell is less than a certain energy threshold E H , which is usually 5% of the original signal [31,32].
Step 6: The target image is obtained by scaling the scattered images u = u/ sin α and stacking them linearly.
Step 7: The 2D ISAL images can be obtained by using the above methods according to the sequence numbers of the range cells.

Imaging Procedure
For a maneuvering target with approximately uniformly accelerated motion, the ISAL imaging algorithm flow is shown in Figure 5.

Experimental Results of Simulation and Real Data
In order to verify the effectiveness of the proposed algorithm, simulation and real data experiments were completed.Some other imaging algorithms were considered for comparison.

Experimental Simulation Results
The ladar and target simulation parameters used in the simulation experiment are provided in Table 1, which refers to the scattering point model in Papers [1,11].The simulation model shown in

Experimental Results of Simulation and Real Data
In order to verify the effectiveness of the proposed algorithm, simulation and real data experiments were completed.Some other imaging algorithms were considered for comparison.

Experimental Simulation Results
The ladar and target simulation parameters used in the simulation experiment are provided in Table 1, which refers to the scattering point model in Papers [1,11].The simulation model shown in Figure 6 is a plane model that contains 52 scattering points.The ladar parameters are typical parameters that can be realized and the range resolution was 0.001 m.We assumed that the positional relationship between the ladar and the target was as shown in Figure 1.The target motion parameters were set as the turntable motion parameters as shown in Table 1. Figure 7 is the smooth pseudo Wigner distribution (SPWVD) time frequency graph of the 128th pulse echo.The ISAL single echo signal is a multicomponent LFM signal with the same modulation frequency, which confirms the analysis of the echo signal in Section 2. Therefore, the compression of all scattering points can be accomplished through one compression of a range.Figure 8 shows the 128th pulse range compression result using the DFT method.From the display results, the direct use of DFT compression results in a serious dispersion effect of the range image and the scattering points of the adjacent resolution units form.Serious mutual interference occurs for scattering points of adjacent resolution units.Figure 9 shows the 128th pulse range compression result using the ICPF-FRFT method proposed in this paper.As seen in the figure, the results show that a better compression effect was achieved and the range dispersion was eliminated.Figure 10 provides the SPWVD time frequency graph of the 79th (the lower edge of the aircraft) and 128th (the aircraft range center) range cell.From the diagram, due to the short ISAL imaging time, the target azimuth echo signals can be approximated as a multi-component LFM signal, even if there is a third order rotational component (angle acceleration rate) in the target.This demonstrates a slope with different slopes on the time frequency graph, which is the same as the previous theoretical analysis.So, dividing the scattering imaging into scattering points on different range cells was necessary.Figure 11a shows the result of traditional DFT azimuthal compression, which highlights that the scattered points from the center of the azimuth were seriously defocused-a poor imaging result.For comparison, we also provide three instantaneous Doppler (RID) imaging results based on STFT, WVD and SPWVD, as shown in Figures 11b-d, respectively.The image results used the 24th frame, which was t = 0.116 s.From the results of STFT in Figure 11b, the time-frequency resolution is affected by the window function and the azimuth defocus was severe.From Figure 11c, WVD can improve upon the time-frequency resolution but since the azimuth echo is a multi-component signal, the imaging result produces a cross term, so the imaging results are poorly readable and cannot identify the target.The SPWVD provides windowing and smoothing of WVD, so it weakens the cross terms but the time-frequency resolution also decreases.From Figure 11d, we can verify that the SPWVD has no cross-scattering point compared to the result of Figure 11c but the resolution is reduced.From the four results in Figure 11, the direct imaging range-Doppler (RD) algorithm and the three rangeinstantaneous Doppler (RID) imaging methods cannot achieve better imaging results and the azimuth defocus still exists.Figure 10 provides the SPWVD time frequency graph of the 79th (the lower edge of the aircraft) and 128th (the aircraft range center) range cell.From the diagram, due to the short ISAL imaging time, the target azimuth echo signals can be approximated as a multi-component LFM signal, even if there is a third order rotational component (angle acceleration rate) in the target.This demonstrates a slope with different slopes on the time frequency graph, which is the same as the previous theoretical analysis.So, dividing the scattering imaging into scattering points on different range cells was necessary.Figure 10 provides the SPWVD time frequency graph of the 79th (the lower edge of the aircraft) and 128th (the aircraft range center) range cell.From the diagram, due to the short ISAL imaging time, the target azimuth echo signals can be approximated as a multi-component LFM signal, even if there is a third order rotational component (angle acceleration rate) in the target.This demonstrates a slope with different slopes on the time frequency graph, which is the same as the previous theoretical analysis.So, dividing the scattering imaging into scattering points on different range cells was necessary.Figure 11a shows the result of traditional DFT azimuthal compression, which highlights that the scattered points from the center of the azimuth were seriously defocused-a poor imaging result.For comparison, we also provide three instantaneous Doppler (RID) imaging results based on STFT, WVD and SPWVD, as shown in Figures 11b-d, respectively.The image results used the 24th frame, which was t = 0.116 s.From the results of STFT in Figure 11b, the time-frequency resolution is affected by the window function and the azimuth defocus was severe.From Figure 11c, WVD can improve upon the time-frequency resolution but since the azimuth echo is a multi-component signal, the imaging result produces a cross term, so the imaging results are poorly readable and cannot identify the target.The SPWVD provides windowing and smoothing of WVD, so it weakens the cross terms but the time-frequency resolution also decreases.From Figure 11d, we can verify that the SPWVD has no cross-scattering point compared to the result of Figure 11c but the resolution is reduced.From the four results in Figure 11, the direct imaging range-Doppler (RD) algorithm and the three rangeinstantaneous Doppler (RID) imaging methods cannot achieve better imaging results and the azimuth defocus still exists.Figure 11a shows the result of traditional DFT azimuthal compression, which highlights that the scattered points from the center of the azimuth were seriously defocused-a poor imaging result.For comparison, we also provide three instantaneous Doppler (RID) imaging results based on STFT, WVD and SPWVD, as shown in Figure 11b-d, respectively.The image results used the 24th frame, which was t = 0.116 s.From the results of STFT in Figure 11b, the time-frequency resolution is affected by the window function and the azimuth defocus was severe.From Figure 11c, WVD can improve upon the time-frequency resolution but since the azimuth echo is a multi-component signal, the imaging result produces a cross term, so the imaging results are poorly readable and cannot identify the target.The SPWVD provides windowing and smoothing of WVD, so it weakens the cross terms but the time-frequency resolution also decreases.From Figure 11d, we can verify that the SPWVD has no cross-scattering point compared to the result of Figure 11c but the resolution is reduced.From the four results in Figure 11, the direct imaging range-Doppler (RD) algorithm and the three range-instantaneous Doppler (RID) imaging methods cannot achieve better imaging results and the azimuth defocus still exists.cells from each other, proving the effectiveness of the square method and further illustrates that the ISAL can be applied in real situations.High precision (millimeter level) imaging of space targets is now being performed.
Figures 13a-c are schematic diagrams of the separation of the first three peaks of the 95 range cell on the lower wing of the aircraft.The left image is the search process of the FRFT peak in the ICPF estimation error, the right side is the Clean processing for the peak point, the red frame is a narrow band filter and the frame content is the transverse focus image.Through this process, the azimuth scattering points with different intensities were separated and appropriately imaged.Figure 12 is the result of azimuth compression using ICPF-FRFT.In Figure 12a, the CLEAN technique is not used, whereas it is applied in Figure 12b.Figure 12a shows that the imaging results of a single scattering point on a certain range cell are good.With multiple scattering points, the strong signal suppresses the separation of the weak signal due to the different signal intensities.Most of the signal is strong and is the signal of its side lobe component, and the weak signal is missing.Figure 12b shows that when the CLEAN technique was used to separate the scattering points on different range cells from strong to weak separation imaging, a better focusing effect was achieved.The diagram demonstrates that this method can effectively separate the scattering points of two range cells from each other, proving the effectiveness of the square method and further illustrates that the ISAL can be applied in real situations.High precision (millimeter level) imaging of space targets is now being performed.
Figure 13a-c are schematic diagrams of the separation of the first three peaks of the 95 range cell on the lower wing of the aircraft.The left image is the search process of the FRFT peak in the ICPF estimation error, the right side is the Clean processing for the peak point, the red frame is a narrow band filter and the frame content is the transverse focus image.Through this process, the azimuth scattering points with different intensities were separated and appropriately imaged.In order to quantitatively evaluate the effectiveness of the proposed ICPF-FRFT algorithm, image entropy, contrast and running time were used to illustrate the imaging quality of the algorithm.Suppose the acquired ISAL image is f (n, k), where n and k are the range and azimuth number of the sampling unit, respectively.The definition of image entropy is: where F is the total energy of the ISAL image.The image entropy is small when the image is well-focused.Conversely, a large image entropy indicates that the compensation effect is worse.The definition of image contrast is: where E(•) represents the average operation.The image contrast is large when the image is well-focused.Conversely, a small image entropy indicates that the compensation effect is worse.
The results of the proposed ICPF-FRFT algorithm compared with the RD algorithm, FRFT algorithm, the ICPF-FRFT algorithm without CLEAN technique and three RID imaging methods based on STFT, WVD and SPWVD (Table 2).From the table, the algorithm proposed in this paper has smaller image entropy and a larger image contrast than the other algorithms, which shows that the image quality of the algorithm proposed in this paper is better.Notably, although all the indexes of the CLEAN technique are better, the loss of the scattering points cannot correctly reflect the distribution of the target scattering point, so the imaging quality was not the best.When the FRFT imaging algorithm with a small step size is performed directly, the result can reach an image entropy and contrast close to that of the proposed ICPF-FRFT algorithm paper but considerable computation time is required, indicating that the proposed algorithm is more efficient.

Experimental Results of Real Data
Since no ISAL data have been published to date, the research on ISAL at this stage is mainly based on simulation data to verify algorithms.However, considering the problem of azimuth Doppler time-varying when imaging a maneuvering target in ISAR, ISAR is consistent with ISAL imaging in the pursuit of azimuth focusing.Therefore, the ISAL azimuth imaging algorithm based on ICPF-FRFT is also suitable for imaging ISAR maneuvering targets but the radar signal bandwidth in ISAR is much smaller than ISAL, so the spread over the range can be ignored.To further validate the effectiveness of the algorithm, the publicly available Boeing B727 ISAR aircraft data from Victor C. Chen of the U.S. Naval Research Laboratory (Washington, DC, USA) was used for experimental verification [33].The data included 256 continuous pulses with a carrier frequency of 9 GHz, a bandwidth of 150 MHz and a pulse repetition rate of 20 kHz.Range compression and the motion compensation for the data were completed.
The imaging results using the imaging method proposed in this paper and other comparison methods described in the previous section are shown in Figure 14.The evaluation indexes of each imaging result are shown in Table 3.As can be seen from Figure 14a, the azimuth defocusing that occurred when using the RD algorithm was severe.From Figure 14b-d, the RID imaging results are related to the time-frequency method used, in which the time-frequency resolution of STFT was the worst and WVD had the highest time-frequency resolution but the cross-term was the most serious and SPWVD was somewhere in between.The result of the notable time-frequency imaging method shows that as the azimuth Doppler dynamically changes, the results displayed at different azimuths are different.In addition, some weak scattering point energy loss occurs, as shown in the wing part of the figure.Some scattering points are missing.It can be seen from Figure 14e-g that all three imaging methods can effectively improve the azimuth focusing effect but the direct FRFT requires a long computation time to achieve the same focusing effect as ICPF-FRFT.However, for the No-CLEAN technique, although all the indicators are superior, this algorithm only focuses on the strong scattering point, resulting in a lack of partial scattering points.Considering the minimum entropy, contrast and running time, the ICPF-FRFT algorithm is optimal, which is consistent with the results of the previous simulation analysis.

where 2 NUFFT
 indicates the NUFFT operation on the variable 2  and t FFT indicates FFT operation on the variable t .

Figure 4 .
Figure 4.The imaging procedure for the azimuth compression.

Figure 11 .Figure 11 .
Figure 11.Azimuth compression via classical methods: (a) direct DFT, (b) RID image based on STFT, (c) RID image based on WVD and (d) RID image based on SPWVD.

Figure 12 .
Figure 12.Azimuth compression via ICPF-FRFT: (a) without the CLEAN technique and (b) with the CLEAN technique.

Figure 13 .Figure 13 .
Figure 13.Separation imaging of the 95th range cell: (a) peak one, (b) peak two and (c) peak three.

Table 2 .
Comparison of simulation aircraft imaging results.