ISAR Imaging of Ship Targets Based on an Integrated Cubic Phase Bilinear Autocorrelation Function

For inverse synthetic aperture radar (ISAR) imaging of a ship target moving with ocean waves, the image constructed with the standard range-Doppler (RD) technique is blurred and the range-instantaneous-Doppler (RID) technique has to be used to improve the image quality. In this paper, azimuth echoes in a range cell of the ship target are modeled as noisy multicomponent cubic phase signals (CPSs) after the motion compensation and a RID ISAR imaging algorithm is proposed based on the integrated cubic phase bilinear autocorrelation function (ICPBAF). The ICPBAF is bilinear and based on the two-dimensionally coherent energy accumulation. Compared to five other estimation algorithms, the ICPBAF can acquire higher cross term suppression and anti-noise performance with a reasonable computational cost. Through simulations and analyses with the synthetic model and real radar data, we verify the effectiveness of the ICPBAF and corresponding RID ISAR imaging algorithm.


Introduction
High-resolution inverse synthetic aperture radar (ISAR) imaging has attracted the attention of radar researchers in the past three decades due to its importance in civil and military applications [1][2][3][4]. Two challenges of the high-resolution ISAR imaging are the motion compensation and scatterer-dependent Doppler spread compensation [3,[5][6][7]. The motion compensation includes the translational range migration compensation, migration through resolution cells compensation and phase focusing [5][6][7]. The scatterer-dependent Doppler spread is corresponding to the coordinate of the scatterer, and the range-instantaneous-Doppler (RID) technique, whose essence is the parameters estimation, is developed to resolve it [3]. After three decades of research, there are mature processing methods for the motion compensation, such as the standard range alignment method for the translational range migration compensation, keystone transform for the migration through resolution cells compensation and phase gradient autofocus method for the phase focusing [5][6][7], while the scatterer-dependent Doppler spread compensation with the RID technique remains a great challenge [8][9][10]. This is because, in each range cell, multicomponent scatterers exist, and the RID technique (or the parameters estimation algorithm) needs overall consideration of the computational cost, cross term suppression, anti-noise performance, etc. [3,[10][11][12][13][14][15][16]. Most researchers currently focus on the scatterer-dependent Doppler spread compensation with the RID technique [8][9][10][11][12], which is also the focus of this paper. Figure 1 shows the geometry for ISAR imaging of the ship target. The X, Y and Z axes of the Cartesian coordinate system overlap with the longitudinal, horizontal and vertical axes of the ship target, respectively. R is the unit vector of the radar line of sight (parallel to the range cell).
roll Ω , pitch Ω and yaw Ω denote the angular rotational vectors of the ship target rotating around the X, Y and Z axes, respectively. Ω is the synthetic vector of roll Ω , pitch Ω and yaw Ω , and can be decomposed into e Ω and R Ω in the plane determined by R and Ω . R Ω does not cause the radial motion and thus has no effect on the phase of echoes, while e Ω causes the change of the Doppler frequency and is called the effective rotational vector. Considering the generic scatterer p with the directional vector p r , we represent its Doppler frequency as: where  and  represent the outer product and inner product, respectively.  denotes the wavelength of the transmitted signal. The effective rotational vector e Ω is usually time-varying for the ship target. Based on the Weierstrass approximation principle [29] and analyses in [12], e Ω of the ship target can be approximated as: α β γ (2) where m t is the slow time (or azimuth time). α , β and γ denote coefficients of the constant-, first-, and second-term of e Ω , respectively (or the angular velocity, acceleration and acceleration rate, respectively). Substituting (2) into (1), we have: With (3), after the motion compensation with the standard range alignment method, keystone transform and phase gradient autofocus method, the azimuth echo of the generic scatterer p [18,19] can be represented as:  Considering the generic scatterer p with the directional vector r p , we represent its Doppler frequency as: where × and • represent the outer product and inner product, respectively. λ denotes the wavelength of the transmitted signal. The effective rotational vector Ω e is usually time-varying for the ship target. Based on the Weierstrass approximation principle [29] and analyses in [12], Ω e of the ship target can be approximated as: where t m is the slow time (or azimuth time). α, β and γ denote coefficients of the constant-, first-, and second-term of Ω e , respectively (or the angular velocity, acceleration and acceleration rate, respectively). Substituting (2) into (1), we have: With (3), after the motion compensation with the standard range alignment method, keystone transform and phase gradient autofocus method, the azimuth echo of the generic scatterer p [18,19] can be represented as: where σ p denotes the amplitude. φ 1,p = −2 r p × α •R/λ, φ 2,p = −2 r p × β •R/λ and φ 3,p = −2 r p × γ •R/λ denote the centroid frequency, chirp rate and quadratic chirp rate, respectively. The azimuth echo of the generic scatterer p takes the form of the CPS. In ISAR imaging, multicomponent scatterers coexist each range cell. Assuming the number of scatterers in the i-th (1 ≤ i ≤ I) range cell is P and taking the additive complex white Gaussian noise n i (t m ) into account, we can represent azimuth echoes in the i-th range cell as noisy multicomponent CPSs: Obviously, scatterers at different coordinates correspond to different centroid frequencies in each range cell and we can use this characteristic to construct ISAR image of the ship target [3]. Nevertheless, the chirp rate and quadratic chirp rate also exist in (5) and will induce the Doppler spread (degrade the cross-range resolution). Thus, now, the task is to estimate these two parameters and compensate the corresponding Doppler spread.

ICPBAF for Multicomponent CPSs
In order to construct a well-focused ISAR image of the ship target, the ICPBAF is proposed for the parameters estimation of the CPS. The principle of the ICPBAF and its cross term characteristic are discussed in this section.

Principle of the ICPBAF
According to the reference [24], the CPBAF of (5) is presented as: where τ m denotes the lag variable. R i,cros (t m , τ m ) and n R,i (t m , τ m ) denote the cross term and noise, respectively. In the auto term (corresponding to the CPS) of the CPBAF, the slow time t m and lag τ m nonlinearly couple with each other. If the first exponent of the auto term does not exist, we can employ the decoupling technique GDT to eliminate the coupling and then, accumulate the signal energy coherently with the Fourier transform along the τ m axis. That is, the first exponent of the auto term disables the GDT and the following Fourier transform-based coherent accumulation. It is easily seen from (6) that, due to the coupling, the Fourier transform along the τ m axis can accumulate the signal energy into the inclined line, which benefits the signal detection and parameters estimation [28]. Here, we employ this characteristic and perform the Fourier transform along the τ m axis. However, the τ m axis corresponds to the nonuniform τ 2 m and the FFT is no longer applicable. Fortunately, we can adopt the NUFFT to speed up the Fourier transform of the nonuniform data without the performance loss. Readers can refer to references [15,16] for more details about the NUFFT: where f τ 2 m is the frequency domain with respect to τ m . δ(•) is the Dirac delta function. G i,cros t m , f τ 2 m and n G,i t m , f τ 2 m denote the cross tern and noise after the NUFFT, respectively. In (7), the energy of the auto term peaks along the inclined line f τ 2 m − φ 2,p − φ 3,p t m = 0. In general, we can accumulate the signal energy along this line by employing the Radon transform or Hough transform. However, this kind of processing methods is incoherent and the two-dimensionally brute-force searching will influence the efficiency [8,10,12]. In order to accumulate the signal energy coherently, we take the complex modulus of G i t m , f τ 2 m and perform the IFFT with respect to f τ 2 m : where b m is the time domain with respect to f τ 2 m .|•| denotes the operation of taking the complex modulus. Q i,cros (t m , b m ) and n Q,i (t m , b m ) denote the cross term and noise, respectively. In realistic applications, the signal length is finite and the sinc function should be used instead of the Dirac function. Under this condition, taking the modulus bends the negative side lobes to the positive and taking the IFFT may induce a different equation than (8). However, the energy of the sinc function concentrates in the main lobe and the effectiveness of negative side lobes is very small. Therefore, to be exact, we use the approximation ≈ in (8).
In (8), the operation of taking the complex modulus eliminates the disturbance back into the form of the exponent. Obviously, we can employ the GDT to eliminate the linear coupling between t m and b m in the auto term of Q i (t m , b m ). According to the reference [22], the GDT takes the form: is the scaled frequency domain with respect to Υ m . b is a known constant. ξ denotes the zoom factor. φ is an unknown parameter.
Based on the coupling in (8) and the GDT in (9), several substitutions are done as follows: With substitutions in (10), we apply the GDT to Q i (t m , b m ): The GDT in (11) can be implemented by using the IFFT-based CZT. More details about the fast implementation of the GDT can be found in the reference [22]. After the GDT, the energy of the auto term is accumulated into the beeline f [b m t m ] − φ 3,p = 0. Now, we perform the FFT along the b m axis and obtain the ICPBAF: where f b m is the frequency domain with respect to b m .
, f b m denote the cross term and noise after the FFT, respectively. For the ICPBAF, the auto term peaks at (φ 2,p ,φ 3,p ) on the chirp rate-quadratic chirp rate domain. Parameters φ 2,p and φ 3,p can be estimated by constructing a cost function to (12) [23]. With these two estimated parameters, other parameters (σ p and φ 1,p ) can be obtained by the dechirping and FFT operations [23].
According to (6)-(8), (11) and (12), we give the abbreviated expression of the ICPBAF: Analyses above indicate that the ICPBAF is a coherent bilinear algorithm and can be implemented by using the complex multiplication, FFT, IFFT and NUFFT. The two-dimensionally coherent accumulation and bilinearity guarantee the low computational cost, high anti-noise performance and cross term suppression of the ICPBAF, which will be demonstrated in Section 4. The cross term characteristic analysis can demonstrate whether the cross term can accumulate as the auto term or not [10,12]. Further, with the cross term characteristic analysis, we can also have a more in-depth understanding of the cross term suppression. Therefore, we analyze the cross term characteristic of the ICPBAF in the following section.

Cross Term Characteristic Analysis
In order to formulate the cross term problem arising from multicomponent CPSs, we consider the noise-free signal with two components, l ∈ [1, P − 1] and q ∈ [l + 1, P], which is denoted as: The CPBAF of s l,q (t m ) can be presented as: where: Performing the NUFFT, operation of taking the complex modulus, IFFT, GDT and FFT on R l,q (t m , τ m ), we have: With careful analyses of Equations (15a) and (15b), we find that, compared to the auto term R l,q,auto (t m , τ m ), the cosine function of R l,q,cros (t m , τ m ) will disturb the NUFFT, operation of taking the complex modulus, IFFT, GDT and FFT. According to the characteristic of the cosine function, when the Equation (17) is satisfied for every t m and τ m , the influence of the cosine function will be eliminated: Obviously, only when φ 1,l = φ 1,q , φ 2,l = φ 2,q and φ 3,l = φ 3,q , Equation (17) equals to zero for every t m and τ m . That is, different from the auto term, the cross term of the ICPBAF cannot be accumulated. Thus, when the signal length is infinite, compared to the auto term, the cross term can be ignored, i.e., the ICPBAF in (12) can be approximated as: This section gives the basic analyses of the ICPBAF, including the principle of the ICPBAF and its cross term characteristic analysis. In order to illustrate how the ICPBAF works under multicomponent CPSs, we give a numerical example in the following. Considering realistic applications, this numerical example includes two situations, multicomponent CPSs with the same amplitude and multicomponent CPSs with the different amplitudes.

Performance Analysis of the ICPBAF
The computational cost, cross term suppression and anti-noise performance play important roles in the parameters estimation algorithm [22][23][24]. In this section, we analyze the ICPBAF from these three aspects, and some comparisons with other five representative algorithms including the CPF, IGCPF, SCFT-based algorithm, GSCFT-based algorithm and GDT-based algorithm will be performed.

Computational Cost Analysis
Assume the length of slow time samples  Table 1, which also gives computational costs of five other representative algorithms, the CPF, IGCPF, SCFT-based algorithm, GSCFT-based algorithm and GDT-based algorithm.
Referring to references [22,23], we know that the CPF only uses the discrete Fourier transform to accumulate the CPBAF along m  and discards the energy along m t . Compared to the ICPBAF, its

Performance Analysis of the ICPBAF
The computational cost, cross term suppression and anti-noise performance play important roles in the parameters estimation algorithm [22][23][24]. In this section, we analyze the ICPBAF from these three aspects, and some comparisons with other five representative algorithms including the CPF, IGCPF, SCFT-based algorithm, GSCFT-based algorithm and GDT-based algorithm will be performed.

Computational Cost Analysis
Assume the length of slow time samples N t m is equal to the length of lag samples N τ m . The ICPBAF implementation needs the CPBAF [O N 2 Table 1. Computational cost.
Referring to references [22,23], we know that the CPF only uses the discrete Fourier transform to accumulate the CPBAF along τ m and discards the energy along t m . Compared to the ICPBAF, its low computational cost is at the cost of the low cross term suppression and anti-noise performance. The IGCPF and SCFT-based algorithm need higher computational costs than the GSCFT-based algorithm, GDT-based algorithm and ICPBAF. This is because the Fourier transform along the non-uniform τ m axis of these two algorithms are not speed up [13,19]. The GSCFT-based algorithm, GDT-based algorithm and ICPBAF need the similar computational cost. In the following two sections, we will demonstrate that the ICPBAF has superiorities in the cross term suppression and anti-noise performance. We can refer to references [12,13,19,[22][23][24] for more details about these five referenced algorithms.

Cross Term Suppression Analysis
Analyses in Section 3.2 demonstrate that the cross term cannot accumulate as the auto term. In this section, we uses the numerical example below to demonstrate the high cross term suppression of the ICPBAF.

Example 2.
We consider three noise-free CPSs denoted by Bu1, Bu2 and Bu3. F t m and N t m are 128 Hz and 256, respectively. The signal parameters are set as follows: σ 1 = 1, φ 1,1 = 10 Hz, φ 2,1 = 2 Hz/s, φ 3,1 = 2 Hz/s 2 for Bu1; Figure 3 gives simulation results of the CPF, IGCPF, SCFT-based algorithm, GSCFT-based algorithm, GDT-based algorithm and ICPBAF. low computational cost is at the cost of the low cross term suppression and anti-noise performance. The IGCPF and SCFT-based algorithm need higher computational costs than the GSCFT-based algorithm, GDT-based algorithm and ICPBAF. This is because the Fourier transform along the nonuniform m  axis of these two algorithms are not speed up [13,19]. The GSCFT-based algorithm, GDT-based algorithm and ICPBAF need the similar computational cost. In the following two sections, we will demonstrate that the ICPBAF has superiorities in the cross term suppression and anti-noise performance. We can refer to references [12,13,19,[22][23][24] for more details about these five referenced algorithms.

Cross Term Suppression Analysis
Analyses in Section 3.2 demonstrate that the cross term cannot accumulate as the auto term. In this section, we uses the numerical example below to demonstrate the high cross term suppression of the ICPBAF.
Hz/s 2 for Bu3. Figure 3 gives simulation results of the CPF, IGCPF, SCFT-based algorithm, GSCFT-based algorithm, GDT-based algorithm and ICPBAF. Obviously, in Figure 3, only the ICPBAF can obtain correct positions of Bu1, Bu2 and Bu3. The CPF is bilinear, while it discards the energy along m t . Although the IGCPF, SCFT-based algorithm, GSCFT-based algorithm and GDT-based algorithm employ the two-dimensional energy accumulation to suppress the cross term, their fourth-order autocorrelation functions influence the cross term suppression seriously. In this example, their fourth-order autocorrelation functions induce the cross term with 78 components. Note that, under some special situations, the cross terms of the CPF and IGCPF can accumulate as their auto terms [19,30]. Compared to the CPF, IGCPF, SCFT- Obviously, in Figure 3, only the ICPBAF can obtain correct positions of Bu1, Bu2 and Bu3. The CPF is bilinear, while it discards the energy along t m . Although the IGCPF, SCFT-based algorithm, GSCFT-based algorithm and GDT-based algorithm employ the two-dimensional energy accumulation to suppress the cross term, their fourth-order autocorrelation functions influence the cross term suppression seriously. In this example, their fourth-order autocorrelation functions induce the cross term with 78 components. Note that, under some special situations, the cross terms of the CPF and IGCPF can accumulate as their auto terms [19,30]. Compared to the CPF, IGCPF, SCFT-based algorithm, GSCFT-based algorithm and GDT-based algorithm, the ICPBAF is bilinear and employs the two-dimensionally coherent energy accumulation. Therefore, as shown in Figure 3, the ICPBAF acquires the highest cross term suppression.

Anti-Noise Performance Analysis
The ICPBAF inevitably suffers from the estimation error under the presence of the noise. In this section, combing with a numerical example, we will demonstrate the high anti-noise performance of the ICPBAF. The input-output signal-to-noise ratio (SNR) (SNR out is listed in (19)) [31] and mean square error (MSE) [32,33] are utilized as measures of the noise resistance: where σ 2 is the power of the complex white Gaussian noise. a 2.p and a 3,p are estimations. based algorithm, GSCFT-based algorithm and GDT-based algorithm, the ICPBAF is bilinear and employs the two-dimensionally coherent energy accumulation. Therefore, as shown in Figure 3, the ICPBAF acquires the highest cross term suppression.

Anti-Noise Performance Analysis
The ICPBAF inevitably suffers from the estimation error under the presence of the noise. In this section, combing with a numerical example, we will demonstrate the high anti-noise performance of the ICPBAF. The input-output signal-to-noise ratio (SNR) ( out SNR is listed in (19)) [31] and mean square error (MSE) [32,33] are utilized as measures of the noise resistance:  In Figure 4, threshold SNRs of the ICPBAF, GDT-based algorithm, GSCFT-based algorithm, SCFT-based algorithm, IGCPF and CPF are −8 dB, −5 dB, −3 dB, −3 dB, −2 dB and −2 dB, respectively. The bilinear CPF discards the energy accumulation along m t and references [23,24] has simulated its threshold SNR. The IGCPF employs the two-dimensional energy accumulation, while the accumulation is incoherent and its autocorrelation function is fourth-order [19]. Thus, its threshold SNR is also −2 dB and no better than that of the CPF. Note that, the CPF and IGCPF estimate the chirp rate and quadratic chirp rate separately, and thus the propagation of the estimation error exists in these two algorithms. Although the GSCFT-based algorithm, SCFT-based algorithm and GDT-based algorithm employ the two-dimensionally coherent energy accumulation and do not encounter the estimation error propagation, their fourth-order autocorrelation functions influence the anti-noise performance seriously. In Figure 4, the superiority of the ICPBAF in the anti-noise performance is obvious. This is because (1) the ICPBAF is bilinear; (2) the estimation error propagation does not exist, and (3) the two-dimensional energy accumulation is coherent. Above analyses also apply to noisy multicomponent CPSs. We consider noisy CPSs with two components. The GDT-based algorithm, GSCFT-based algorithm, SCFT-based algorithm and IGCPF are based on fourth-order autocorrelation functions, and the number of noise items is 65. However, the ICPBAF is bilinear and the number of noise items is 5. Although the CPF is bilinear and the number of noise items is 5, CPF is based on one-dimensional energy accumulation. In Figure 4, threshold SNRs of the ICPBAF, GDT-based algorithm, GSCFT-based algorithm, SCFT-based algorithm, IGCPF and CPF are −8 dB, −5 dB, −3 dB, −3 dB, −2 dB and −2 dB, respectively. The bilinear CPF discards the energy accumulation along t m and references [23,24] has simulated its threshold SNR. The IGCPF employs the two-dimensional energy accumulation, while the accumulation is incoherent and its autocorrelation function is fourth-order [19]. Thus, its threshold SNR is also −2 dB and no better than that of the CPF. Note that, the CPF and IGCPF estimate the chirp rate and quadratic chirp rate separately, and thus the propagation of the estimation error exists in these two algorithms. Although the GSCFT-based algorithm, SCFT-based algorithm and GDT-based algorithm employ the two-dimensionally coherent energy accumulation and do not encounter the estimation error propagation, their fourth-order autocorrelation functions influence the anti-noise performance seriously. In Figure 4, the superiority of the ICPBAF in the anti-noise performance is obvious. This is because (1) the ICPBAF is bilinear; (2) the estimation error propagation does not exist, and (3) the two-dimensional energy accumulation is coherent. Above analyses also apply to noisy multicomponent CPSs. We consider noisy CPSs with two components. The GDT-based algorithm, GSCFT-based algorithm, SCFT-based algorithm and IGCPF are based on fourth-order autocorrelation functions, and the number of noise items is 65. However, the ICPBAF is bilinear and the number of noise items is 5. Although the CPF is bilinear and the number of noise items is 5, CPF is based on one-dimensional energy accumulation.
As functions of the input SNR, the observed MSEs for the chirp rate and quadratic chirp rate estimations are plotted in Figure 5a,b, respectively. The corresponding Cramer-Rao bounds (CRBs) are also shown in solid lines and their expressions can be found in [24,31,32].
As expected, the observed MSEs of the chirp rate and quadratic chirp rate estimations are inversely proportional to the input SNRs in Figure 5. MSEs of the chirp rate and the quadratic chirp rate estimations are close to CRB when SNR ≥ −8 dB. Results shown in Figure 5 verify the high anti-noise performance of the ICPBAF and also validate the result shown in Figure 4. As functions of the input SNR, the observed MSEs for the chirp rate and quadratic chirp rate estimations are plotted in Figures 5a,b, respectively. The corresponding Cramer-Rao bounds (CRBs) are also shown in solid lines and their expressions can be found in [24,31,32]. As expected, the observed MSEs of the chirp rate and quadratic chirp rate estimations are inversely proportional to the input SNRs in Figure 5. MSEs of the chirp rate and the quadratic chirp rate estimations are close to CRB when SNR ≥ −8 dB. Results shown in Figure 5 verify the high antinoise performance of the ICPBAF and also validate the result shown in Figure 4.

RID ISAR Imaging Algorithm Based on the ICPBAF
Analyses and simulations in Section 4 demonstrate that, compared to the CPF, IGCPF, SCFTbased algorithm, GSCFT-based algorithm and GDT-based algorithm, the ICPBAF is more suitable for the parameters estimation of noisy multicomponent CPSs. In this section, based on the ICPBAF, a RID ISAR imaging algorithm is presented for the ship target. Detailed implementation procedures are given as follows.
Step 1: Complete the range compression for radar echoes with the matched filter (where t , s T and  denote the fast time, pulse width and frequency modulation rate, respectively).
Step 2: Employ the standard range alignment method [5], keystone transform [6] and phase gradient autofocus method [7] to complete the motion compensation.
Step 3: Extract the data   i m s t in the i-th range cell.
Step 4: Calculate the energy of the extracted data . If the energy is smaller than the set threshold s E [8,12], set where ' p  and ' 1, p  denote estimations of the centroid frequency and amplitude for the pth CPS, respectively. ' A denotes the peak value after the FFT.
Step 7: Eliminate the estimated p-th CPS from the original signal  

RID ISAR Imaging Algorithm Based on the ICPBAF
Analyses and simulations in Section 4 demonstrate that, compared to the CPF, IGCPF, SCFT-based algorithm, GSCFT-based algorithm and GDT-based algorithm, the ICPBAF is more suitable for the parameters estimation of noisy multicomponent CPSs. In this section, based on the ICPBAF, a RID ISAR imaging algorithm is presented for the ship target. Detailed implementation procedures are given as follows.
Step 1: Complete the range compression for radar echoes with the matched filter H t = rect t /T s exp jπγt 2 (wheret, T s and γ denote the fast time, pulse width and frequency modulation rate, respectively).
Step 2: Employ the standard range alignment method [5], keystone transform [6] and phase gradient autofocus method [7] to complete the motion compensation.
Step 3: Extract the data s i (t m ) in the i-th range cell.
Step 4: Calculate the energy of the extracted data s i (t m ). If the energy is smaller than the set threshold E s [8,12], set i = i + 1 and repeat Step 3 until i = I.
Step 5: Apply the ICPBAF to s i (t m ) and estimate φ 2,p and φ 3,p with the peak detection technique [23,24].
Step 6: Dechirp s i (t m ) with exp −j2π where σ p and φ 1,p denote estimations of the centroid frequency and amplitude for the pth CPS, respectively. A denotes the peak value after the FFT.
Step 7: Eliminate the estimated p-th CPS from the original signal s i (t m ) denotes the narrowband filter with the bandwidth 2d (can be determined based on the resolution).
Step 8: Repeat Steps 5-7 until the residual signal energy E of the i-th range cell is less than E H (saying 5% of the original signal [18,19]), which is an energy threshold.
Step 9: If i < I, set i = i +1 and repeat Steps 3-8 until i = I.
denotes the narrowband filter with the bandwidth 2d (can be determined based on the resolution).
Step 8: Repeat Steps 5-7 until the residual signal energy E of the i-th range cell is less than H E (saying 5% of the original signal [18,19]), which is an energy threshold.
Step 9: If Above is the proposed RID ISAR imaging algorithm for the ship target based on the ICPBAF. With this algorithm, we can construct a well-focused ISAR image for the ship target. In order to provide insight into the working of the RID ISAR imaging algorithm, the flow chart of the proposed RID ISAR imaging algorithm is shown in Figure 6.

Verification of the ICPBAF-Based RID ISAR Imaging Algorithm
Section 4 verifies high performances of the ICPBAF and Section 5 presents a RID ISAR imaging algorithm based on the ICPBAF. In this section, we employ the synthetic model (Section 6.1) and real radar data (Section 6.2) to verify the superiority and practicability of the ICPBAF-based RID ISAR imaging algorithm.

Verification with the Synthetic Model
In this section, referring to [3,8,9,[17][18][19][20][21][22], we model the ship target shown in Figure 7a as a set of ideal scatterers. Table 2 gives radar and motion parameters. After the motion compensation, Figure  7b shows the image constructed with the standard RD technique under the absence of the scattererdependent Doppler spread, and Figure 7c shows the image constructed with the standard RD technique under the existence of the scatterer-dependent Doppler spread. Obviously, due to the scatterer-dependent Doppler spread, the image quality is degraded in Figure 7c. Note that, as described in Section 1, this paper focuses on the scatterer-dependent Doppler spread compensation and we can refer to [5][6][7] for more details about the motion compensation.

Verification of the ICPBAF-Based RID ISAR Imaging Algorithm
Section 4 verifies high performances of the ICPBAF and Section 5 presents a RID ISAR imaging algorithm based on the ICPBAF. In this section, we employ the synthetic model (Section 6.1) and real radar data (Section 6.2) to verify the superiority and practicability of the ICPBAF-based RID ISAR imaging algorithm.

Verification with the Synthetic Model
In this section, referring to [3,8,9,[17][18][19][20][21][22], we model the ship target shown in Figure 7a as a set of ideal scatterers. Table 2 gives radar and motion parameters. After the motion compensation, Figure 7b shows the image constructed with the standard RD technique under the absence of the scatterer-dependent Doppler spread, and Figure 7c shows the image constructed with the standard RD technique under the existence of the scatterer-dependent Doppler spread. Obviously, due to the scatterer-dependent Doppler spread, the image quality is degraded in Figure 7c. Note that, as described in Section 1, this paper focuses on the scatterer-dependent Doppler spread compensation and we can refer to [5][6][7] for more details about the motion compensation.  We contaminate echoes of the ship target with the additive complex white Gaussian noise and the SNRin equals to −5 dB after the motion compensation. Here, RID ISAR imaging algorithms, which use the CPF, IGCPF, SCFT-based algorithm, GSCFT-based algorithm and GDT-based algorithm, are adopted to compare with the ICPBAF-based RID ISAR imaging algorithm. Below, images constructed with these six RID ISAR imaging algorithms are normalized and shown in Figure 8. The entropy of (22) is used as a criterion to measure the quality of the image   , X i n in Table 3   We contaminate echoes of the ship target with the additive complex white Gaussian noise and the SNR in equals to −5 dB after the motion compensation. Here, RID ISAR imaging algorithms, which use the CPF, IGCPF, SCFT-based algorithm, GSCFT-based algorithm and GDT-based algorithm, are adopted to compare with the ICPBAF-based RID ISAR imaging algorithm. Below, images constructed with these six RID ISAR imaging algorithms are normalized and shown in Figure 8. The entropy of (22) is used as a criterion to measure the quality of the image X(i, n) in Table 3 [12]:  We contaminate echoes of the ship target with the additive complex white Gaussian noise and the SNRin equals to −5 dB after the motion compensation. Here, RID ISAR imaging algorithms, which use the CPF, IGCPF, SCFT-based algorithm, GSCFT-based algorithm and GDT-based algorithm, are adopted to compare with the ICPBAF-based RID ISAR imaging algorithm. Below, images constructed with these six RID ISAR imaging algorithms are normalized and shown in Figure 8. The entropy of (22) is used as a criterion to measure the quality of the image   , X i n in Table 3 Table 3. Entropies of ISAR images in Figure 8. The superiority of the ICPBAF-based RID ISAR imaging algorithm is obvious in Figure 8. All scatterers are reconstructed correctly and very few spurious scatterers appear. This is because the ICPBAF has the highest cross term suppression and anti-noise performance, which has been analyzed and demonstrated in Section 4. The better focused ISAR image results in the smaller entropy [12,13,[18][19][20][21][22]. Results in Table 3 demonstrate the high quality of Figure 8f and validate the superiority of the ICPBAF-based RID ISAR imaging algorithm also.
Actually, the additive complex white Gaussian noise is random. Thus, only with the experiment above, it may not be persuasive to determine the superiority of the ICPBAF-based RID ISAR imaging algorithm. It is known that the result of the Monte Carlo experiment has the generality and representability. Thus, here, a Monte Carlo experiment is used to determine the superiority of the ICPBAF-based RID ISAR imaging algorithm. The data in the 15th range cell, which exists six scatterers and is marked with the red ellipse in Figure 7a, is extracted. We contaminate the extracted data with the additive complex white Gaussian noise. The tested input SNRs are SNRin = [−11:1:0] and 100 trials are performed for each SNRin. Aforementioned six RID ISAR imaging algorithms are adopted to perform on the extracted data. In a well-focused ISAR image, most scatterers can be reconstructed correctly and fewest artifacts appear. Thus, the ratio between the number of the correctly reconstructed scatterers and the number of all reconstructed scatterers is used as a measure. Figure 9 shows the simulation result.
The high cross term suppression and anti-noise performance of the ICPBAF guarantee the superiority of the ICPBAF-based RID ISAR imaging algorithm. In Figure 9, even under lower SNRs, the ICPBAF-based RID ISAR imaging algorithm can reconstruct most scatterers correctly and fewest artifacts appear. For example, when SNRin = [−7 dB, −6 dB, −5 dB], the ratio equals to [0.7747, 0.8739, 0.944], which means that the average number of the artifacts is smaller than [1.7449, 0.8658, 0.3559]. This experiment further demonstrates that, compared to RID ISAR imaging algorithms using the CPF, IGCPF, SCFT-based algorithm, GSCFT-based algorithm and GDT-based algorithm, the ICPBAF-based RID ISAR imaging algorithm is more suitable for realistic applications.  Table 3. Entropies of ISAR images in Figure 8. The superiority of the ICPBAF-based RID ISAR imaging algorithm is obvious in Figure 8. All scatterers are reconstructed correctly and very few spurious scatterers appear. This is because the ICPBAF has the highest cross term suppression and anti-noise performance, which has been analyzed and demonstrated in Section 4. The better focused ISAR image results in the smaller entropy [12,13,[18][19][20][21][22]. Results in Table 3 demonstrate the high quality of Figure 8f and validate the superiority of the ICPBAF-based RID ISAR imaging algorithm also.
Actually, the additive complex white Gaussian noise is random. Thus, only with the experiment above, it may not be persuasive to determine the superiority of the ICPBAF-based RID ISAR imaging algorithm. It is known that the result of the Monte Carlo experiment has the generality and representability. Thus, here, a Monte Carlo experiment is used to determine the superiority of the ICPBAF-based RID ISAR imaging algorithm. The data in the 15th range cell, which exists six scatterers and is marked with the red ellipse in Figure 7a, is extracted. We contaminate the extracted data with the additive complex white Gaussian noise. The tested input SNRs are SNR in = [−11:1:0] and 100 trials are performed for each SNR in . Aforementioned six RID ISAR imaging algorithms are adopted to perform on the extracted data. In a well-focused ISAR image, most scatterers can be reconstructed correctly and fewest artifacts appear. Thus, the ratio between the number of the correctly reconstructed scatterers and the number of all reconstructed scatterers is used as a measure. Figure 9 shows the simulation result.
The high cross term suppression and anti-noise performance of the ICPBAF guarantee the superiority of the ICPBAF-based RID ISAR imaging algorithm. In Figure 9, even under lower SNRs, the ICPBAF-based RID ISAR imaging algorithm can reconstruct most scatterers correctly and fewest artifacts appear. For example, when SNR in = [−7 dB, −6 dB, −5 dB], the ratio equals to [0.7747, 0.8739, 0.944], which means that the average number of the artifacts is smaller than [1.7449, 0.8658, 0.3559]. This experiment further demonstrates that, compared to RID ISAR imaging algorithms using the CPF, IGCPF, SCFT-based algorithm, GSCFT-based algorithm and GDT-based algorithm, the ICPBAF-based RID ISAR imaging algorithm is more suitable for realistic applications.

Verification with Real Radar Data
The real radar data used here is received by a shore-based radar, which works in Ku band with a bandwidth of 240 MHz and a PRF of 125 Hz. The imaged ship target is moving away from the shore-based radar with a velocity of about 10 m/s in the harsh ocean environment. In the real radar data, the number of echo pulses m t N is 250 and the number of the slant range cells is 400. Figure 10a gives the result after the motion compensation with the standard range alignment method, keystone transform and phase gradient autofocus method. The Wigner-Ville distribution of the 191th range cell is given in Figure 10b. According to analyses and simulations in [12,13,21,22], the curve in Figure  10b demonstrates that azimuth echoes of the ship target should be modeled as noisy multicomponent CPSs. In order to give an unequivocal evidence that the ICPBAF-based RID ISAR imaging algorithm benefits the imaging quality, we use it to process the extracted data in the 191th range cell. Figure 11 gives the simulation results. The signal energy is accumulated in Figure 11a. With the peak detection technique, the chirp rate and quadratic chirp rate are estimated as −1 Hz/s and −8 Hz/s 2 , respectively. By compensating the Doppler spread pertaining to the estimated parameters and performing an FFT, we complete the energy accumulation in Figure 11b, where the result of the standard RD technique is also shown. Obviously, due to the Doppler spread, the standard RD technique cannot focus the signal energy into the correct Doppler cell. Actually, several scatterers may exist in this range cell. Thus, combing with the Clean technique, we still need to relocate other potential scatterers of this range cell.

Verification with Real Radar Data
The real radar data used here is received by a shore-based radar, which works in Ku band with a bandwidth of 240 MHz and a PRF of 125 Hz. The imaged ship target is moving away from the shore-based radar with a velocity of about 10 m/s in the harsh ocean environment. In the real radar data, the number of echo pulses N t m is 250 and the number of the slant range cells is 400. Figure 10a gives the result after the motion compensation with the standard range alignment method, keystone transform and phase gradient autofocus method. The Wigner-Ville distribution of the 191th range cell is given in Figure 10b. According to analyses and simulations in [12,13,21,22], the curve in Figure 10b demonstrates that azimuth echoes of the ship target should be modeled as noisy multicomponent CPSs.

Verification with Real Radar Data
The real radar data used here is received by a shore-based radar, which works in Ku band with a bandwidth of 240 MHz and a PRF of 125 Hz. The imaged ship target is moving away from the shore-based radar with a velocity of about 10 m/s in the harsh ocean environment. In the real radar data, the number of echo pulses m t N is 250 and the number of the slant range cells is 400. Figure 10a gives the result after the motion compensation with the standard range alignment method, keystone transform and phase gradient autofocus method. The Wigner-Ville distribution of the 191th range cell is given in Figure 10b. According to analyses and simulations in [12,13,21,22], the curve in Figure  10b demonstrates that azimuth echoes of the ship target should be modeled as noisy multicomponent CPSs. In order to give an unequivocal evidence that the ICPBAF-based RID ISAR imaging algorithm benefits the imaging quality, we use it to process the extracted data in the 191th range cell. Figure 11 gives the simulation results. The signal energy is accumulated in Figure 11a. With the peak detection technique, the chirp rate and quadratic chirp rate are estimated as −1 Hz/s and −8 Hz/s 2 , respectively. By compensating the Doppler spread pertaining to the estimated parameters and performing an FFT, we complete the energy accumulation in Figure 11b, where the result of the standard RD technique is also shown. Obviously, due to the Doppler spread, the standard RD technique cannot focus the signal energy into the correct Doppler cell. Actually, several scatterers may exist in this range cell. Thus, combing with the Clean technique, we still need to relocate other potential scatterers of this range cell. In order to give an unequivocal evidence that the ICPBAF-based RID ISAR imaging algorithm benefits the imaging quality, we use it to process the extracted data in the 191th range cell. Figure 11 gives the simulation results. The signal energy is accumulated in Figure 11a. With the peak detection technique, the chirp rate and quadratic chirp rate are estimated as −1 Hz/s and −8 Hz/s 2 , respectively. By compensating the Doppler spread pertaining to the estimated parameters and performing an FFT, we complete the energy accumulation in Figure 11b, where the result of the standard RD technique is also shown. Obviously, due to the Doppler spread, the standard RD technique cannot focus the signal energy into the correct Doppler cell. Actually, several scatterers may exist in this range cell. Thus, combing with the Clean technique, we still need to relocate other potential scatterers of this range cell.   In Figure 12, RID ISAR imaging algorithms, which use the CPF, IGCPF, SCFT-based algorithm, GSCFT-based algorithm and GDT-based algorithm, are adopted to compare with the ICPBAF-based RID ISAR imaging algorithm. The result of the standard RD technique is also shown. These images are normalized and corresponding entropies are given in Table 4. Table 4. Entropies OF ISAR images in Figure 12.  Figure 12a is the result of the standard RD technique. Obviously, due to the Doppler spread induced by the chirp rate and quadratic chirp rate, it cannot construct a well-focused image for the ship target. The CPF, IGCPF, SCFT-based algorithm, GSCFT-based algorithm, GDT-based algorithm and ICPBAF can estimate parameters of noisy multicomponent CPSs. Thus, images shown in Figure 12b-g are better than the image shown in Figure 12a, which can also be demonstrated with the entropies listed in Table 4. In Figure 12d-g, it is easy to find that, compared to the SCFT-based algorithm, GSCFT-based algorithm and GDT-based algorithm, the advantage of the ICPBAF-based algorithm is not so obvious. This is because the realistic environment is unknown and we cannot control characteristics of the real radar data, such as the SNR and the number of scatterers in each range cell. Actually, with simulations and analyses in Sections 4 and 6.1, we know that, if a harsher realistic environment is considered, the advantage of the ICPBAF-based algorithm will be much more obvious. Images shown in Figure 12 and entropies listed in Table 4 verify the practicability of the ICPBAF-based RID ISAR imaging algorithm.

Conclusions
In this paper, a bilinear coherent estimation algorithm, known as the ICPBAF, is proposed for the CPS. The bilinear ICPBAF can complete the two-dimensionally coherent energy accumulation. The principle, cross term characteristic, computational cost, cross term suppression and anti-noise performance are analyzed for the ICPBAF. Comparisons with five other representative estimation algorithms demonstrate that the ICPBAF can acquire the higher cross term suppression and anti-noise performance with a moderate computational cost. It is worthwhile noting that, for the CPS, the ICPBAF can be seen as the first bilinear coherent algorithm, which has a high practicability. Thereafter, the ICPBAF-based RID ISAR imaging algorithm is presented for ship targets, and we use the synthetic model and real radar data to verify its effectiveness.