Micro-Doppler Effect Removal in ISAR Imaging by Promoting Joint Sparsity in Time-Frequency Domain

For micromotion scatterers with small rotating radii, the micro-Doppler (m-D) effect interferes with cross-range compression in inverse synthetic aperture radar (ISAR) imaging and leads to a blurred main body image. In this paper, a novel method is proposed to remove the m-D effect by promoting the joint sparsity in the time-frequency domain. Firstly, to obtain the time-frequency representations of the limited measurements, the short-time Fourier transform (STFT) was modelled by an underdetermined equation. Then, a new objective function was used to measure the joint sparsity of the STFT entries so that the joint sparse recovery problem could be formulated as a constrained minimization problem. Similar to the smoothed l0 (SL0) algorithm, a steepest descend approach was used to minimize the new objective function, where the projection step was tailored to make it suitable for m-D effect removal. Finally, we utilized the recovered STFT entries to obtain the main body echoes, based on which cross-range compression could be realized without m-D interference. After all contaminated range cells were processed by the proposed method, a clear main body image could be achieved. Experiments using both the point-scattering model and electromagnetic (EM) computation validated the performance of the proposed method.


Introduction
In inverse synthetic aperture radar (ISAR) imaging, mechanical rotations or vibrations of structures on a target may introduce additional frequency modulations on the returned signal, known as the micro-Doppler (m-D) effect [1][2][3]. The main body image is usually blurred because of the interference of the rotating or vibrating scatterers, which are also called micromotion scatterers [3]. Therefore, the m-D effect should be removed in order to obtain a clear image of the main body.
Micromotion targets with rotating structures can often be found in real world scenarios, e.g., a ship with scanning antennas or an aircraft in flight equipped with turbofans. A micromotion scatterer with a large rotating radius will generate a m-D signal which can be depicted in the form of sinusoidal modulations in the spectrogram after range compression [4,5]. In contrast, the main body signal in the spectrogram has the shape of straight lines [5]. Based on the different shapes in the spectrogram, many methods [4][5][6][7] have been proposed to eliminate the m-D signal. Li and Ling [4] put forward an adaptive chirplet decomposition algorithm in which the returned signal is decomposed into a family of chirplet functions. The main body signal and the micromotion signal subsequently are separated according to their distinct chirp rates. However, this algorithm has high computation cost because of the large chirplet dictionary. To remove the sinusoidal m-D interference, the method in Q. Zhang et al. [6] extracts the straight lines in the spectrogram using the Hough transform. Similarly, the methods in H.C. Liu et al. and L. Sun et al. [5,7] only recover the straight lines via the methods in M. Bevacqua et al., Y. Liu et al., and S. Sun et al. [18,19,21], we propose a new objective function in this paper which can achieve tight approximation to the original l 2 /l 0 -norm with proper parameters. The new objective function was minimized, similar to the smoothed l 0 (SL0) algorithm approach [23] which consisted of two steps, i.e., the steepest descend step and the projection step. Since the m-D signal was suppressed in the steepest descend step, in this paper, the projection step of SL0 was tailored to make it suitable for m-D effect removal. By minimizing the new objective function, the joint sparsity in the time-frequency domain was promoted and the m-D interference could be effectively removed. Subsequently, we utilized the recovered STFT entries to obtain the main body echoes free of m-D effect. Finally, the cross-range compression [24] of the main body echoes was realized through the sparse Bayesian learning (SBL) algorithm [25]. When the m-D effect in all of the contaminated range cells was removed by the proposed method, a clear main body image could be achieved.
The rest of this paper is organized as follows. Section 2 presents the signal model of ISAR imaging for micromotion targets with rotating parts. Section 3 demonstrates the joint sparsity of the main body signal in the time-frequency domain. The proposed ISAR imaging method is elaborated in Section 4, where we introduce the joint sparse recovery problem based on the l 2 /l 0 -norm and briefly review the SL0 algorithm. In Section 5, we discuss experiments based on the point-scattering model and electromagnetic (EM) computation which were conducted to validate the effectiveness of the proposed method. Finally, conclusions are drawn in Section 6.
We introduce the following notations in this paper. Bold uppercase letters and bold lowercase letters are reserved for matrices and vectors, respectively. A ∈ C M×N denotes a matrix of size M × N with complex elements. A F denotes the Frobenius norm of A. For a vector x, its ith element is denoted by x i . x 0 , x 1 , and x 2 represent the l 0 -norm, l 1 -norm and l 2 -norm of x, respectively. (·) T and (·) H stand for the transpose and the conjugate transpose of a matrix or a vector.

ISAR Imaging Model
To simplify the analysis, we consider the point-scattering model [4] to illustrate the m-D effect. As shown in Figure 1, the radar is located at the origin O of the coordinate system XOY. The target center u(x u , y u ) is located at the Y-axis which indicates the line of sight (LOS). Without loss of generality, we assume that the target moves within the 2-dimensional (2-D) imaging plane XOY with velocity v and the motion compensation [26] has been accomplished. The projection of v along the direction of the X-axis is denoted by v x which generates the aspect angle variation utilized in ISAR imaging. p x p , y p and q x q , y q denote the main body scatterer and the micromotion scatterer, respectively. O (x o , y o ) is the rotating center of the micromotion part. q rotates around O with radius r q , angle frequency ω q , and initial rotation angle θ q . At the initial processing time, the distances from p, q, and u to the radar are R p , R q , and R u , respectively.
The transmitted chirp signal is wheret and t m denote the fast time and the slow time, respectively. T is the pulse duration, f c is the carrier frequency, and γ is the chirp rate. rect(·) represents the rectangular function, which is defined as To reduce the received effective bandwidth, the dechirp method is applied. Taking the target center as the reference point in the dechirping process, the reference signal is given by where T re f is the duration of the reference signal with T re f > T, and c denotes the speed of light. After dechirping and removing the residual video phase, the echo can be expressed as where σ p and σ q denote the backscattering coefficients of the pth main body scatterer and the qth micromotion scatterer, respectively. After the Fourier transform to Equation (4) along thet dimension, the range compression is achieved and the spectrogram can be written as wheref is the fast frequency and λ is the wavelength. A p = σ p T and A q = σ q T denote the complex coefficients of scatterers p and q, respectively. It is assumed that the target has a constant velocity v x in a short coherent processing interval (CPI). Thus, the aspect angle variation of the target is calculated to be v x t m /R u . Because the target is far from the radar, the instantaneous distances from p and q to the reference point can be approximately expressed as Substituting Equations (6) and (7) into Equation (5) yields It can be seen from Equation (8) that the peaks of the spectrogram are located at and the range resolution is ρ r = c/(2γT). For the targets with large rotors, the rotating radius r q might be much larger than the range resolution. Equation (8) indicates that the signal of the micromotion scatterer migrates through the range cells in the spectrogram and has the shape of a sinusiod. In contrast, the position of the main body signal is slow time invariant, which has the shape of a straight line. The sinusoid in the spectrogram will seriously affect the cross-range compression, i.e., the integration in the slow time domain, and lead to a smeared ISAR image. To remove the  [5,7] only recover the straight lines where the sinusoids are eliminated. However, when the target has small rotors, r q might be smaller than half of the range resolution. In this condition, the signal of the micromotion scatterer is located in a constant range cell and also has the shape of a straight line. The spectrogram can be approximately expressed as Because the main body signal and the signal of the micromotion scatterer have the same straight line shapes in the spectrogram, the methods in H.C. Liu et al., Q. Zhang et al., and L. Sun et al. [5][6][7] are invalid to remove the m-D effect. According to Equation (11), the Doppler frequencies of the main body scatterer and the micromotion scatterer are It is obvious from Equations (12) and (13) that the main body scatterer has a constant Doppler frequency while the Doppler frequency of the micromotion scatterer is slow time variant. Therefore, it is possible to separate the main body signal from the micromotion signal in the time-frequency domain.

Joint Sparsity of Main Body Signal in Time-Frequency Domain
Typical time-frequency transforms include the Wigner-Ville distribution (WVD) [27] and the STFT. Although WVD has better time-frequency resolution than STFT, it suffers from cross-terms for multiple signal components. Thus, we chose the STFT to achieve the time-frequency representations of the echoes in the contaminated range cell. In order to obtain the STFT entries, the echoes in the contaminated range cell were first segmented into narrow time intervals so that the signal in each segment could be considered stationary. Then, the Fourier transform was carried out for the windowed signal in each segment to obtain the local spectral representations. In this paper, we chose the rectangular window function and the segmentation of the echoes could be realized by a sliding window with proper time step.
To demonstrate the joint sparsity of the main body signal in the time-frequency domain, we conducted the following simulation. The target model was composed of nine main body scatterers and two micromotion scatterers, as described in Figure 2. The micromotion scatterers 1 q and 2 q Figure 1. Inverse synthetic aperture radar (ISAR) imaging geometry with rotating scatterers.

Joint Sparsity of Main Body Signal in Time-Frequency Domain
Typical time-frequency transforms include the Wigner-Ville distribution (WVD) [27] and the STFT. Although WVD has better time-frequency resolution than STFT, it suffers from cross-terms for multiple signal components. Thus, we chose the STFT to achieve the time-frequency representations of the echoes in the contaminated range cell. In order to obtain the STFT entries, the echoes in the contaminated range cell were first segmented into narrow time intervals so that the signal in each segment could be considered stationary. Then, the Fourier transform was carried out for the windowed signal in each segment to obtain the local spectral representations. In this paper, we chose the rectangular window function and the segmentation of the echoes could be realized by a sliding window with proper time step.
To demonstrate the joint sparsity of the main body signal in the time-frequency domain, we conducted the following simulation. The target model was composed of nine main body scatterers and two micromotion scatterers, as described in Figure 2. The micromotion scatterers q1 and q2 rotated around the origin counterclockwise with radii of 0.1 m and 0.2 m, and with rotating frequencies of 10 Hz and 3 Hz, respectively. The initial rotating phases were π/2 rad and 0 rad, respectively. The target moved along the cross-range direction with a velocity of 300 m/s. The radar carrier frequency was 10 GHz, the bandwidth was 600 MHz, the pulse repetition interval (PRI) was 1 ms, and the CPI was 1 s.

Joint Sparsity of Main Body Signal in Time-Frequency Domain
Typical time-frequency transforms include the Wigner-Ville distribution (WVD) [27] and the STFT. Although WVD has better time-frequency resolution than STFT, it suffers from cross-terms for multiple signal components. Thus, we chose the STFT to achieve the time-frequency representations of the echoes in the contaminated range cell. In order to obtain the STFT entries, the echoes in the contaminated range cell were first segmented into narrow time intervals so that the signal in each segment could be considered stationary. Then, the Fourier transform was carried out for the windowed signal in each segment to obtain the local spectral representations. In this paper, we chose the rectangular window function and the segmentation of the echoes could be realized by a sliding window with proper time step.
To demonstrate the joint sparsity of the main body signal in the time-frequency domain, we conducted the following simulation. The target model was composed of nine main body scatterers and two micromotion scatterers, as described in Figure 2. The micromotion scatterers 1 q and 2 q rotated around the origin counterclockwise with radii of 0.1 m and 0.2 m, and with rotating frequencies of 10 Hz and 3 Hz, respectively. The initial rotating phases were 2 π rad and 0 rad, respectively. The target moved along the cross-range direction with a velocity of 300 m/s. The radar carrier frequency was 10 GHz, the bandwidth was 600 MHz, the pulse repetition interval (PRI) was 1 ms, and the CPI was 1 s.  According to the radar system parameters, the range resolution was calculated to be 0.25 m, which was larger than twice the rotating radii of both q1 and q2. Thus, the micromotion signal did not migrate through the range cells in the spectrogram. As shown in Figure 3a, the micromotion signal was located in the 40th range cell and had the same straight line shape as the main body signal. In this situation, the methods in H.C. Liu et al., Q. Zhang et al., and L. Sun et al. [5][6][7] were invalid to eliminate the micromotion components. As a result from the m-D effect, the application of a fast Fourier transform (FFT) to the spectrogram along the slow time dimension would lead to a blurred main body image, as depicted in Figure 3b. We selected 50 cross-range cells where the target was located in the imaging result to display. According to the radar system parameters, the range resolution was calculated to be 0.25 m, which was larger than twice the rotating radii of both 1 q and 2 q . Thus, the micromotion signal did not migrate through the range cells in the spectrogram. As shown in Figure 3a, the micromotion signal was located in the 40th range cell and had the same straight line shape as the main body signal. In this situation, the methods in H.C. Liu et al., Q. Zhang et al., and L. Sun et al. [5][6][7] were invalid to eliminate the micromotion components. As a result from the m-D effect, the application of a fast Fourier transform (FFT) to the spectrogram along the slow time dimension would lead to a blurred main body image, as depicted in Figure 3b. We selected 50 cross-range cells where the target was located in the imaging result to display. To obtain the time-frequency representations, we applied the STFT to the echoes in the 40th range cell of the spectrogram. The time step was 10 PRI and the window width was 30 PRI. Thus, there were 30 frequency bins in the time-frequency domain. As shown in Figure 4, the main body signal was located at the 15th and 16th frequency bins. Consequently, the support of the main body signal in the time-frequency domain was slow time invariant which indicated the joint sparsity of the main body signal. In contrast, the micromotion scatterers' Doppler frequencies were slow time variant and the micromotion signal had the shape of sinusoids, as illustrated in Figure 4. Because of the distinct patterns, the micromotion signal could be removed by promoting the joint sparsity in the  Figure 4, the main body signal was located at the 15th and 16th frequency bins. Consequently, the support of the main body signal in the time-frequency domain was slow time invariant which indicated the joint sparsity of the main body signal. In contrast, the micromotion scatterers' Doppler frequencies were slow time variant and the micromotion signal had the shape of sinusoids, as illustrated in Figure 4. Because of the distinct patterns, the micromotion signal could be removed by promoting the joint sparsity in the time-frequency domain, which is detailed in the next section. To obtain the time-frequency representations, we applied the STFT to the echoes in the 40th range cell of the spectrogram. The time step was 10 PRI and the window width was 30 PRI. Thus, there were 30 frequency bins in the time-frequency domain. As shown in Figure 4, the main body signal was located at the 15th and 16th frequency bins. Consequently, the support of the main body signal in the time-frequency domain was slow time invariant which indicated the joint sparsity of the main body signal. In contrast, the micromotion scatterers' Doppler frequencies were slow time variant and the micromotion signal had the shape of sinusoids, as illustrated in Figure 4. Because of the distinct patterns, the micromotion signal could be removed by promoting the joint sparsity in the time-frequency domain, which is detailed in the next section.

Problem Formulation for Joint Sparse Recovery
The echoes in the range cell contaminated by the m-D effect are denoted by s, and we assume that the number of the received pulses is V . The window width and the time step of the STFT are

Problem Formulation for Joint Sparse Recovery
The echoes in the range cell contaminated by the m-D effect are denoted by s, and we assume that the number of the received pulses is V. The window width and the time step of the STFT are denoted by N and G, respectively. Thus, the number of the segments in the STFT is M = (V − N)/G + 1, where · denotes rounding to the nearest integer towards infinity.
The discrete Fourier matrix for each segment is where W = exp(j2π/N). Denoting the STFT entries in Figure 4 by X, the corresponding STFT entry vector is x, which is generated by stacking the columns of X into a single vector. Then, the STFT of s can be written in a compact matrix form where H ∈ C K×V denotes the STFT matrix and K = MN is the number of the STFT entries in x. 0 N×G denotes a matrix of size N × G, whose elements are equal to 0. The echoes in the contaminated range cell can be modeled as follows: where φ ∈ C V×1 denotes the noise, and D = H H H −1 H H is the Moore-Penrose inverse of the STFT matrix H. However, the data samples of the echoes might be randomly missing in real world scenarios where strong electromagnetic interference or sensor failure prevents effective observation. The data samples under strong electromagnetic interference must be discarded and the remaining L samples constitute the measurement vector y ∈ C L×1 , which can be expressed as where n ∈ C L×1 is the noise. A ∈ C L×K is formed by omitting the rows of D which correspond to the missing samples. Because the remaining samples are fewer than the STFT entries to recover, i.e., L < K, Equation (18) is underdetermined.
To facilitate the notation, we rearrange the STFT entries in the ith frequency bin into a vector According to the analysis in Section 3, the micromotion components in the time-frequency domain can be eliminated by promoting the joint sparsity of x. The joint sparse recovery problem is usually formulated by the following l 2 /l 0 minimization problem where ε bounds the l 2 -norm of the noise in the measurements. x 2,0 stands for the l 2 /l 0 -norm of x defined as [16] x 2,0 = β 0 where β = [ b 1 2 , b 2 2 , . . . , b N 2 ] T . It is obvious that x 2,0 is equivalent to the number of frequency bins which have nonzero STFT entries. Therefore, by solving the problem in Equation (20), we can find the solutionx which occupies the fewest frequency bins and the joint sparsity in the time-frequency domain can be promoted. However, the problem in Equation (20) is essentially a combinatorial optimization problem [17], which is difficult to solve. Inspired by the smoothed l 0 (SL0) algorithm, we propose a novel objective function to replace x 2,0 , so that the new problem is tractable. In the following subsections, the original SL0 algorithm is briefly reviewed at first. Then, we elaborate the proposed ISAR imaging method.

Brief Review of SL0 Algorithm
The basic idea of the SL0 algorithm is to approximate the l 0 -norm of x with the following function [23] where and It can be seen from Equations (22) and (24) that x 0 ≈ MN − F σ (x) for small values of σ. Furthermore, this approximation becomes equality as σ → 0 . The SL0 algorithm obtains the sparse solution of Equation (18) by solving the following problem To avoid trapping into local maxima, a decreasing sequence of σ is utilized, i.e., σ 1 , . . . , σ η , where σ i > σ j for i > j. η is the iteration number and the current maximizer of F σ (x) is used as the starting point of the next iteration. To maximize F σ (x) for a fixed σ, there are two steps in each iteration: the unconstrained maximization step x ← x + µ∇F σ (x) and the projection to the feasible set where A † = A H AA H −1 , and ∇F σ (x) is the gradient of F σ (x).

Proposed Method
Although it has been shown in H. Mohimani et al. [23] that the SL0 algorithm is effective to obtain the sparse solution of Equation (18), the joint sparsity of x cannot be promoted. It can be seen from Equation (22) that the SL0 algorithm only promotes the sparsity of the individual entries in x, and (MN − F σ (x)) is merely able to approximate x 0 rather than x 2,0 , which measures the joint sparsity of x. As a result, the SL0 algorithm is inappropriate to remove the m-D effect in ISAR imaging for micromotion targets. In this paper, we propose a smoothed l 2 /l 0 (SL2L0) algorithm, where x 2,0 in Equation (20) is approximated by a new objective function Since x 2,0 is equivalent to the number of frequency bins which have nonzero STFT entries, it is obvious from Equations (19), (24) and (27) that U σ (x) = x 2,0 when σ → 0 . Therefore, the joint sparsity can be promoted by minimizing the new objective function U σ (x) with small values of σ. The joint sparse recovery problem can be formulated as Similar to the SL0 algorithm, a decreasing sequence of σ is used to avoid trapping into local minimum and a steepest descend approach is employed to minimize U σ (x) in each iteration, i.e., where µ is the parameter controlling the step-size 2µσ 2 . The gradient ∇U σ (x) in Equation (29) is calculated to be ∇U σ (x) = 1 N ⊗ α x (30) where 1 N is a vector of length N with its entries equal to 1, ⊗ is the Kronecker product, is the Hadamard product, and α ∈ C N×1 is Since the micromotion components are eliminated by promoting the joint sparsity of x, we modify the projection step in Equation (26) as follows: where λ bounds the l 2 -norm of the micromotion components in the echoes. Equation (32) can be rewritten into an unconstrained optimization problem where z = x − x and y = y − Ax . τ is the regularization parameter related to λ. g(z) = z 1 is called the regularizer which indicates that the optimization problem in Equation (33) is convex. We use the fast iterative shrinkage-thresholding algorithm (FISTA) [28,29] to solve the problem in Equation (33) efficiently. The key concept in FISTA is known as Moreau's proximal operator, or proximal operator for short, which is expressed as and θ in Equation (34) is a positive constant with θ ≥ 1/ A 2 F . In summary, the successive steps of the proposed SL2L0 algorithm are given in Algorithm 1.

Iteration:
For i = 1, . . . , η By promoting the joint sparsity in the time-frequency domain, the output x of the SL2L0 algorithm only contains the main body components. Based on the recovered STFT entries, we can obtain the main body echoes free of m-D effect as y = Ax. Assuming that the residual noise in y is e ∈ C L×1 , y can be rewritten as where p ∈ C V×1 is the cross-range compression result of y, and Ω ∈ C L×V is a partial Fourier matrix with its (l, v)th element given by exp(−j2π f v t l ). f v is the vth Doppler frequency bin, and t l is the lth sampling point of the slow time. Subsequently, we used the sparse Bayesian learning (SBL) [25] algorithm to achieve the cross-range compression of y by recovering p in Equation (35). Although the discussions in this section are focused on the echoes in a single range cell, the proposed method is also suitable for other contaminated range cells. After all of the contaminated range cells are processed by the proposed method, a clear main body image can be achieved.

Experimental Results and Performance Comparisons
In this section, we discuss the experiments conducted using the echoes generated by the point-scattering model and EM computation, respectively. To validate the effectiveness, the imaging results of the proposed method and other existing methods are compared. The parameters of the proposed method used in the experiments are set as follows: . . , η, ez = 10 −4 , and ex = 10 −4 . Since ε is related to the noise level, this parameter is tuned manually in each experiment.

Experiments Using the Point-Scattering Model
The target model is depicted in Figure 2, and the radar system parameters are the same as in Section 3. We obtained the limited measurements by choosing 500 pulses randomly in the slow time domain, and the signal-to-noise ratio (SNR) was 20 dB. The number of cross-range cells to be recovered was 1000. The window width of the STFT was 30 PRI and the time step was 10 PRI. According to the above parameters, the size of A was calculated to be 500 × 2940. Similar to Figure 3b, we selected 50 cross-range cells to display in the imaging results.
To eliminate the micromotion components of the echoes in the contaminated range cell, both the method in R. Zhang et al. [9] and the L-statistics-based method utilize the time-frequency representations obtained by the STCS [11]. However, part of the main body signal in the time-frequency domain is missing due to the discarded pulses, as illustrated in Figure 5. Additionally, there is a strong m-D signal which has the shape of sinusoids. As a result, the histogram analysis in R. Zhang et al. [9] could not be used to remove the m-D signal, which might have been mistaken for the main body signal. Figure 6a shows the imaging result produced by the method in R. Zhang et al. [9]. Many spurious points exist in the 40th range cell due to the m-D interference and the main body scatterers are defocused. Because the measurements within the sliding window are much fewer than the frequency bins in STCS, the frequency resolution in the time-frequency domain is lower than the reciprocal of the CPI. Therefore, the L-statistics-based method cannot obtain the accurate support of the main body scatterers in the cross-range domain. As depicted in Figure 6b, some spurious points exist near the main body scatterers. In addition, one of the main body scatterers in the 40th range cell cannot be recovered because part of the main body signal is missing in Figure 5. In contrast to STCS, the proposed SL2L0 algorithm can remove the m-D signal by promoting the joint sparsity in the time-frequency domain, where only the main body signal with constant Doppler frequencies are preserved, as shown in Figure 7a. Based on the recovered STFT entries, the proposed method achieves a clear main body image in Figure 7b.
In addition to the visual results, we introduce the normalized mean square error (NMSE) and the image contrast (IC) as two metrics to compare the performance of different methods. A lower NMSE means the image has less recovery error, and a higher IC indicates the image is better focused. The NMSE is computed by [30] where U = 100 independent trials are conducted, andÎ k ∈ C R×V denotes the recovered result of the true main body image I in the kth trial. The IC is computed by X. Z. Gao et al. [24] whereÎ 2 k represents the intensity matrix of the recovered main body image with its (r, v)th element Ave(·) denotes the mean operator, which is defined as Figure 8 gives the NMSEs and ICs obtained by different methods against the undersampling factor [31], which is calculated as ψ = L/V, i.e., the number of measurements in slow time dimension divided by the number of cross-range cells to recover. For each ψ, the IC is averaged over 100 independent trials. It can be observed from Figure 8 that the proposed method achieves the lowest NMSEs and the highest ICs among all of the methods, which indicate the least recovery error and the best focused quality, respectively.

Experiments Using EM Computation
To further testify the performance of the proposed method, we used the echoes of a ship target generated by the graphic electromagnetic computing (GRECO) technique [32], which combines physical optical (PO) theory and physical theory of diffraction (PTD). The target model is described in Figure 9a, where the dashed circle indicates the antenna rotating with a frequency of 1 Hz. Figure 9b depicts the imaging geometry. The ship was 6.9 km from the airborne radar and located at the origin. The radar operated at 5 GHz, the bandwidth was 1 GHz, and the radar moved along the opposite direction of the X -axis with a constant velocity of 300 m/s. 256 pulses were collected in a CPI of 2.3 s and the SNR was 20 dB. We obtained the limited measurements by choosing 128 pulses randomly. Within each pulse, there were 128 complex range samples. The window width and the time step of the STFT were 20 PRI and 2 PRI, respectively. According to the aforementioned parameters, the size of A was calculated to be

Experiments Using EM Computation
To further testify the performance of the proposed method, we used the echoes of a ship target generated by the graphic electromagnetic computing (GRECO) technique [32], which combines physical optical (PO) theory and physical theory of diffraction (PTD). The target model is described in Figure 9a, where the dashed circle indicates the antenna rotating with a frequency of 1 Hz. Figure 9b depicts the imaging geometry. The ship was 6.9 km from the airborne radar and located at the origin. The radar operated at 5 GHz, the bandwidth was 1 GHz, and the radar moved along the opposite direction of the X -axis with a constant velocity of 300 m/s. 256 pulses were collected in a CPI of 2.3 s and the SNR was 20 dB. We obtained the limited measurements by choosing 128 pulses randomly. Within each pulse, there were 128 complex range samples. The window width and the time step of the STFT were 20 PRI and 2 PRI, respectively. According to the aforementioned parameters, the size of A was calculated to be

Experiments Using EM Computation
To further testify the performance of the proposed method, we used the echoes of a ship target generated by the graphic electromagnetic computing (GRECO) technique [32], which combines physical optical (PO) theory and physical theory of diffraction (PTD). The target model is described in Figure 9a, where the dashed circle indicates the antenna rotating with a frequency of 1 Hz. Figure 9b depicts the imaging geometry. The ship was 6.9 km from the airborne radar and located at the origin. The radar operated at 5 GHz, the bandwidth was 1 GHz, and the radar moved along the opposite direction of the X-axis with a constant velocity of 300 m/s. 256 pulses were collected in a CPI of 2.3 s and the SNR was 20 dB. We obtained the limited measurements by choosing 128 pulses randomly. Within each pulse, there were 128 complex range samples. The window width and the time step of the STFT were 20 PRI and 2 PRI, respectively. According to the aforementioned parameters, the size of A was calculated to be 128 × 2380.
The imaging results by different methods are shown in Figure 10. It can be seen from Figure 10a that the m-D effect led to a blurred main body image by the FFT-based method. Though the method in R. Zhang et al. [9] and the L-statistics-based method suppressed the m-D effect to some degree, there were residuals and some spurious points exist in Figure 10b,c. In contrast, the proposed method could remove the m-D effect through the SL2L0 algorithm and a clear image with the least spurious points is achieved in Figure 10d. the origin. The radar operated at 5 GHz, the bandwidth was 1 GHz, and the radar moved along the opposite direction of the X -axis with a constant velocity of 300 m/s. 256 pulses were collected in a CPI of 2.3 s and the SNR was 20 dB. We obtained the limited measurements by choosing 128 pulses randomly. Within each pulse, there were 128 complex range samples. The window width and the time step of the STFT were 20 PRI and 2 PRI, respectively. According to the aforementioned parameters, the size of A was calculated to be  The imaging results by different methods are shown in Figure 10. It can be seen from Figure 10a that the m-D effect led to a blurred main body image by the FFT-based method. Though the method in R. Zhang et al. [9] and the L-statistics-based method suppressed the m-D effect to some degree, there were residuals and some spurious points exist in Figure 10b,c. In contrast, the proposed method could remove the m-D effect through the SL2L0 algorithm and a clear image with the least spurious points is achieved in Figure 10d. To quantitatively evaluate the performance of different methods, we used the image entropy (IE). Generally, a smaller IE indicates a better focused quality of the image. The IE is defined by W. Qiu et al. [33] The IEs of the imaging results are given in Table 1. The proposed method achieved the minimum IE which verifies the superiority of the proposed method to other methods.
The computational time of different methods was recorded based on the MATLAB code and summarized in Table 2, where the results are given by the averages of 50 independent trials conducted on an Intel Core i7 machine at 3.50 GHz. The computational time of the FFT-based method was not recorded since it could not suppress the m-D effect, as shown in Figure 10a. It can be seen To quantitatively evaluate the performance of different methods, we used the image entropy (IE). Generally, a smaller IE indicates a better focused quality of the image. The IE is defined by W. Qiu et al. [33] where I is the power normalized image computed by The IEs of the imaging results are given in Table 1. The proposed method achieved the minimum IE which verifies the superiority of the proposed method to other methods.
The computational time of different methods was recorded based on the MATLAB code and summarized in Table 2, where the results are given by the averages of 50 independent trials conducted on an Intel Core i7 machine at 3.50 GHz. The computational time of the FFT-based method was not recorded since it could not suppress the m-D effect, as shown in Figure 10a. It can be seen from Table 2 that the proposed method is more computationally efficient than the method in R. Zhang et al. [9] and the L-statistics-based method. Table 1. Performance Evaluation by Image Entropy.

Method Image Entropy
The FFT-based method 5.977 The method in [9] 4.865 The L-statistics-based method 4.803 The proposed method 4.204 Table 2. Computational Time.

Method Computational Time (s)
The method in [9] 5.529 The L-statistics-based method 3.406 The proposed method 1.811

Conclusions
In ISAR imaging, the m-D effect generated by the micromotion scatterers with small rotating radii leads to a blurred main body image. It is claimed that the Doppler frequency of the micromotion scatterer is slow time variant while the main body scatterer has a constant Doppler frequency, indicating the joint sparsity property in the time-frequency domain. To remove the m-D effect, this paper proposes a novel algorithm named SL2L0, where a new objective function is used to measure the joint sparsity and is minimized in a way similar to the SL0 algorithm. We have tailored the projection step in the SL0 algorithm to make it suitable for m-D effect removal. Experimental results showed that the proposed SL2L0 algorithm can effectively remove the m-D effect by promoting the joint sparsity in the time-frequency domain and a clear main body image can be achieved after cross-range compression. Diverse metrics, including the NMSE, IC, and IE, were employed to validate the superiority of the proposed method to other existing methods.