Home | About | Journals | Submit | Contact Us | Français |

**|**Sensors (Basel)**|**v.17(3); 2017 March**|**PMC5375784

Formats

Article sections

- Abstract
- 1. Introduction
- 2. ISAR Imaging Geometry for the Ship Target
- 3. ICPBAF for Multicomponent CPSs
- 4. Performance Analysis of the ICPBAF
- 5. RID ISAR Imaging Algorithm Based on the ICPBAF
- 6. Verification of the ICPBAF-Based RID ISAR Imaging Algorithm
- 7. Conclusions
- References

Authors

Related links

Sensors (Basel). 2017 March; 17(3): 498.

Published online 2017 March 3. doi: 10.3390/s17030498

PMCID: PMC5375784

Ram M. Narayanan, Academic Editor

Received 2017 January 20; Accepted 2017 February 24.

Copyright © 2017 by the authors.

Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

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.

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.

In general, after the motion compensation, azimuth echoes in a range cell of the ship target can be modeled as noisy multicomponent linear frequency-modulated (LFM) signals [8,9,10]. Nevertheless, in the harsh ocean environment, the LFM signal model is not suitable for ship targets. In references [12,13], simulations and analyses with the real radar data demonstrate that, for the ship targets in the harsh ocean environment, azimuth echoes in a range cell have to be modeled as noisy multicomponent cubic phase signals (CPSs). For ISAR imaging with the CPS model, the chirp rate and quadratic chirp rate induce the Doppler spread and defocus the ISAR image [14]. Therefore, now, the task of the RID technique is to estimate these two parameters and compensate the Doppler spread.

Until now, many researchers have studied the CPS and proposed several successful estimation algorithms. These estimation algorithms can be generally divided into three categories: linear algorithms [14], multilinear algorithms [12,13,17,18,19,20,21,22] and bilinear algorithms [23,24,25,26]. The linear algorithms, such as the modified discrete chirp Fourier transform for the CPS [14], employ three-dimensional brute-force searching to obtain a high anti-noise performance without the cross term, but this is at the cost of a high computational burden [10,12]. Multilinear algorithms employ the fourth-order autocorrelation function to reduce the phase order of the CPS, and then complete the two-dimensional energy accumulation with several operations, such as the fast Fourier transform (FFT), non-uniform FFT (NUFFT) [15,16] and decoupling techniques. Compared with linear algorithms, the multilinear algorithms have a lower computational cost [12]. This is the reason why many researchers currently focus on them. Representative multilinear algorithms include the product high-order ambiguity function [17], product generalized cubic phase function [18], integrated generalized cubic phase function (IGCPF) [19], scaled Fourier transform (SCFT)-based algorithm [13], noise-resistant parameter estimation algorithm [20], coherently IGCPF [21], generalized SCFT (GSCFT)-based algorithm [12] and generalized decoupling technique (GDT)-based algorithm [22]. However, the fourth-order autocorrelation function influences the cross term suppression and anti-noise performance seriously. Drawing lessons from LFM signal research, we know that the bilinear algorithm, which is based on the two-dimensionally coherent energy accumulation, can resolve the problems of the linear and multilinear algorithms [8,9,10]. Several bilinear algorithms have been developed for the CPS, such as the cubic phase function (CPF) [23,24], non-uniform sampled CPF [25] and local polynomial Wigner distribution [26]. Nevertheless, due to the quadratic chirp rate, these bilinear algorithms can only accumulate the energy along the delay axis and discard the energy along the slow time axis. Thus, spurious peaks can be easily induced and the anti-noise performance is poor. Several researchers have employed two-dimensional brute-force searching to improve these bilinear algorithms [27], although the two-dimensional brute-force searching makes the bilinear algorithm lose its inherent advantages [8]. In [28], we used the parameter space switching method to speed up the accumulation of the auto term. However, the accumulation is incoherent, and the resolution, cross term suppression and anti-noise performance are still poor.

In this paper, by employing the cubic phase bilinear autocorrelation function (CPBAF), NUFFT, operation of taking the complex modulus, inverse FFT (IFFT), GDT and FFT, a novel estimation algorithm, known as the integrated CPBAF (ICPBAF) is proposed for the CPS. The ICPBAF is bilinear and based on the two-dimensionally coherent energy accumulation. The bilinearity and two-dimensionally coherent accumulation of the ICPBAF guarantee the high cross term suppression and anti-noise performance. Compared to five other representative estimation algorithms including the CPF, IGCPF, SCFT-based algorithm, GSCFT-based algorithm and GDT-based algorithm, the ICPBAF can acquire higher cross term suppression and anti-noise performance with a moderate computational cost. Thereafter, with the ICPBAF, we present a RID ISAR imaging algorithm for the ship target. Through simulations and analyses with the synthetic model and the real radar data, we verify the effectiveness of the ICPBAF and corresponding RID ISAR imaging algorithm.

According to references [5,6,7,12,13], the ISAR imaging geometry described here is based on the assumption that the motion compensation has been implemented and only the scatterer-dependent Doppler spread is considered.

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). ${\mathbf{\Omega}}_{\mathbf{roll}}$, ${\mathbf{\Omega}}_{\mathbf{pitch}}$ and ${\mathbf{\Omega}}_{\mathbf{yaw}}$ denote the angular rotational vectors of the ship target rotating around the *X*, *Y* and *Z* axes, respectively. $\mathbf{\Omega}$ is the synthetic vector of ${\mathbf{\Omega}}_{\mathbf{roll}}$, ${\mathbf{\Omega}}_{\mathbf{pitch}}$ and ${\mathbf{\Omega}}_{\mathbf{yaw}}$, and can be decomposed into ${\mathbf{\Omega}}_{\mathbf{e}}$ and ${\mathbf{\Omega}}_{\mathbf{R}}$ in the plane determined by **R** and $\mathbf{\Omega}$. ${\mathbf{\Omega}}_{\mathbf{R}}$ does not cause the radial motion and thus has no effect on the phase of echoes, while ${\mathbf{\Omega}}_{\mathbf{e}}$ causes the change of the Doppler frequency and is called the effective rotational vector.

Considering the generic scatterer *p* with the directional vector ${\mathbf{r}}_{\mathbf{p}}$, we represent its Doppler frequency as:

$${f}_{D}=-\frac{2}{\lambda}\left[\left({\mathbf{r}}_{\mathbf{p}}\times {\mathbf{\Omega}}_{\mathbf{e}}\right)\u2022\mathbf{R}\right]$$

(1)

where × and • represent the outer product and inner product, respectively. $\lambda $ denotes the wavelength of the transmitted signal.

The effective rotational vector ${\mathbf{\Omega}}_{\mathbf{e}}$ is usually time-varying for the ship target. Based on the Weierstrass approximation principle [29] and analyses in [12], ${\mathbf{\Omega}}_{\mathbf{e}}$ of the ship target can be approximated as:

$${\mathbf{\Omega}}_{\mathbf{e}}\left({t}_{m}\right)=\mathsf{\alpha}+\mathsf{\beta}{t}_{m}+\frac{1}{2}\mathsf{\gamma}{t}_{m}^{2}$$

(2)

where ${t}_{m}$ is the slow time (or azimuth time). $\mathsf{\alpha}$, $\mathsf{\beta}$ and $\mathsf{\gamma}$ denote coefficients of the constant-, first-, and second-term of ${\mathbf{\Omega}}_{\mathbf{e}}$, respectively (or the angular velocity, acceleration and acceleration rate, respectively).

Substituting (2) into (1), we have:

$${f}_{D}=-\frac{2}{\lambda}\left\{\left({\mathbf{r}}_{\mathbf{p}}\times \mathsf{\alpha}\right)\u2022\mathbf{R}+\left[\left({\mathbf{r}}_{\mathbf{p}}\times \mathsf{\beta}\right)\u2022\mathbf{R}\right]{t}_{m}+\frac{1}{2}\left[\left({\mathbf{r}}_{\mathbf{p}}\times \mathsf{\gamma}\right)\u2022\mathbf{R}\right]{t}_{m}^{2}\right\}$$

(3)

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:

$${s}_{p}\left({t}_{m}\right)={\sigma}_{p}\mathrm{exp}\left[j2\pi \left({\varphi}_{1,p}{t}_{m}+\frac{1}{2}{\varphi}_{2,p}{t}_{m}^{2}+\frac{1}{6}{\varphi}_{3,p}{t}_{m}^{3}\right)\right]$$

(4)

where ${\sigma}_{p}$ denotes the amplitude. ${\varphi}_{1,p}=-2\left({\mathbf{r}}_{\mathbf{p}}\times \mathsf{\alpha}\right)\u2022\mathbf{R}/\lambda $, ${\varphi}_{2,p}=-2\left({\mathbf{r}}_{\mathbf{p}}\times \mathsf{\beta}\right)\u2022\mathbf{R}/\lambda $ and ${\varphi}_{3,p}=-2\left({\mathbf{r}}_{\mathbf{p}}\times \mathsf{\gamma}\right)\u2022\mathbf{R}/\lambda $ 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\le i\le \mathrm{I}$) range cell is P and taking the additive complex white Gaussian noise ${n}_{i}\left({t}_{m}\right)$ into account, we can represent azimuth echoes in the *i*-th range cell as noisy multicomponent CPSs:

$${s}_{i}\left({t}_{m}\right)={\displaystyle \sum _{p=1}^{\mathrm{P}}{s}_{p}\left({t}_{m}\right)}+{n}_{i}\left({t}_{m}\right)={\displaystyle \sum _{p=1}^{\mathrm{P}}{\sigma}_{p}\mathrm{exp}\left[j2\pi \left({\varphi}_{1,p}{t}_{m}+\frac{1}{2}{\varphi}_{2,p}{t}_{m}^{2}+\frac{1}{6}{\varphi}_{3,p}{t}_{m}^{3}\right)\right]}+{n}_{i}\left({t}_{m}\right)$$

(5)

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.

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.

According to the reference [24], the CPBAF of (5) is presented as:

$$\begin{array}{ll}{R}_{i}\left({t}_{m},{\tau}_{m}\right)& ={s}_{i}\left({t}_{m}-{\tau}_{m}\right){s}_{i}\left({t}_{m}+{\tau}_{m}\right)\\ & =\underset{\text{the auto tern}}{\underbrace{{\displaystyle \sum _{p=1}^{\mathrm{P}}{\sigma}_{p}^{2}\mathrm{exp}\left[j2\pi \left(2{\varphi}_{1,p}{t}_{m}+{\varphi}_{2,p}{t}_{m}^{2}+\frac{1}{3}{\varphi}_{3,p}{t}_{m}^{3}\right)\right]\mathrm{exp}\left[j2\pi \left({\varphi}_{2,p}+{\varphi}_{3,p}{t}_{m}\right){\tau}_{m}^{2}\right]}}}+{R}_{i,cros}\left({t}_{m},{\tau}_{m}\right)+{n}_{R,i}\left({t}_{m},{\tau}_{m}\right)\end{array}$$

(6)

where ${\tau}_{m}$ denotes the lag variable. ${R}_{i,cros}\left({t}_{m},{\tau}_{m}\right)$ and ${n}_{R,i}\left({t}_{m},{\tau}_{m}\right)$ 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 ${\tau}_{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 ${\tau}_{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 ${\tau}_{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 ${\tau}_{m}$ axis. However, the ${\tau}_{m}$ axis corresponds to the nonuniform ${\tau}_{m}^{2}$ 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:

$$\begin{array}{ll}{G}_{i}\left({t}_{m},{f}_{{\tau}_{m}^{2}}\right)& ={\mathrm{NUFFT}}_{{\tau}_{m}^{2}}\left[{R}_{i}\left({t}_{m},{\tau}_{m}\right)\right]\\ & =\underset{\text{the auto tern}}{\underbrace{{\displaystyle \sum _{p=1}^{\mathrm{P}}{\sigma}_{p}^{2}\mathrm{exp}\left[j2\pi \left(2{\varphi}_{1,p}{t}_{m}+{\varphi}_{2,p}{t}_{m}^{2}+\frac{1}{3}{\varphi}_{3,p}{t}_{m}^{3}\right)\right]\delta \left({f}_{{\tau}_{m}^{2}}-{\varphi}_{2,p}-{\varphi}_{3,p}{t}_{m}\right)}}}+{G}_{i,cros}\left({t}_{m},{f}_{{\tau}_{m}^{2}}\right)+{n}_{G,i}\left({t}_{m},{f}_{{\tau}_{m}^{2}}\right)\end{array}$$

(7)

where ${f}_{{\tau}_{m}^{2}}$ is the frequency domain with respect to ${\tau}_{m}$. $\delta (\u2022)$ is the Dirac delta function. ${G}_{i,cros}\left({t}_{m},{f}_{{\tau}_{m}^{2}}\right)$ and ${n}_{G,i}\left({t}_{m},{f}_{{\tau}_{m}^{2}}\right)$ denote the cross tern and noise after the NUFFT, respectively.

In (7), the energy of the auto term peaks along the inclined line ${f}_{{\tau}_{m}^{2}}-{\varphi}_{2,p}-{\varphi}_{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}\left({t}_{m},{f}_{{\tau}_{m}^{2}}\right)$ and perform the IFFT with respect to ${f}_{{\tau}_{m}^{2}}$:

$$\begin{array}{ll}{Q}_{i}\left({t}_{m},{b}_{m}\right)& ={\mathrm{IFFT}}_{{f}_{{\tau}_{m}^{2}}}\left[\left|{G}_{i}\left({t}_{m},{f}_{{\tau}_{m}^{2}}\right)\right|\right]\\ & \approx \underset{\text{the auto term}}{\underbrace{{\displaystyle \sum _{p=1}^{\mathrm{P}}{\sigma}_{p}^{2}\mathrm{exp}\left[j2\pi \left({\varphi}_{2,p}+{\varphi}_{3,p}{t}_{m}\right){b}_{m}\right]}}}+{Q}_{i,cros}\left({t}_{m},{b}_{m}\right)+{n}_{Q,i}\left({t}_{m},{b}_{m}\right)\end{array}$$

(8)

where ${b}_{m}$ is the time domain with respect to ${f}_{{\tau}_{m}^{2}}$.$|\u2022|$ denotes the operation of taking the complex modulus. ${Q}_{i,cros}\left({t}_{m},{b}_{m}\right)$ and ${n}_{Q,i}\left({t}_{m},{b}_{m}\right)$ 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 $\mathrm{exp}\left[j2\pi \left(2{\varphi}_{1,p}{t}_{m}+{\varphi}_{2,p}{t}_{m}^{2}+\left({\varphi}_{3,p}/3\right){t}_{m}^{3}\right)\right]$ and the IFFT transforms $\left|{G}_{i}\left({t}_{m},{f}_{{\tau}_{m}^{2}}\right)\right|$ 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}\left({t}_{m},{b}_{m}\right)$. According to the reference [22], the GDT takes the form:

$$\begin{array}{cc}\hfill D\left({f}_{\left[g\left({\Delta}_{m}\right){\mathsf{\Upsilon}}_{m}^{b}\right]},{\Delta}_{m}\right)& ={\displaystyle \underset{{\mathsf{\Upsilon}}_{m}}{\int}\left\{u\left({\Delta}_{m}\right)\mathrm{exp}\left[j2\pi \varphi g\left({\Delta}_{m}\right){\mathsf{\Upsilon}}_{m}^{b}\right]\right\}}\times \mathrm{exp}\left[-j2\pi \xi g\left({\Delta}_{m}\right){\mathsf{\Upsilon}}_{m}^{b}{f}_{\left[g\left({\Delta}_{m}\right){\mathsf{\Upsilon}}_{m}^{b}\right]}\right]d\left[{\mathsf{\Upsilon}}_{m}\right]\hfill \\ & =u\left({\Delta}_{m}\right)\delta \left({f}_{\left[g\left({\Delta}_{m}\right){\mathsf{\Upsilon}}_{m}^{b}\right]}-\frac{\varphi}{\xi}\right)\hfill \end{array}$$

(9)

where $g\left({\Delta}_{m}\right)$ and $u\left({\Delta}_{m}\right)$ are both functions of the variable ${\Delta}_{m}$. ${f}_{\left[g\left({\Delta}_{m}\right){\mathsf{\Upsilon}}_{m}^{b}\right]}$ is the scaled frequency domain with respect to ${\mathsf{\Upsilon}}_{m}$. $b$ is a known constant. $\xi $ denotes the zoom factor. $\varphi $ is an unknown parameter.

Based on the coupling in (8) and the GDT in (9), several substitutions are done as follows:

$$\{\begin{array}{c}g\left({\Delta}_{m}\right)={b}_{m},{\mathsf{\Upsilon}}_{m}={t}_{m},b=1,\xi =1\\ u\left({\Delta}_{m}\right)\mathrm{exp}\left[j2\pi \varphi g\left({\Delta}_{m}\right){\mathsf{\Upsilon}}_{m}^{b}\right]={Q}_{i}\left({t}_{m},{b}_{m}\right)\end{array}$$

(10)

With substitutions in (10), we apply the GDT to ${Q}_{i}\left({t}_{m},{b}_{m}\right)$:

$$\begin{array}{ll}{\Gamma}_{i}\left({f}_{\left[{b}_{m}{t}_{m}\right]},{b}_{m}\right)& ={\displaystyle \underset{{t}_{m}}{\int}{Q}_{i}\left({t}_{m},{b}_{m}\right)}\mathrm{exp}\left[-j2\pi {f}_{\left[{b}_{m}{t}_{m}\right]}{b}_{m}{t}_{m}\right]d\left({t}_{m}\right)\\ & =\underset{\text{the auto tern}}{\underbrace{{\displaystyle \sum _{p=1}^{\mathrm{P}}{\sigma}_{p}^{2}\mathrm{exp}\left(j2\pi {\varphi}_{2,p}{b}_{m}\right)\delta \left({f}_{\left[{b}_{m}{t}_{m}\right]}-{\varphi}_{3,p}\right)}}}+{\Gamma}_{i,cros}\left({f}_{\left[{b}_{m}{t}_{m}\right]},{b}_{m}\right)+{n}_{\Gamma ,i}\left({f}_{\left[{b}_{m}{t}_{m}\right]},{b}_{m}\right)\end{array}$$

(11)

where ${f}_{\left[{b}_{m}{t}_{m}\right]}$ is the scaled frequency domain with respect to ${t}_{m}$. ${\Gamma}_{i,cros}\left({f}_{\left[{b}_{m}{t}_{m}\right]},{b}_{m}\right)$ and ${n}_{\Gamma ,i}\left({f}_{\left[{b}_{m}{t}_{m}\right]},{b}_{m}\right)$ denote the cross term and noise after the GDT, respectively.

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}_{\left[{b}_{m}{t}_{m}\right]}-{\varphi}_{3,p}=0$. Now, we perform the FFT along the ${b}_{m}$ axis and obtain the ICPBAF:

$$\begin{array}{ll}{\Psi}_{i}\left({f}_{\left[{b}_{m}{t}_{m}\right]},{f}_{{b}_{m}}\right)& ={\mathrm{FFT}}_{{b}_{m}}\left[{\Gamma}_{i}\left({f}_{\left[{b}_{m}{t}_{m}\right]},{b}_{m}\right)\right]\\ & =\underset{\text{the auto tern}}{\underbrace{{\displaystyle \sum _{p=1}^{\mathrm{P}}{\sigma}_{p}^{2}\delta \left({f}_{{b}_{m}}-{\varphi}_{2,p}\right)\delta \left({f}_{\left[{b}_{m}{t}_{m}\right]}-{\varphi}_{3,p}\right)}}}+{\Psi}_{i,cros}\left({f}_{\left[{b}_{m}{t}_{m}\right]},{f}_{{b}_{m}}\right)+{n}_{\Psi ,i}\left({f}_{\left[{b}_{m}{t}_{m}\right]},{f}_{{b}_{m}}\right)\end{array}$$

(12)

where ${f}_{{b}_{m}}$ is the frequency domain with respect to ${b}_{m}$. ${\Psi}_{i,cros}\left({f}_{\left[{b}_{m}{t}_{m}\right]},{f}_{{b}_{m}}\right)$ and ${n}_{\Psi ,i}\left({f}_{\left[{b}_{m}{t}_{m}\right]},{f}_{{b}_{m}}\right)$ denote the cross term and noise after the FFT, respectively.

For the ICPBAF, the auto term peaks at (${\varphi}_{2,p}$,${\varphi}_{3,p}$) on the chirp rate-quadratic chirp rate domain. Parameters ${\varphi}_{2,p}$ and ${\varphi}_{3,p}$ can be estimated by constructing a cost function to (12) [23]. With these two estimated parameters, other parameters (${\sigma}_{p}$ and ${\varphi}_{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:

$${\Psi}_{i}\left({f}_{\left[{b}_{m}{t}_{m}\right]},{f}_{{b}_{m}}\right)={\mathrm{FFT}}_{{b}_{m}}\left\{{\mathrm{GDT}}_{{b}_{m}{t}_{m}}\left\{{\mathrm{IFFT}}_{{f}_{{\tau}_{m}^{2}}}\left\{\left|{\mathrm{NUFFT}}_{{\tau}_{m}^{2}}\left[\mathrm{CPBAF}\left[{s}_{i}\left({t}_{m}\right)\right]\right]\right|\right\}\right\}\right\}$$

(13)

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.

In order to formulate the cross term problem arising from multicomponent CPSs, we consider the noise-free signal with two components, $l\in \left[1,\mathrm{P}-1\right]$ and $q\in \left[l+1,\mathrm{P}\right]$, which is denoted as:

$${s}_{l,q}\left({t}_{m}\right)={\sigma}_{l}\mathrm{exp}\left[j2\pi \left({\varphi}_{1,l}{t}_{m}+\frac{1}{2}{\varphi}_{2,l}{t}_{m}^{2}+\frac{1}{6}{\varphi}_{3,l}{t}_{m}^{3}\right)\right]+{\sigma}_{q}\mathrm{exp}\left[j2\pi \left({\varphi}_{1,q}{t}_{m}+\frac{1}{2}{\varphi}_{2,q}{t}_{m}^{2}+\frac{1}{6}{\varphi}_{3,q}{t}_{m}^{3}\right)\right]$$

(14)

The CPBAF of ${s}_{l,q}\left({t}_{m}\right)$ can be presented as:

$${R}_{l,q}\left({t}_{m},{\tau}_{m}\right)={R}_{l,q,auto}\left({t}_{m},{\tau}_{m}\right)+{R}_{l,q,cros}\left({t}_{m},{\tau}_{m}\right)$$

(15)

where:

$$\begin{array}{cc}\hfill {R}_{l,q,auto}\left({t}_{m},{\tau}_{m}\right)& ={\sigma}_{l}^{2}\mathrm{exp}\left[j2\pi \left(2{\varphi}_{1,l}{t}_{m}+{\varphi}_{2,l}{t}_{m}^{2}+\frac{1}{3}{\varphi}_{3,l}{t}_{m}^{3}\right)\right]\mathrm{exp}\left[j2\pi \left({\varphi}_{2,l}+{\varphi}_{3,l}{t}_{m}\right){\tau}_{m}^{2}\right]\hfill \\ & +{\sigma}_{q}^{2}\mathrm{exp}\left[j2\pi \left(2{\varphi}_{1,q}{t}_{m}+{\varphi}_{2,q}{t}_{m}^{2}+\frac{1}{3}{\varphi}_{3,q}{t}_{m}^{3}\right)\right]\mathrm{exp}\left[j2\pi \left({\varphi}_{2,q}+{\varphi}_{3,q}{t}_{m}\right){\tau}_{m}^{2}\right]\hfill \end{array}$$

(15a)

$$\begin{array}{cc}\hfill {R}_{l,q,cros}\left({t}_{m},{\tau}_{m}\right)& =2{\sigma}_{l}{\sigma}_{q}\mathrm{exp}\left\{j2\pi \left[\left({\varphi}_{1,l}+{\varphi}_{1,q}\right){t}_{m}+\frac{1}{2}\left({\varphi}_{2,l}+{\varphi}_{2,q}\right){t}_{m}^{2}+\frac{1}{6}\left({\varphi}_{3,l}+{\varphi}_{3,q}\right){t}_{m}^{3}\right]\right\}\mathrm{cos}\{2\pi [({\varphi}_{1,l}-{\varphi}_{1,q}+\left({\varphi}_{2,l}-{\varphi}_{2,q}\right){t}_{m}\hfill \\ & +\frac{1}{2}\left({\varphi}_{3,l}-{\varphi}_{3,q}\right){t}_{m}^{2}){\tau}_{m}+\frac{1}{6}\left({\varphi}_{3,l}-{\varphi}_{3,q}\right){\tau}_{m}^{3}]\}\mathrm{exp}\left\{j2\pi \left[\frac{1}{2}\left({\varphi}_{2,l}+{\varphi}_{2,q}\right)+\frac{1}{2}\left({\varphi}_{3,l}+{\varphi}_{3,q}\right){t}_{m}\right]{\tau}_{m}^{2}\right\}\hfill \end{array}$$

(15b)

Performing the NUFFT, operation of taking the complex modulus, IFFT, GDT and FFT on ${R}_{l,q}\left({t}_{m},{\tau}_{m}\right)$, we have:

$${\Psi}_{l,q}\left({f}_{\left[{b}_{m}{t}_{m}\right]},{f}_{{b}_{m}}\right)={\mathrm{FFT}}_{{b}_{m}}\left\{{\mathrm{GDT}}_{{b}_{m}{t}_{m}}\left\{{\mathrm{IFFT}}_{{f}_{{\tau}_{m}^{2}}}\left\{\left|{\mathrm{NUFFT}}_{{\tau}_{m}^{2}}\left[{R}_{l,q,auto}\left({t}_{m},{\tau}_{m}\right)+{R}_{l,q,cros}\left({t}_{m},{\tau}_{m}\right)\right]\right|\right\}\right\}\right\}$$

(16)

With careful analyses of Equations (15a) and (15b), we find that, compared to the auto term ${R}_{l,q,auto}\left({t}_{m},{\tau}_{m}\right)$, the cosine function of ${R}_{l,q,cros}\left({t}_{m},{\tau}_{m}\right)$ 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 ${\tau}_{m}$, the influence of the cosine function will be eliminated:

$$\left({\varphi}_{1,l}-{\varphi}_{1,q}+\left({\varphi}_{2,l}-{\varphi}_{2,q}\right){t}_{m}+\frac{1}{2}\left({\varphi}_{3,l}-{\varphi}_{3,q}\right){t}_{m}^{2}\right){\tau}_{m}+\frac{1}{6}\left({\varphi}_{3,l}-{\varphi}_{3,q}\right){\tau}_{m}^{3}=n\pi ,n=\cdots -2,-1,0,1,2,\cdots $$

(17)

Obviously, only when ${\varphi}_{1,l}={\varphi}_{1,q}$, ${\varphi}_{2,l}={\varphi}_{2,q}$ and ${\varphi}_{3,l}={\varphi}_{3,q}$, Equation (17) equals to zero for every ${t}_{m}$ and ${\tau}_{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:

$${\Psi}_{i}\left({f}_{\left[{b}_{m}{t}_{m}\right]},{f}_{{b}_{m}}\right)=\underset{\text{the auto tern}}{\underbrace{{\displaystyle \sum _{p=1}^{\mathrm{P}}{\sigma}_{p}^{2}\delta \left({f}_{{b}_{m}}-{\varphi}_{2,p}\right)\delta \left({f}_{\left[{b}_{m}{t}_{m}\right]}-{\varphi}_{3,p}\right)}}}+{n}_{\Psi ,i}\left({f}_{\left[{b}_{m}{t}_{m}\right]},{f}_{{b}_{m}}\right)$$

(18)

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.

*We consider three noise-free CPSs denoted by Au1, Au2 and Au3. The sampling frequency*
${F}_{{t}_{m}}$
*[same as the pulse repetition frequency (PRF)] is 256 Hz, and the signal length*
${N}_{{t}_{m}}$
*(same as the echo pulses) is equal to 512. The signal parameters are set as follows:*
${\varphi}_{1,1}$
*= 100 Hz,*
${\varphi}_{2,1}$
*= 84 Hz/s,*
${\varphi}_{3,1}$
*= 80 Hz/s ^{2} for Au1;*
${\varphi}_{1,2}$

By performing the NUFFT along the ${\tau}_{m}$ axis of the CPBAF in (6), we obtain the slow time-chirp rate distribution in Figure 2a. The bilinearity of the CPBAF causes the situation that the auto term and cross term coexist. In Figure 2a, we can find that, due to the coupling between ${t}_{m}$ and ${\tau}_{m}$, the auto term takes the form of three inclined lines. In order to accumulate the auto term, we take the complex modulus, and then, perform the IFFT, GDT and FFT sequentially. Figure 2b shows the contour of the ICPBAF and its stereogram is shown in Figure 2c. Obviously, as analyzed above, the auto term accumulates into the ideal spread function, while the cross term cannot accumulate and even can be ignored. In Figure 2c, with the peak detection technique [22,23], (${\varphi}_{2,1}$, ${\varphi}_{3,1}$), (${\varphi}_{2,2}$, ${\varphi}_{3,2}$) and (${\varphi}_{2,3}$, ${\varphi}_{3,3}$) are estimated as (84 Hz/s, 80 Hz/s^{2}), (12 Hz/s, 10 Hz/s^{2}) and (−64 Hz/s, −50 Hz/s^{2}), respectively. Thereafter, compensating the phase term pertaining to the estimated parameters and performing the FFT, we estimate (${\sigma}_{1}$, ${\varphi}_{1,1}$), (${\sigma}_{2}$, ${\varphi}_{1,2}$) and (${\sigma}_{3}$, ${\varphi}_{1,3}$) as (1, 100 Hz), (1, 20 Hz) and (1, −80 Hz), respectively.

Figure 2a–c consider multicomponent CPSs with the same amplitude. However, the amplitudes are always different in realistic applications. Figure 2d considers multicomponent CPSs with different amplitudes. With the result shown in Figure 2d, we know that, when the differences between amplitudes are large, the auto terms of weak CPSs may be submerged in the residual cross terms generated by the strong CPSs and the Clean technique [10,12,13] should be employed to separate the strong and weak CPSs.

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.

Assume the length of slow time samples ${N}_{{t}_{m}}$ is equal to the length of lag samples ${N}_{{\tau}_{m}}$. The ICPBAF implementation needs the CPBAF [$O\left({N}_{{t}_{m}}^{2}\right)$], NUFFT along ${\tau}_{m}$[$O\left({N}_{{t}_{m}}^{2}{\mathrm{log}}_{2}{N}_{{t}_{m}}\right)$], operation of taking the complex modulus [$O\left({N}_{{t}_{m}}^{2}\right)$], IFFT along ${f}_{{\tau}_{m}^{2}}$[$O\left({N}_{{t}_{m}}^{2}{\mathrm{log}}_{2}{N}_{{t}_{m}}\right)$], GDT [$O\left({N}_{{t}_{m}}^{2}{\mathrm{log}}_{2}{N}_{{t}_{m}}\right)$] and FFT along ${b}_{m}$[$O\left({N}_{{t}_{m}}^{2}{\mathrm{log}}_{2}{N}_{{t}_{m}}\right)$]. Therefore, the overall computational cost of the ICPBAF is in the order of $O\left({N}_{{t}_{m}}^{2}{\mathrm{log}}_{2}{N}_{{t}_{m}}\right)$ and listed in 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 ${\tau}_{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 ${\tau}_{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.

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.

*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:*
${\sigma}_{1}=1$*,*
${\varphi}_{1,1}$
*= 10 Hz,*
${\varphi}_{2,1}$
*= 2 Hz/s,*
${\varphi}_{3,1}$
*= 2 Hz/s ^{2} for Bu1;*
${\sigma}_{2}=1$

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.

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) ($SN{R}_{out}$ is listed in (19)) [31] and mean square error (MSE) [32,33] are utilized as measures of the noise resistance:

$${\mathrm{SNR}}_{out}=10{\mathrm{log}}_{10}\frac{{\sigma}_{p}^{2}}{{N}_{{t}_{m}}{\sigma}^{2}}{\left\{{\left|{\displaystyle \sum _{m=-\frac{{N}_{{t}_{m}}}{2}}^{\frac{{N}_{{t}_{m}}}{2}-1}{s}_{i}\left(m\right)}\mathrm{exp}\left[-j\pi {a}_{2.p}^{\prime}{\left(\frac{m}{{F}_{{t}_{m}}}\right)}^{2}-j\pi \frac{{a}_{3,p}^{\prime}}{3}{\left(\frac{m}{{F}_{{t}_{m}}}\right)}^{3}\right]\right|}_{\mathrm{max}}\right\}}^{2}$$

(19)

where ${\sigma}^{2}$ is the power of the complex white Gaussian noise. ${a}_{2.p}^{\prime}$ and ${a}_{3,p}^{\prime}$ are estimations.

*We consider a CPS denoted by Cu.*
${F}_{{t}_{m}}$
*and*
${N}_{{t}_{m}}$
*are 256 Hz and 256, respectively. The signal parameters are set as follows:*
${\varphi}_{1,1}$
*= 106 Hz,*
${\varphi}_{2,1}$
*= 100 Hz/s,*
${\varphi}_{3,1}$
*= 80 Hz/s ^{2} for Cu. The tested input SNRs are SNR_{in} = [−11:1:0] and 200 trials have been performed for each SNR_{in} under the ten-time interpolation. Figure 4 and Figure 5 gives the input-output SNR and MSE, respectively.*

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.

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\left(\widehat{t}\right)=\mathrm{rect}\left(\widehat{t}/{T}_{s}\right)\mathrm{exp}\left(j\pi \gamma {\widehat{t}}^{2}\right)$ (where $\widehat{t}$, ${T}_{s}$ and $\gamma $ 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}\left({t}_{m}\right)$ in the*i*-th range cell.*Step 4*: Calculate the energy of the extracted data ${s}_{i}\left({t}_{m}\right)$. If the energy is smaller than the set threshold ${E}_{s}$ [8,12], set $i=i+1$ and repeat Step 3 until $i=\mathrm{I}$.*Step 5*: Apply the ICPBAF to ${s}_{i}\left({t}_{m}\right)$ and estimate ${\varphi}_{2,p}$ and ${\varphi}_{3,p}$ with the peak detection technique [23,24].*Step 6*: Dechirp ${s}_{i}\left({t}_{m}\right)$ with $\mathrm{exp}\left[-j2\pi \left(\frac{{\varphi}_{2,p}^{\prime}}{2}{t}_{m}^{2}+\frac{{\varphi}_{3,p}^{\prime}}{6}{t}_{m}^{3}\right)\right]$, and then estimate ${\varphi}_{1,p}$ and ${\sigma}_{p}$ via the FFT and peak detection technique.where ${\sigma}_{p}^{\prime}$ and ${\varphi}_{1,p}^{\prime}$ denote estimations of the centroid frequency and amplitude for the$$\left({\sigma}_{p}^{\prime}=\frac{{A}^{\prime}}{{N}_{{t}_{m}}},{\varphi}_{1,p}^{\prime}={f}_{{t}_{m}}^{\prime}\right)=\begin{array}{c}\mathrm{arg}\mathrm{max}\\ \left(A,{f}_{{t}_{m}}\right)\end{array}\left|{\mathrm{FFT}}_{{t}_{m}}\left({s}_{i}\left({t}_{m}\right)\mathrm{exp}\left[j2\pi \left(-\frac{{\varphi}_{2,p}^{\prime}}{2}{t}_{m}^{2}-\frac{{\varphi}_{3,p}^{\prime}}{6}{t}_{m}^{3}\right)\right]\right)\right|$$(20)*p*th CPS, respectively. ${A}^{\prime}$ denotes the peak value after the FFT.*Step 7*: Eliminate the estimated*p*-th CPS from the original signal ${s}_{i}\left({t}_{m}\right)$where ${\mathrm{win}}_{p}\left({f}_{{t}_{m}}\right)=\{\begin{array}{c}0,{\varphi}_{1,p}^{\prime}-d\le {f}_{{t}_{m}}\le {\varphi}_{1,p}^{\prime}+d\\ 1,else\end{array}$ denotes the narrowband filter with the$${s}_{i}\left({t}_{m}\right)={\mathrm{FFT}}_{{t}_{m}}\left({s}_{i}\left({t}_{m}\right)\mathrm{exp}\left[j2\pi \left(-\frac{{\varphi}_{2,p}^{\prime}}{2}{t}_{m}^{2}-\frac{{\varphi}_{3,p}^{\prime}}{6}{t}_{m}^{3}\right)\right]\right){\mathrm{win}}_{p}\left({f}_{{t}_{m}}\right)\mathrm{exp}\left[j2\pi \left(\frac{{\varphi}_{2,p}^{\prime}}{2}{t}_{m}^{2}+\frac{{\varphi}_{3,p}^{\prime}}{6}{t}_{m}^{3}\right)\right]$$(21)*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<\mathrm{I}$, set*i*=*i*+1 and repeat Steps 3–8 until*i*= I.

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.

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.

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.

Ship target model and results of the standard RD technique. (**a**) Ship target model; (**b**) Result of the standard RD technique under the absence of the scatterer-dependent Doppler spread; (**c**) Result of the standard RD technique under the existence of the **...**

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\left(i,n\right)$ in Table 3 [12]:

$$ENT=-{\displaystyle \sum _{i=1}^{\mathrm{I}}{\displaystyle \sum _{n=1}^{{N}_{{t}_{m}}}{\left|X\left(i,n\right)\right|}^{2}\mathrm{ln}{\left|X\left(i,n\right)\right|}^{2}}}$$

(22)

ISAR imaging results of the synthetic model. (**a**) RID ISAR imaging algorithm using the CPF; (**b**) RID ISAR imaging algorithm using the IGCPF; (**c**) RID ISAR imaging algorithm using the SCFT-based algorithm; (**d**) RID ISAR imaging algorithm using the GSCFT-based **...**

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.

Ratio between the number of the correctly reconstructed scatterers and the number of all reconstructed scatterers.

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.

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.

Processing results of the real radar data. (**a**) Azimuth echoes after the motion compensation; (**b**) Wigner-Ville distribution of the 191th 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.

Processing results of the extracted data in the 191th range cell. (**a**) ICPBAF of the extracted data; (**b**) Results of the standard RD technique and ICPBAF-based RID ISAR imaging algorithm.

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.

ISAR imaging results of the real radar data. (**a**) Standard RD technique; (**b**) RID ISAR imaging algorithm using the CPF; (**c**) RID ISAR imaging algorithm using the IGCPF; (**d**) RID ISAR imaging algorithm using the SCFT-based algorithm; (**e**) RID ISAR imaging algorithm **...**

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 Section 4 and Section 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.

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.

This work was supported in part by the National Natural Science Foundation of China (grant 61601341), in part by the Project Funded by China Postdoctoral Science Foundation (grant 2015M582615 and 2016T90891) and in part by the Program for the National Science Fund for Distinguished Young Scholars (grant 61525105).

Author Contributions

Jibin Zheng proposed the ISAR imaging algorithm and wrote the paper; Hongwei Liu, Zheng Liu and Qing Huo Liu helped revise the paper.

Conflicts of Interest

The authors declare no conflict of interest.

1. Prickett M.J., Chen C.C. Principles of inverse synthetic aperture radar/ ISAR/imaging; Proceedings of the Electronics and Aerospace Systems Conference; Arlington, VA, USA. 29 September–1 October 1980.

2. Soumekh M. Synthetic Aperture Radar Signal Processing. Wiley; New York, NY, USA: 1999.

3. Berizzi F., Mese E., Diania M., Martorella M. High-resolution ISAR imaging of maneuvering targets by means of the range instantaneous Doppler technique: Modeling and performance analysis. IEEE Trans. Image Process. 2001;10:1880–1890. doi: 10.1109/83.974573. [PubMed] [Cross Ref]

4. Wang Y., Abdelkader A., Zhao B., Wang J. ISAR imaging of maneuvering targets based on the modified discrete polynomial-phase transform. Sensors. 2015;15:22401–22418. doi: 10.3390/s150922401. [PMC free article] [PubMed] [Cross Ref]

5. Wang J., Kasilingam D. Global range alignment for ISAR. IEEE Trans. Aerosp. Electron. Syst. 2003;39:351–357. doi: 10.1109/TAES.2003.1188917. [Cross Ref]

6. Xing M., Wu R., Lan J., Bao Z. Migration through resolution cell compensation in ISAR imaging. IEEE Geosci. Remote Sens. Lett. 2004;1:141–144. doi: 10.1109/LGRS.2004.824766. [Cross Ref]

7. Wahl D.E., Eichel P.H., Ghigetia D.C. Phase gradient autofocus-a robust tool for high resolution SAR phase correction. IEEE Trans. Aerosp. Electron. Syst. 1994;30:827–835. doi: 10.1109/7.303752. [Cross Ref]

8. Xing M., Wu R., Li Y., Bao Z. New ISAR imaging algorithm based on modified Wigner Ville distribution. IET Radar Sonar Navig. 2008;3:70–80. doi: 10.1049/iet-rsn:20080003. [Cross Ref]

9. Lv X., Xing M., Wang C., Zhang S. ISAR imaging of maneuvering targets based on the range centroid Doppler technique. IEEE Trans. Image Process. 2010;19:141–153. [PubMed]

10. Lv X., Bi G., Wang C., Xing M. Lv’s distribution: Principle, implementation, properties, and performance. IEEE Trans. Signal Process. 2011;59:3576–3591. doi: 10.1109/TSP.2011.2155651. [Cross Ref]

11. Guo X., Sun H., Wang S., Liu G. Comments on Discrete Chirp-Fourier Transform and its Application to Chirp Rate Estimation. IEEE Trans. Signal Process. 2002;50:3115.

12. Zheng J., Su T., Zhu W., Zhang L., Liu Z., Liu Q. ISAR imaging of nonuniformly rotating target based on a fast parameter estimation algorithm of cubic phase signal. IEEE Trans. Geosci. Remote Sens. 2015;53:4727–4740. doi: 10.1109/TGRS.2015.2408350. [Cross Ref]

13. Bai X., Tao R., Wang Z., Wang Y. ISAR imaging of a ship target based on parameter estimation of multicomponent quadratic frequency-modulated signals. IEEE Trans. Geosci. Remote Sens. 2014;52:1418–1429. doi: 10.1109/TGRS.2013.2251348. [Cross Ref]

14. Wu L., Wei X., Yang D., Wang H., Li X. ISAR imaging of targets with complex motion based on discrete chirp Fourier transform for cubic chirps. IEEE Trans. Geosci. Remote Sens. 2012;50:4201–4212. doi: 10.1109/TGRS.2012.2189220. [Cross Ref]

15. Nguyen N., Liu Q.H. The regular Fourier matrices and nonuniform fast Fourier transforms. SIAM J. Sci. Comput. 1999;21:283–293. doi: 10.1137/S1064827597325712. [Cross Ref]

16. Liu Q.H., Nguyen N., Tang X.Y. Accurate algorithms for nonuniform fast forward and inverse Fourier transforms and their applications. Proc. IEEE Geosci. Remote Sens. Symp. 1998;1:288–290.

17. Barbarossa S., Scaglione A., Giannakis G.B. Product high-order ambiguity function for multicomponent polynomial-phase signal modeling. IEEE Trans. Signal Process. 1998;46:691–708. doi: 10.1109/78.661336. [Cross Ref]

18. Wang Y., Jiang Y. Inverse synthetic aperture radar imaging of maneuvering target based on the product generalized cubic phase function. IEEE Geosci. Remote Sens. Lett. 2011;8:958–962. doi: 10.1109/LGRS.2011.2143387. [Cross Ref]

19. Wang Y., Lin Y. ISAR imaging of non-uniformly rotating target via range-instantaneous-Doppler-derivatives algorithm. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2014;7:167–176. doi: 10.1109/JSTARS.2013.2257699. [Cross Ref]

20. Zheng J., Liu H., Liao G., Su T., Liu Z., Liu Q.H. ISAR imaging of targets with complex motions based on a noise-resistant parameter estimation algorithm without nonuniform axis. IEEE Sens. J. 2016;16:2509–2518. doi: 10.1109/JSEN.2016.2516040. [Cross Ref]

21. Li D., Gui X., Liu H., Su J., Xiong H. An ISAR imaging algorithm for maneuvering targets with low SNR based on parameter estimation of multicomponent quadratic FM signals and nonuniform FFT. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2016;9:5688–5702. doi: 10.1109/JSTARS.2016.2543233. [Cross Ref]

22. Zheng J., Liu H., Liao G., Su T., Liu Z., Liu Q.H. ISAR Imaging of Nonuniformly Rotating Targets Based on Generalized Decoupling Technique. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2016;9:520–532. doi: 10.1109/JSTARS.2015.2509169. [Cross Ref]

23. O’Shea P. A new technique for instantaneous frequency rate estimation. IEEE Signal Process. Lett. 2002;9:251–252. doi: 10.1109/LSP.2002.803003. [Cross Ref]

24. O’Shea P. A fast algorithm for estimating the parameters of a quadratic FM signal. IEEE Trans. Signal Process. 2004;52:385–393. doi: 10.1109/TSP.2003.821097. [Cross Ref]

25. Simeunović M., Djurović I. Non-uniform sampled cubic phase function. Signal Process. 2014;101:99–103. doi: 10.1016/j.sigpro.2014.02.005. [Cross Ref]

26. Wang Y., Kang J., Jiang Y. ISAR imaging of maneuvering target based on the local polynomial Wigner distribution and integrated high-order ambiguity function for cubic phase signal model. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2014;7:2971–2991. doi: 10.1109/JSTARS.2014.2301158. [Cross Ref]

27. Farquharson M., O’Shea P. Extending the performance of the cubic phase function algorithm. IEEE Trans. Signal Process. 2007;55:4767–4774. doi: 10.1109/TSP.2007.896085. [Cross Ref]

28. Zheng J., Su T., Liao G., Liu H., Liu Z., Liu Q.H. ISAR imaging for fluctuating ships based on a fast bilinear parameter estimation algorithm. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2015;8:3954–3966. doi: 10.1109/JSTARS.2015.2440911. [Cross Ref]

29. Chen X., Guan J., Liu N., He Y. Maneuvering target detection via radon-Fractional Fourier transform-based long-time coherent integration. IEEE Trans. Signal Process. 2014;62:939–953. doi: 10.1109/TSP.2013.2297682. [Cross Ref]

30. Wang P., Li H., Djurović I., Himed B. Integrated cubic phase function for linear FM signal analysis. IEEE Trans. Aerosp. Electron. Syst. 2010;46:963–977. doi: 10.1109/TAES.2010.5545167. [Cross Ref]

31. Zheng J., Su T., Liu Q.H., Zhang L., Zhu W. Fast parameter estimation algorithm for cubic phase signal based on quantifying effects of Doppler frequency shift. Prog. Electromagn. Res. 2013;142:57–74. doi: 10.2528/PIER13061008. [Cross Ref]

32. Peleg S., Porat B. The Cramer-Rao lower bound for signals with constant amplitude and polynomial phase. IEEE Trans. Aerosp. Electron. Syst. 1998;34:1070–1084. doi: 10.1109/78.80864. [Cross Ref]

33. Peleg S., Porat B. Linear FM signal parameter estimation from discrete-time observations. IEEE Trans. Aerosp. Electron. Syst. 1991;27:607–616. doi: 10.1109/7.85033. [Cross Ref]

Articles from Sensors (Basel, Switzerland) are provided here courtesy of **Multidisciplinary Digital Publishing Institute (MDPI)**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |