A Novel Bistatic SAR Maritime Ship Target Imaging Algorithm Based on Cubic Phase Time-Scaled Transformation

Abstract

Due to the advantages of flexible configuration, bistatic synthetic aperture radar (BiSAR) has the ability to effectively observe from various visual angles, such as forward view area and squint area, and has good anti-jamming characteristics. It can be applied to the surveillance of ship targets on the sea and is gradually gaining an increasing amount of attention. However, for ship targets with complex motions on the sea surface, such as maneuvering targets or ship targets under high sea conditions, the high-order Doppler frequency of the scattering points is always spatial variation (related to the spatial position of scattering points), which poses a considerable challenge for the imaging of maritime ship targets in BiSAR. To resolve this problem, a BiSAR maritime ship target imaging algorithm based on cubic phase time-scaled transformation is proposed in this paper. First, through pre-processing of echo such as Doppler prefiltering and keystone transform, the translation compensation of the BiSAR maritime ship target is completed, and the scattering point energy is corrected to within one range unit. Then, the azimuth signal is modeled as a multi-component cubic phase signal. Based on the proposed cubic phase time-scaled transformation, the Doppler centroid, frequency rate, and third-order frequency of scattering points are estimated. Eventually, the BiSAR imaging of maritime ship targets is realized. This algorithm has excellent noise immunity and low cross-terms. The simulation leads to the verification of the validity of the proposed algorithm.

1. Introduction

Synthetic aperture radar (SAR) has been used widely in civil and military applications due to its all weather and all time advantage [1,2,3]. The bistatic SAR (in which the transmitter and the receiver are mounted on different platforms), thanks to its flexible configuration, has practical observation in the front view area, squint area, and other visual angles. This makes up for its shortcomings, as the mono-static SAR cannot realize forward-looking imaging. In addition, the bistatic SAR has a good anti-jamming ability due to the muted reception of the receiver. In the past years, the bistatic SAR has gradually received more and more attention [4,5].

For stationary targets and scenes, bistatic SAR can obtain their two-dimensional high-resolution images, in which the azimuthal resolution is obtained by the perspective variation caused by the movements of the transmitter and receiver, and range resolution is obtained through the transmitted large bandwidth signal. In the existing literature, $\omega -$ k [6,7], nonlinear chirp scaling (NLCS) [8,9] and other algorithms [10,11,12,13,14] can achieve high-resolution imaging of stationary targets and scenes. In addition, high-resolution imaging of ground/sea moving targets has become a long-standing hot topic in order to meet the growing need for monitoring and surveillance of moving targets. Therefore, by combining the advantages of bistatic SAR, it will be of great significance to apply it to moving target surveillance.

However, for moving targets, the high azimuthal resolution is obtained from the perspective variation between the bistatic SAR platform’s movements and the target movements. Therefore, in order to obtain high-resolution images of moving targets, it is necessary to accurately estimate the target motion parameters from the echoes in bistatic SAR [15,16]. In fact, many bistatic SAR experiments have been carried out. The University of Electronic Science and Technology (UESTC) carried out bistatic side- and forward-looking SAR experiments in 2007 and 2012, respectively [17,18], and carried out the imaging experiment of bistatic forward-looking SAR ground moving targets in 2020 [19], both of which achieved good imaging results. For bistatic SAR ground moving targets, an effective and simple technology is to estimate the Doppler frequency modulation rate and other Doppler parameters of the target after the range cell migration (RCM) correction of the target is achieved, and to adopt the azimuth matching filtering method, which can achieve good focusing of ground moving targets, so as to obtain the high-resolution image of the target [20,21,22], such as mismatched compression [20], and spectrum matching method [21]. However, a maritime ship target has time-varying and unknown rotation under the disturbance of the sea wave, resulting in the Doppler frequency of the target scattering points always being time-varying, and the high-order Doppler frequency always being related to the position of the scattering points, that is, its spatial variation. Therefore, the above imaging methods for ground moving targets will cause serious defocusing of maritime ship target images.

For the time-varying and spatial variation characteristics of target Doppler frequency, researchers have proposed the range-instantaneous-Doppler (RID) technology [23,24,25,26,27,28,29,30,31,32,33,34]. In RID technology, after translation compensation, the azimuth signal of scattering points in a range gate can be modeled as a linear frequency modulation (LFM) signal or a cubic phase signal (CPS). Specifically, the LFM signal model can be used to accurately estimate the Doppler centroid and frequency rate of scattering points when the ship target is moving stably, such as cruising in low sea conditions. These methods, based on LFM [23,24,25,26,27,28], can obtain well-focused ship target images. For BiSAR ship target imaging, Li et al. [28] proposed the hybrid SAR/ISAR image formation via joint FrFT-WVD processing based on the LFM model. In this method, scattering points of ship targets were separated by fractional Fourier transform (FrFT), and the high-resolution time-frequency distribution was obtained by Wigner–Ville distribution (WVD). In practice, when the target is in complex motion, such as maneuver motion or disturbance under high sea conditions, the LFM model is no longer applicable, and a higher precision CPS model is required. The methods based on the CPS model are mainly divided into two categories: non-correlated algorithms and correlated algorithms. For uncorrelated algorithms, maximum likelihood (ML) [29] and discrete chirp Fourier transform (DCFT) [30] require the multi-dimensional violent search to obtain the target’s Doppler centroid, frequency modulation rate, third-order term coefficient, and other parameters. It has the characteristics of heavy computational burden and low processing efficiency. Therefore, to improve the calculation efficiency, these algorithms [31,32,33,34] mainly include the scaled Fourier transform (SCFT)-based algorithm [31], the generalized scaled Fourier transform (GSCFT) [33], the coherently integrated generalized CPF (CIGCPF) [34], the chirp rate–quadratic chirp rate distribution (CRQCRD) [32], etc., through the autocorrelation processing of azimuth signals. However, these methods either require high SNR or have serious cross-terms, which limits the high precision imaging of maritime ship targets.

In order to solve the above problems, a bistatic SAR maritime ship target imaging algorithm based on cubic phase time scale transformation is proposed, which has excellent noise immunity and low cross-terms. First, through the pre-processing of echoes such as Doppler prefiltering and keystone transform [35], the translation compensation of the BiSAR maritime ship target is completed, and the scattering point energy is corrected to within one range unit. Then, the azimuth signal is modeled as a multi-component cubic phase signal. Based on the proposed cubic phase time-scaled transformation, the coupling term relationship between azimuth time variables and lag time variables is changed to hyperbola. In addition, after a one-dimensional search of second-order term coefficients, two-dimensional (2-D) coherent integration can be achieved by phase compensation and 2-D fast Fourier transform (FFT), so as to estimate the Doppler centroid, Doppler frequency rate, and the third-order frequency of scattering points. After all the range gates are processed, the high-resolution image of the maritime ship target in BiSAR can be obtained.

This paper is structured as follows. Section 2 gives the imaging geometry and echo model. In Section 3, the proposed algorithm is introduced in detail, including pre-processing, cubic phase time-scaled transformation, and the BiSAR maritime ship target imaging method. Section 4 gives the performance analysis of parameter estimation based on cubic phase time scale transform, mainly from the cross-term and SNR. Section 5 gives the simulation of BiSAR maritime ship target imaging. Finally, Section 6 gives the conclusion.

2. Imaging Geometry and Data Model

Figure 1 illustrates the BiSAR geometry for the maritime ship target, which is represented in the Cartesian coordinate system O-$XYZ$, and the origin O is the projection point of the receiver on the ground plane. The receiver and transmitter travel at speeds $V_R$ and $V_T$, respectively, in the Cartesian coordinate system. In Figure 1, o-$xyz$ is the body-fixed reference frame whose origin o is located at the rotation center of the target. The x-axis and y-axis are the longitudinal axis (directed from aft to fore) and the transversal axis (directed to larboard), respectively. The z-axis is the normal axis (directed from bottom to top), which is perpendicular to the x-axis and y-axis. The body-fixed reference frame o-$xyz$ moves and oscillates with the target, which is used to describe the positions of the scatterers relative to the rotation center of the target o. The reference coordinate $O^{\prime }$-$X^{\prime }{Y}^{\prime }{Z}^{\prime }$ is a translation of O-$XYZ$ with the origin $O^{\prime }$ located at rotation center of the target. $\phi$ represents the heading angle of the maritime ship target in the reference coordinate $O^{\prime }$-$X^{\prime }{Y}^{\prime }{Z}^{\prime }$ (the included angle with the $X^{\prime }$-axis).

When the target rotates with an instantaneous angle of roll ${\theta }_x$, pitch ${\theta }_y$ and yaw ${\theta }_z$, the instantaneous total rotational matrix $\Re ({\theta }_x,{\theta }_y,{\theta }_z;\xi )$ can be derived as

(1)$\begin{array}{cc}\hfill \Re ({\theta }_x,{\theta }_y,{\theta }_z;\xi )& ={\Re }_x\left(\xi \right){\Re }_y\left(\xi \right){\Re }_z\left(\xi \right)\hfill \\ \hfill & =\left[\begin{array}{ccc}1& 0& 0\\ 0& cos{\theta }_x& -sin{\theta }_x\\ 0& sin{\theta }_x& cos{\theta }_x\end{array}\right]\left[\begin{array}{ccc}cos{\theta }_y& 0& sin{\theta }_y\\ 0& 1& 0\\ -sin{\theta }_y& 0& cos{\theta }_y\end{array}\right]\left[\begin{array}{ccc}cos{\theta }_z& -sin{\theta }_z& 0\\ sin{\theta }_z& cos{\theta }_z& 0\\ 0& 0& 1\end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{ccc}cos{\theta }_{y}cos{\theta }_z& -cos{\theta }_{y}sin{\theta }_z& sin{\theta }_y\\ sin{\theta }_{y}sin{\theta }_{x}cos{\theta }_z+cos{\theta }_{x}sin{\theta }_z& cos{\theta }_{x}cos{\theta }_z-sin{\theta }_{y}sin{\theta }_{x}sin{\theta }_z& -cos{\theta }_{y}sin{\theta }_x\\ sin{\theta }_{x}sin{\theta }_z-sin{\theta }_{y}cos{\theta }_{x}cos{\theta }_z& sin{\theta }_{y}cos{\theta }_{x}sin{\theta }_z+sin{\theta }_{x}cos{\theta }_z& cos{\theta }_{y}cos{\theta }_x\end{array}\right]\hfill \end{array}$

where ${\Re }_x\left(\xi \right)$, ${\Re }_y\left(\xi \right)$ and ${\Re }_z\left(\xi \right)$ are the instantaneous roll, pitch, and yaw rotational matrix, respectively.

According to the Rodrigues’ formula [36] and instantaneous total rotational matrix $\Re ({\theta }_x,{\theta }_y,{\theta }_z;\xi )$, the synthetic rotation angle ${\theta }_s$ at the time $\xi$ can be given by

(2)${\theta }_s=arccos\left(\frac{\mathrm{tr}(\Re ({\theta }_x,{\theta }_y,{\theta }_z;\xi ))-1}{2}\right)$

where $\mathrm{tr}(•)$ denote the trace of a matrix.

Moreover, the three-dimensional vector $\mathit{\omega }$, which represents the counterclockwise rotation of the ship target around this direction, and whose magnitude represents the rotational speed of the ship target at $\xi$, can be given by

(3)$\mathit{\omega }=\frac{1}{2sin{\theta }_s}{\left[r_{32}-r_{23}{r}_{13}-r_{31}{r}_{21}-r_{12}\right]}^T$

where, $r_{ij}$ denote the ith-row and jth-col element of the instantaneous total rotational matrix $\Re ({\theta }_x,{\theta }_y,{\theta }_z;\xi )$.

For an arbitrary scattering point P at position ${\mathit{r}}_p=(x_P,y_P,z_P)$, its Doppler frequency can be derived by

(4)$\begin{array}{c}\hfill f_d=\frac{1}{\lambda }\left[v_r+\left(\mathit{\omega }\times {\mathit{r}}_p\right)·{\mathit{i}}_R\right]+\frac{1}{\lambda }\left[v_t+\left(\mathit{\omega }\times {\mathit{r}}_p\right)·{\mathit{i}}_T\right]\end{array}$

where × and · represent respectively the outer and inner product operations. $v_r$ and $v_t$ represent the radial velocity of between the receiver or transmitter and the target, respectively. ${\mathit{i}}_R$ and ${\mathit{i}}_T$ indicate the unit vector of radar line of sight (RLOS) of the receiver and transmitter, respectively.

Due to BiSAR platforms’ movements and the disturbance of sea waves, the magnitude and direction of the target’s rotation axis vary with time, i.e., $v_r$, $v_t$ and $\mathit{\omega }$ in (4) are time-varying, so they can be approximately expressed as

(5)$\begin{array}{cc}\hfill v_r& \approx v_{r0}+v_r^{\prime }\xi +0.5v_r^{″}{\xi }^2\hfill \\ \hfill v_t& \approx v_{t0}+v_t^{\prime }\xi +0.5v_t^{″}{\xi }^2\hfill \\ \hfill \mathit{\omega }& \approx {\mathit{\omega }}_0+{\mathit{\omega }}^{\prime }\xi +0.5{\mathit{\omega }}^{″}{\xi }^2\hfill \end{array}$

where ${[\ast ]}_0$, ${[\ast ]}^{\prime }$, and ${[\ast ]}^{″}$ indicate the constant term and first and second order term coefficients of $[\ast ]$, respectively. (Note: $[\ast ]$ denotes $v_r$ or $v_t$ or $\mathit{\omega }$ in (5).)

Then, the instantaneous two-way range history $R\left(\xi \right)$ of the arbitrary scattering point P can be denoted as

(6)$\begin{array}{cc}\hfill R(\xi ;{\mathit{r}}_p)=& R\left(t_0\right)+\frac{\lambda }{2}\int f_{d}d\xi \hfill \\ \hfill =& R\left(t_0\right)+(v_{r0}+v_{t0}+{\mathit{\omega }}_0·\mathit{\mu })\xi +\frac{1}{2}(v_r^{\prime }+v_t^{\prime }+{\mathit{\omega }}^{\prime }·\mathit{\mu }){\xi }^2+\frac{1}{6}(v_r^{″}+v_t^{″}+{\mathit{\omega }}^{″}·\mathit{\mu }){\xi }^3\hfill \end{array}$

where $\mathit{\mu }={\mathit{r}}_p\times \left({\mathit{i}}_R+{\mathit{i}}_T\right)$.

Assuming that the radar transmits the LFM signal as follows:

(7)$s\left(t\right)=\mathrm{rect}\left(\frac{t}{T_p}\right)exp\left\{j2\pi \left(f_{c}t+\frac{1}{2}{K}_r{t}^2\right)\right\}$

where $\mathrm{rect}\left(t\right)=\left\{\begin{array}{c}1,\left|t\right|\le 1/2\hfill \\ 0,\left|t\right|\ge 1/2\hfill \end{array}\right.$, t is the fast time, $T_p$ denotes the width of the signal pulse. $f_c$ and $K_r$ indicate the carrier frequency and the chirp rate, respectively.

Therefore, the echo signal of any scattering point P in the maritime ship target in BiSAR is

(8)$\begin{array}{cc}\hfill s_p(t,\xi )={\sigma }_p\left(\xi \right)\mathrm{rect}& \left(\frac{t-R\left(\xi ;{\mathit{r}}_p\right)/c}{T_p}\right)\hfill \\ \hfill & \times exp\left\{j2\pi \left[f_c\left(t-\frac{R\left(\xi ;{\mathit{r}}_p\right)}{c}\right)+\frac{1}{2}{K}_r{\left(t-\frac{R\left(\xi ;{\mathit{r}}_p\right)}{c}\right)}^2\right]\right\}\hfill \end{array}$

where ${\sigma }_p\left(\xi \right)$ is the reflection coefficient of the scattering point P at the slow time $\xi$ and c is the speed of the light.

The echo signal of the arbitrary scattering point P is transformed to the baseband and range compression, which can be expressed as

(9)$\begin{array}{c}\hfill s_p(t,\xi )={\sigma }_p\left(\xi \right)\sin \mathrm{c}\left\{B_r\left[t-\frac{R\left(\xi ;{\mathit{r}}_p\right)}{c}\right]\right\}\times exp\left\{-j\frac{2\pi }{\lambda }R\left(\xi ;{\mathit{r}}_p\right)\right\}\end{array}$

where $B_r=K_r{T}_p$ and $\lambda =c/f_c$ denote the LFM signal bandwidth and the carrier wavelength, respectively.

3. Proposed Algorithm

3.1. Pre-Processing

For BiSAR, the first-order range cell migration (RCM) occupies the main part of the RCM, and the spatial variation is mainly reflected in the first order term. For traditional range alignment methods, such as maximum-correlation method [37], minimum entropy method [38], and global range alignment (GRA) method, these methods can only achieve spatial invariant RCM alignment. Fortunately, keystone transform can achieve spatial variation first-order RCM correction.

However, due to the flexible configuration and motion of BiSAR, it is very easy for marine ship targets to generate large Doppler centroid and Doppler frequency rate, which will cause the problem of Doppler spectrum ambiguity and truncation, resulting in the failure of keystone transform. Therefore, Doppler prefiltering is required to remove Doppler ambiguity and truncation before keystone transform. Using the prior information such as the position and velocity of BiSAR platforms measured by the position system, the Doppler centroid and Doppler frequency rate induced by the BiSAR platform’s motion can be obtained by

(10)$\begin{array}{cc}\hfill f_{dcP}& =-\frac{v_r+v_t}{\lambda }\hfill \end{array}$

(11)$\begin{array}{cc}\hfill f_{drP}& =-\frac{{∥V_T∥}^2-v_t^2}{\lambda ∥r_{To}∥}-\frac{{∥V_R∥}^2-v_r^2}{\lambda ∥r_{Ro}∥}\hfill \end{array}$

where, $r_{Ro}$ and $r_{To}$ denote the distance between the receiver or transmitter and the target at slow time $\xi =0$, respectively.

Therefore, the compensation function can be constructed as

(12)$\begin{array}{c}\hfill H_0\left(f_t,\xi \right)=exp\left\{j\frac{2\pi \left(f_c+f_t\right)}{c}{R}_1\left(\xi \right)\right\}\end{array}$

where

(13)$\begin{array}{c}\hfill R_1\left(\xi \right)=-\lambda f_{dcP}\xi -\frac{1}{2}\lambda f_{drP}{\xi }^2\end{array}$

The compensated echo can be expressed as

(14)$\begin{array}{cc}\hfill S_{c1}\left(f_t,\xi \right)& ={\mathrm{FFT}}_t\left[s_p\left(t,\xi \right)\right]H_0\left(f_t,\xi \right)\hfill \\ \hfill & =\sigma \mathrm{rect}\left[\frac{\xi }{T_a}\right]\mathrm{rect}\left[\frac{f_t}{B_r}\right]exp\left\{-j\frac{2\pi \left(f_c+f_t\right)}{c}\left(R\left(\xi ;{\mathit{r}}_p\right)-R_1\left(\xi \right)\right)\right\}\hfill \end{array}$

After removing the Doppler centroid and Doppler frequency rate induced by BiSAR platforms, the Doppler spectrum has been located in one PRF, and it has been narrowed without truncation. Then, the first-order range migration of the marine ship target can be corrected by using the keystone transformation, which is given by

(15)$\xi =\frac{f_c}{f_c+f_t}{t}_m$

After the keystone transform, the signal becomes

(16)$\begin{array}{c}\hfill S_{c2}\left(f_t,t_m\right)=\sigma \mathrm{rect}\left[\frac{t_m}{T_a}\right]\mathrm{rect}\left[\frac{f_t}{B_r}\right]exp\left\{-j\mathsf{\Phi }\left(f_t,t_m\right)\right\}\end{array}$

where

(17)$\begin{array}{cc}\hfill & \mathsf{\Phi }\left(f_t,t_m\right)=2\pi \frac{f_t+f_c}{c}{R}_0+2\pi \frac{v_{r0}+v_{t0}+{\mathit{\omega }}_0·\mathit{\mu }+\lambda f_{dcP}}{\lambda }{t}_m\hfill \\ \hfill & +2\pi \frac{f_c}{\left(f_c+f_t\right)}\frac{v_r^{\prime }+v_t^{\prime }+{\mathit{\omega }}^{\prime }·\mathit{\mu }+\lambda f_{drP}}{2\lambda }{t}_m^2+2\pi \frac{f_c^2}{{\left(f_c+f_t\right)}^2}\frac{v_r^{″}+v_t^{″}+{\mathit{\omega }}^{″}·\mathit{\mu }}{6\lambda }{t}_m^3\hfill \end{array}$

For an actual BiSAR system, there is usually $f_c\gg f_t$, so it has the following approximation:

(18)$\left\{\begin{array}{c}\frac{f_c}{f_c+f_t}\approx 1-\frac{f_t}{f_c}\hfill \\ {\left(\frac{f_c}{f_c+f_t}\right)}^2\approx 1-2\frac{f_t}{f_c}.\hfill \end{array}\right.$

Substituting (18) into (17) yields

(19)$\begin{array}{cc}\hfill \mathsf{\Phi }\left(f_t,t_m\right)\approx & \frac{2\pi }{\lambda }\left(R\left(t_m;{\mathbf{r}}_p\right)-R_1\left(t_m\right)\right)\hfill \\ \hfill & +\frac{2\pi f_t}{c}\left(R_0-\frac{v_r^{\prime }+v_t^{\prime }+{\mathit{\omega }}^{\prime }·\mathit{\mu }+\lambda f_{drP}}{2}{t}_m^2-\frac{v_r^{″}+v_t^{″}+{\mathit{\omega }}^{″}·\mathit{\mu }}{3\lambda }{t}_m^3\right)\hfill \end{array}$

Compared (14) with (19), it can be seen that the product-term of $t_m$ and $f_t$ has been removed. Therefore, spatially variant first-order range migration of maritime ship targets has been eliminated.

After removing the spatial variation of first-order range migration with keystone transform, there may still be high-order range migration in the echo signal, but its spatial variation is usually less than a range resolution unit, which can be ignored. Therefore, traditional range alignment algorithms can be used to correct spatial invariant high-order range migration. After envelope alignment is completed, the spatial invariant phase error introduced by translation movement will also cause the image to become defocused. Therefore, the translation phase compensation can be realized by traditional methods, such as phase gradient autofocus (PGA) [39].

3.2. Cubic Phase Time-Scaled Transformation

After the pre-processing steps such as Doppler prefiltering, keystone transform, and translation compensation, the signal of each scattering point has been corrected to within one range gate, and combined with the analysis in Section 2, the azimuth signal within one range gate can be represented as a multi-component CPSs. Thus, the azimuth signal within the nth range gate can be denoted as

(20)$\begin{array}{cc}\hfill s_n\left(t_m\right)& =\sum _{i=1}^M{\sigma }_{i}exp\left\{j2\pi \left(b_{i0}+b_{i1}{t}_m+a_{i2}{t}_m^2+b_{i3}{t}_m^3\right)\right\}\hfill \end{array}$

where M is the number of scattering points in the nth range gate. ${\sigma }_i$, $b_{i0}$, $b_{i1}$, $b_{i2}$, and $b_{i3}$ respectively determine the amplitude, constant, first, second, and third order coefficients of the ith scattering point.

To better demonstrate the proposed algorithm, we first consider a mono-CPS within one range gate, which can be expressed as

(21)$s_n^{\left(1\right)}\left(t_m\right)={\sigma }_{1}exp\left\{j2\pi \left(b_{10}+b_{11}{t}_m+b_{12}{t}_m^2+b_{13}{t}_m^3\right)\right\}$

Then, the instantaneous autocorrelation function in the cubic phase function (CPF) [40] of (21) is

(22)$R_s\left(t_m,{\tau }_m\right)=s_n^{\left(1\right)}\left(t_m-{\tau }_m\right)s_n^{\left(1\right)}\left(t_m+{\tau }_m\right)$

where ${\tau }_m$ is the lag time variable. Substituting (21) into (22), we have

(23)$R_s\left(t_m,{\tau }_m\right)=A_{1}exp\left\{j4\pi \left[b_{11}{t}_m+b_{12}(t_m^2+{\tau }_m^2)+b_{13}\left(t_m^3+3t_m{\tau }_m^2\right)\right]\right\}$

where $A_1={\sigma }_1^{2}exp\left\{j4\pi b_{10}\right\}$.

The cubic phase time-scaled transform is defined as

(24)$t_m{\tau }_m^2=\eta \Rightarrow {\tau }_m=\sqrt[]{\frac{\eta }{t_m}}$

where $\eta$ denotes the new time variable after the cubic phase time-scaled transform.

Then, $R_s\left(t_m,{\tau }_m\right)$ can be transformed as

(25)$\begin{array}{c}\hfill R_s\left(t_m,\eta \right)=A_{1}exp\left\{j4\pi \left[b_{11}{t}_m+b_{13}\left(t_m^3+3\eta \right)+b_{12}\left(t_m^2+\eta /t_m\right)\right]\right\}\end{array}$

To remove the coupling between $t_m$ and $\eta$, the compensated phase function is given by

(26)$H_1\left(t_m,\eta ;\zeta \right)=exp\left\{j4\pi \zeta \left(t_m^2+\frac{\eta }{t_m}\right)\right\}$

After multiplying (25) by (26), when $\zeta$ is equal to $b_{12}$, the signal is given by

(27)$\begin{array}{cc}\hfill R_s^{C1}\left(t_m,\eta \right)=& R_s^C\left(t_m,\eta \right)H_1^{\ast }\left(t_m,\eta ;\zeta \right)\hfill \\ \hfill =& exp\left\{j4\pi \left[b_{10}+t_m{b}_{11}+\left(t_m^3+3\eta \right)b_{13}\right]\right\}\hfill \end{array}$

where ∗ denotes the complex conjugation.

Then, the FFT along $\eta$-axis yields

(28)$\begin{array}{cc}\hfill R_s^{C1}\left(t_m,f_{\eta }\right)& ={\mathrm{FFT}}_{\eta }\left[R_s^{C1}\left(t_m,\eta \right)\right]\hfill \\ & =A_0\delta \left(f_{\eta }-6b_{13}\right)exp\left\{j4\pi \left(t_m{b}_{11}+t_m^3{b}_{13}\right)\right\}\hfill \end{array}$

where $f_{\eta }$ denotes the frequency domain relative to $\eta$. ${\mathrm{FFT}}_{\eta }[•]$ denotes the FFT operation along the $\eta$-axis. $\delta (•)$ denotes the Dirac delta function.

In the $t_m-f_{\eta }$ domain, $R_s^{C1}\left(t_m,f_{\eta }\right)$ in (28) is the straight line along $t_m$, so we can accumulate signals along this straight line. However, due to the existence of the cubic phase in (28), if we accumulate directly along $t_m$, it will lead to peak defocusing. To cancel the influence of cubic phase, the sampling properties of the Dirac delta function are used, i.e.,

(29)$\begin{array}{c}\hfill \delta \left(f_{\eta }-f_1\right)g\left(f_{\eta }\right)=g\left(f_1\right)\delta \left(f_{\eta }-f_1\right)\end{array}$

where $g\left(f_{\eta }\right)$ represents a general function of the variable $f_{\eta }$ and $f_1$ indicates a constant values.

According to (29), an appropriate phase-term function $H_2(t_m,f_{\eta })$ is designed to eliminate the cubic phase in (28). Therefore, the phase-term function $H_2(t_m,f_{\eta })$ can be constructed as

(30)$H_2\left(t_m,f_{\eta }\right)=exp\left\{-j\frac{2\pi }{3}{f}_{\eta }{t}_m^3\right\}$

After multiplying (28) by (30), and

(31)$\begin{array}{cc}\hfill R_s^{C2}\left(t_m,f_{\eta }\right)& =R_s^{C1}\left(t_m,f_{\eta }\right)H_2\left(t_m,f_{\eta }\right)\hfill \\ \hfill & ={\left.g\left(f_{\eta }\right)\right|}_{f_{\eta }=6b_{13}}\delta \left(f_{\eta }-6b_{13}\right)exp\left\{j4\pi t_m{b}_{11}\right\}\hfill \end{array}$

where $g\left(f_{\eta }\right)=exp\left\{j2\pi \left(f_{\eta }-6b_{13}\right)t_m^3/3\right\}$.

Performing the FFT with respect to $t_m$, the signal will be integrated as a well-focused peak in $f_{\eta }$ and $f_{t_m}$.

(32)$\begin{array}{cc}\hfill R_s^{C2}\left(f_{t_m},f_{\eta }\right)& ={\mathrm{FFT}}_{t_m}\left[R_s^{C2}\left(t_m,f_{\eta }\right)\right]\hfill \\ \hfill & ={\left.g\left(f_{\eta }\right)\right|}_{f_{\eta }=6b_{13}}\delta \left(f_{\eta }-6b_{13}\right)\delta \left(f_{t_m}-2b_{11}\right)\hfill \end{array}$

where $f_{t_m}$ denotes the frequency domain relative to $t_m$.

Only when $\zeta$ accurately matches the real second-order coefficient (i.e., $\zeta =b_{12}$), can the energy of $R_s^{C2}\left(f_{t_m},f_{\eta }\right)$ get the maximum accumulation. When $\zeta$ is mismatched (i.e., $\zeta \ne b_{12}$), the energy of $R_s^{C2}\left(f_{t_m},f_{\eta }\right)$ will be dispersed. Therefore, the estimated value of $b_{12}$ can be given by

(33)${\widehat{b}}_{12}=arg\underset{\zeta }{max}\left|R_s^{C2}(f_{t_m},f_{\eta })\right|$

Furthermore, based on the peak location in (32), $b_{11}$ and $b_{13}$ can be evaluated as

(34)${\widehat{b}}_{11}={\widehat{f}}_{t_m}/2,{\widehat{b}}_{13}={\widehat{f}}_{\eta }/6$

3.3. Main Procedure of the Proposed Algorithm

Based on the estimation results of Doppler parameter in Section 3.2, a bistatic SAR maritime ship target imaging method based on cubic phase time-scaled transformation is established, and its detailed flow is illustrated in the Figure 2. The main procedure is as follows:

Step 1: After the range compression echo of the BiSAR maritime ship target, the Doppler prefiltering and the keystone transform are applied to the translation compensation to complete the echo pre-processing, and get the signal $s_n$ of the nth range gate of the echo.

Step 2: Initialize $n=1$ and $i=1$.

Step 3: Calculate the cubic phase time-scaled transformation of signal $s_n$, and obtain $b_{i1}$, $b_{i2}$, and $b_{i3}$ of signal $s_n^{\left(i\right)}$ through phase compensation and FFT processing.

Step 4: Construct the ith component signal $s_n^{\left(i\right)}=\widehat{A}exp\left(j2\pi \left(b_{i1}{t}_m+b_{i2}{t}_m^2+b_{i3}{t}_m^3\right)\right)$, where, $\widehat{A}=max\left(\mathrm{FFT}\left(s_n·exp\left(-j2\pi \left(b_{i1}{t}_m+b_{i2}{t}_m^2+b_{i3}{t}_m^3\right)\right)\right)\right)/M$, and update $s_n=s_n-s_n^{\left(i\right)}$ and $i=i+1$.

Step 5: Repeat Steps 3 and 4 up to the energy signal of $s_n$, which is under the threshold ${\zeta }_{th}$ (generally 0.1 of the original signal). We then update $n=n+1$.

Step 6: Repeat Steps 3–5 until all range gate signals are completed.

By using the proposed algorithm to process all range gate signals, the imaging results of BiSAR maritime ship targets can be obtained.

4. Performance Analysis of Parameter Estimation

In this simulation, mono-CPS and multi-CPS are considered, and the pulse repetition frequency and the azimuth pulse number are 300 Hz and 1024, respectively. The signal parameters of CPS${}_1$ and CPS${}_2$ are shown in Table 1.

4.1. Mono-CPS

First, the mono-CPS with noise free is considered. Figure 3a,b show respectively the amplitude and phase of the auto-correlation result $R_s\left(t_m,{\tau }_m\right)$ in the $t_m-{\tau }_m$ domain for the CPS${}_1$. After the cubic phase time-scaled transform, the amplitude and phase of $R_s\left(t_m,\eta \right)$ are respectively shown in Figure 3c,d. With the FFT operation along the $\eta$-axis of $R_s\left(t_m,\eta \right)$, Figure 3e gives the time-frequency distribution, and can be expressed by

(35)$\begin{array}{c}\hfill R_s\left(t_m,f_{\eta }\right)=A_1\delta \left(f_{\eta }-\left(6b_{13}+\frac{2b_{12}}{t_m}\right)\right)exp\left\{j4\pi \left[b_{11}{t}_m+b_{12}{t}_m^2+b_{13}{t}_m^3\right]\right\}\end{array}$

Due to the coupling term $\eta /t_m$ between the $t_m$ and the $\eta$ in $R_s^C(t_m,\eta )$, $t_m$ is inversely proportional to $f_{\eta }$, as shown in Figure 3e and Equation (35). Figure 3f shows the searching result of $b_{12}$ according to (33). Obviously, depending on the peak location, the estimated value ${\widehat{b}}_{12}$ can be obtained. In Figure 3f, ${\widehat{b}}_{12}$ is equal to 4 Hz/s. After the first phase compensation, the coupling term is completely removed. Thus, with the FFT operation with respect to $\eta$ of $R_s^{C1}(t_m,\eta )$, the rectangular hyperbolas are corrected into vertical straight lines, as shown in Figure 4a. After the phase-term function $H_2(t_m,f_{\eta })$ and performing the FFT with respect to $t_m$ of Figure 4a, we obtain Figure 4b, where the signal accumulates into a peak in $f_{\eta }$ and $f_{t_m}$. Figure 4c gives the stereogram of the final focusing result. In Figure 4c, by the peak detection technique, ${\widehat{b}}_{11}$ and ${\widehat{b}}_{13}$ are equal to 6 Hz and 2 ${\mathrm{Hz}/\mathrm{s}}^2$, respectively.

Next, a mono-CPS with Gaussian noise condition is considered, and we set the signal-to-noise rate (SNR) to −5 dB, −7 dB, and −10 dB, respectively. Figure 5a–c shows a stereogram of the final focusing result using the SCFT-based algorithm [31] under −5 dB, −7 dB, and −10 dB of SNR. It can be seen from Figure 5a–c, that the CPS is obviously submerged in noise due to the fourth-order nonlinearity of its algorithm. Figure 6a–c shows the stereogram of the final focusing result $R_s^{C2}\left(f_{t_m},f_{\eta }\right)$ by the proposed algorithm under −5 dB, −7 dB, and −10 dB of SNR. As shown in Figure 6a–c, the CPS can still be finely focused. Clearly, it can be seen that the proposed algorithm has a better noise immunity than the SCFT-based algorithm.

4.2. Multi-CPS

For multi-component CPS, since this algorithm uses the bilinear transform of the auto-correlation function in (22), cross-terms will inevitably occur. Next, we will analyze the impact of cross-terms on parameter estimation in detail. After substituting the multi-CPS in (20) into (22) and performing the cubic phase time-scaled transform of (24), we have

(36)$\begin{array}{l}{R}_{multi}\left(t_m,\eta \right)=R_{auto}\left(t_m,\eta \right)+R_{cross}\left(t_m,\eta \right),\\ R_{auto}\left(t_m,\eta \right)=\sum _{i=1}^{K}exp\left\{j4\pi \left[b_{i1}{t}_m+b_{i3}\left(t_m^3+3\eta \right)+b_{i2}\left(t_m^2+\eta /t_m\right)\right]\right\},\\ R_{cross}\left(t_m,\eta \right)=\sum _{i,j=1,i\ne j}^{K}exp\left\{j2\pi \right[b_{i0}+b_{j0}+\left(b_{i1}+b_{j1}\right)t_m+\left(b_{i2}+b_{j2}\right)\left(t_m^2+\eta /t_m\right)+\left(b_{i3}+b_{j3}\right)\left(t_m^3+3\eta \right)\\ \left.\left.+\left(b_{i1}-b_{j1}+2\left(b_{i2}-b_{j2}\right)t_m+3\left(b_{i3}-b_{j3}\right)t_m^2\right)\sqrt{\frac{\eta }{t_m}}+(b_{i3}-b_{j3}){\left(\frac{\eta }{t_m}\right)}^{3/2}\right]\right\}.\end{array}$

With the proper second-order phase compensation and FFT operation with respect to $\eta$ of $R_{multi}\left(t_m,\eta \right)$, the auto-term in (36) is the two straight lines and symmetric $t_m=0$ axis. However, for the cross-term in (36), in addition to the residual coupling term $\eta /t_m$, there are also coupling terms $\sqrt{\eta /t_m}$ and ${\left(\eta /t_m\right)}^{3/2}$, so the cross-term is dispersed in the two-dimensional time-frequency in the form of curves. Then, after the second phase compensation and Fourier transform along the direction of $t_m$, the two straight lines of the auto-term accumulate into a point, while the cross-term does not achieve effective accumulation and is still dispersed in the two-dimensional frequency domain. Therefore, the cross-term will not affect the peak detection of the auto-term. Therefore, the cross-terms do not impact the accumulation of auto-term with peak detection.

Next, the parameter estimation of multi-component CPS is simulated and analyzed. We consider noise-free CPS${}_1$ and CPS${}_2$ in the Table 1. With the FFT operation along $\eta$-axis of $R_s\left(t_m,\eta \right)$, Figure 7a gives the time-frequency distribution of CPS $s_1\left(t_m\right)$ and $s_2\left(t_m\right)$. As shown in Figure 7, the two-dimensional time-frequency distribution results contain the auto-term and cross-term, and it can be seen from $R_{multi}\left(t_m,\eta \right)$ in (36) that the auto-term is rectangular hyperbola in the $t_m-f_{\eta }$ domain, while the cross-term is a curved curve, as shown in Figure 7a.

By searching the $b_{12}$, the results are shown in Figure 7b. The searched coefficients ${\widehat{b}}_{12}$ are used to compensate the cubic phase time-scaled signal. The auto-term of the compensated signal becomes two straight lines in the $t_m-f_{\eta }$ domain and is symmetric about the $t_m$ axis. The other cross-terms are still scattered in the form of curves, as shown in Figure 8a. Then, through phase compensation and FFT processing along $t_m$-axis, the two lines of the auto-term are accumulated into a peak, so that the CPS signal can be accurately focused in the two-dimensional frequency domain, as shown in Figure 8b. In order to highlight the details, the lower left corner with a cyan border is a partial enlarged view at the peak location. As for the cross-term, because it is scattered in the two-dimensional time-frequency domain and is not symmetrical about $t_m$, the cross-term cannot achieve effective accumulation in the two-dimensional frequency domain. The accumulation result of the CPS${}_1$ is shown in Figure 8b,c. After removing the CPS${}_1$ signal with CLEAN technology [41], the results of phase compensation and FFT processing on CPS${}_2$ are shown in Figure 8d–f. As shown in Figure 8d, the signal energy of CPS${}_1$ has been basically removed, and the auto-term has become two straight lines after compensating the second-order term coefficient and FFT along the direction of $t_m$. Then, through phase compensation and FFT processing along $t_m$-axis, the two lines of the auto-term are accumulated into a peak, so that the CPS${}_2$ can be accurately focused in the two-dimensional frequency domain, as shown in Figure 8e,f.

By comparing Figure 8c,f, it can be found that their differences are small, which also means that the cross-terms in the multi-component CPS will not be effectively accumulated by the proposed algorithm, and the interference to signal accumulation and peak detection is relatively small.

4.3. Computational Complexity

The SCFT-based algorithm [31] can be computed by the 2-D chirp-z transform (CZT), which only needs complex multiplications and the fast Fourier transforms (FFTs) [42], and its computational complexity is $O\left(N^{2}logN\right)$, where, N represents the signal length and number of delay variables.

In the proposed algorithm, in addition to the phase compensation and the fast Fourier transform (FFT), there are cubic phase time-scaled transformation. Among them, cubic phase time-scaled transformation can be implemented by the scaling principle [35], and only complex multiplication and FFT/IFFT operations are required. Therefore, in a search process, it only needs complex multiplication and FFT/IFFT operations, and its computational complexity is $O\left(N^{2}logN\right)$. Due to the presence of one-dimensional search, its computation time is longer than that of the SCFT-based algorithm. However, the image quality comparison from the simulation analysis in Section 5 shows that its imaging results are significantly better than the SCFT-based algorithm.

Moreover, compared with the multi-dimensional brute-force search algorithm with the computational complexity of $O\left(N^3\right)$, such as ML [29] and DCFT [30], the proposed algorithm only needs the one-dimensional search of the second-order term to realize the coherent accumulation of the cubic phase signal in the two-dimensional frequency domain, which ensures the accuracy of the parameter estimation and also has a small computational complexity, which is convenient for engineering applications.

5. Simulation

To verify the validity of the algorithm proposed in this paper, a maritime ship target is simulated in this section to emphasize the performance of the algorithm. As shown in Figure 9, a maritime ship target is assumed to contain several scattering points with a width and length of 20 and 120 m. The relevant simulation parameters of BiSAR and the maritime ship target are shown in Table 2.

In Figure 10 and Figure 11, the effectiveness of the pre-processing using the Doppler prefiltering, keystone transform, and translation compensation is first verified. Figure 10 illustrates the range-compressed echo, its 2-D spectrum, and result using the RD algorithm without pre-processing. Figure 10a shows the range-compressed echo of BiSAR maritime ship targets, from which it can be seen that there is a large RCM, and the RCM is spatial variation. Its two-dimensional spectrum is shown in Figure 10b. It can be seen that its Doppler spectrum has problems with Doppler ambiguity and truncation. When keystone transform is used to eliminate the spatial variation of RCM, Doppler prefiltering is required in advance to remove the Doppler ambiguity and truncation. The image obtained using the RD algorithm is severely defocused along the range and azimuth direction, as shown in Figure 10c.

In Figure 11, the pre-processing including the Doppler prefiltering, keystone transform, and translation compensation are performed. It can be seen in Figure 11a that the echo of each scattering point can be corrected to its corresponding one range gate. Additionally, the Doppler spectrum has been corrected to a PRF after Doppler prefiltering, as shown in Figure 11b. The image obtained using the RD algorithm only defocused along the azimuth direction, as shown in Figure 11c. Therefore, it can be concluded that it is effective to compensate the migration of bistatic SAR maritime ship targets through this pre-processing step.

Next, the comparison results of different methods will be given. The maritime ship target imaging results obtained by the RD algorithm are depicted in Figure 12, due to the spatially-variant Doppler parameter including Doppler frequency rate and third-order Doppler frequency rate caused by the three-dimensional rotation of the target, the scattering points of the ship target in the RD algorithm result have defocus, and the degree of defocus of the different scattering points is different, as shown in Figure 12.

The imaging result obtained by the SCFT-based algorithm is shown in Figure 13 for three SNRs of SNR = 0, −5 and −10 dB. Compared with the RD algorithm, the image quality is improved, but the details in the ship target are smeared due to the influence of cross items. In addition, with the reduction of SNR, the imaging quality of the maritime ship target is degraded due to noise interference. In the case of SNR = −10 dB, the complete shape of the ship can no longer be obtained, as shown in Figure 13c. By contrast, the imaging result obtained by proposed algorithm is shown in Figure 14 for three SNRs of SNR = 0, −5, and −10 dB. When SNR = 0, −5, and −10 dB, the complete shape of the maritime ship target and good details can be obtained, with better imaging quality.

Image entropy (IE) and image sharpness (IS) [43,44] are used to quantitatively assess the imaging quality, which can be expressed respectively as

(37)$\begin{array}{cc}\hfill IE& =-\sum _{m=1}^M\sum _{n=1}^N\frac{{\left|I\left(m,n\right)\right|}^2}{S}\ln \frac{{\left|I\left(m,n\right)\right|}^2}{S}\hfill \\ \hfill IS& =\sum _{m=1}^M\sum _{n=1}^N\left[-\ln \left({\left|T\left(m,n\right)\right|}^2+1\right)\right]\hfill \\ \hfill T\left(m,n\right)& =I\left(m,n\right)\sqrt{\frac{MN}{S}},S=\sum _{m=1}^M\sum _{n=1}^N{\left|I\left(m,n\right)\right|}^2\hfill \end{array}$

where $I\left(m,n\right)$ denotes the images, m and n are the index of azimuth and range, respectively. M and N are the number of azimuth and range, respectively.

When the IE is smaller or the IS is larger, better quality and clearer image can be obtained. The entropies and sharpness of different images obtained from different methods are listed respectively in Table 3 and Table 4 to demonstrate the superior performance of the proposed algorithm. This improvement is especially evident when the SNR is low, as shown in the last row of Table 3 and Table 4. Furthermore, compared with the RD algorithm and SCFT-based algorithm, under the same SNR, the ship target image obtained by the proposed algorithm has the largest sharpness, lowest entropy, and the best image quality. In the case of different SNR, the difference of image sharpness of the ship targets obtained by the proposed algorithm is the smallest, indicating that the proposed algorithm is less affected by noise.

6. Conclusions

In this paper, an novel BiSAR maritime ship target imaging algorithm based on cubic phase time-scaled transformation is proposed. First, the Doppler prefiltering is performed on the range-compressed echo to remove the Doppler ambiguity and truncation. Second, keystone transform and GRA method are performed to remove the spatially-variant first-order range migration and accomplish the translation compensation. Then, the azimuth signal is modeled as a multi-component cubic phase signal. Based on the proposed cubic phase time-scaled transformation, the coupling term relationship between time and lag time variables is changed to a hyperbola to make cross-terms more diffuse and reduce its impact on the peak of auto-term. After a one-dimensional search of second-order term coefficients, two-dimensional coherent integration can be achieved by phase compensation and 2-D FFT, so as to estimate the Doppler centroid, frequency rate, and third-order frequency of scattering points. After all the range gates are processed, the high resolution image of the maritime ship target in BiSAR can be obtained. Verified by the simulated data, the algorithm has excellent noise immunity and low cross-terms, and can achieve fine focusing maritime ship target in BiSAR compared with tradition RD algorithm and the SCFT-based algorithm.

Author Contributions

Conceptualization, Q.Y. and Z.L.; methodology, Q.Y. and Z.L.; writing—original draft preparation, Q.Y. and J.L.; writing—review and editing, Q.Y. and H.A.; project administration, Y.P. and J.Y.; funding acquisition, H.A., Z.L. and J.W. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Geometry of the BiSAR platforms and a maritime ship target.

Figure 2. Flowchart of the proposed imaging algorithm.

Figure 3. Results of a noise—free mono—CPS. (a) Amplitude of [Forumla omitted. See PDF.]. (b) Phase of [Forumla omitted. See PDF.]. (c) Amplitude of [Forumla omitted. See PDF.] in [Forumla omitted. See PDF.] domain. (d) Phase of [Forumla omitted. See PDF.] in [Forumla omitted. See PDF.] domain. (e) Amplitude of [Forumla omitted. See PDF.] in [Forumla omitted. See PDF.] domain. (f) Searching result of [Forumla omitted. See PDF.].

Figure 4. Results of a noise–free mono–CPS. (a) [Forumla omitted. See PDF.] in [Forumla omitted. See PDF.] domain after cubic phase time-scaled transform and phase compensation. (b) Final focusing result [Forumla omitted. See PDF.]. (c) Stereogram of the final focusing result [Forumla omitted. See PDF.].

Figure 5. Stereogram of the final focusing result of the SCFT-based algorithm. (a) SNR = −5 dB. (b) SNR = −7 dB. (c) SNR = −10 dB.

Figure 6. Stereogram of the final focusing result of the proposed algorithm. (a) SNR = −5 dB. (b) SNR = −7 dB. (c) SNR = −10 dB.

Figure 7. Results of a multi-CPS. (a) [Forumla omitted. See PDF.] in [Forumla omitted. See PDF.] domain with cubic phase time-scaled transform. (b) The searching result of [Forumla omitted. See PDF.].

Figure 8. Results of a multi-CPS. (a) [Forumla omitted. See PDF.] in [Forumla omitted. See PDF.] domain with cubic phase time-scaled transform and phase compensation. (b) Final focusing result [Forumla omitted. See PDF.]. (c) Stereogram of the final focusing result [Forumla omitted. See PDF.]. (d) [Forumla omitted. See PDF.] in [Forumla omitted. See PDF.] domain after removing CPS[Forumla omitted. See PDF.]. (e) Final focusing result [Forumla omitted. See PDF.] after removing CPS[Forumla omitted. See PDF.]. (f) Stereogram of the final focusing result [Forumla omitted. See PDF.] after removing CPS[Forumla omitted. See PDF.].

Figure 9. Model of the maritime ship target.

Figure 10. Result without pre-processing. (a) Range-compressed echo. (b) 2-D spectrum of range-compressed echo. (c) RD result.

Figure 11. Result with pre-processing. (a) Range-compressed echo. (b) 2-D spectrum of range-compressed echo. (c) RD result.

Figure 12. Maritime ship target imaging results obtained by the RD algorithm. (a) SNR = 0 dB. (b) SNR = −5 dB. (c) SNR = −10 dB.

Figure 13. Maritime ship target imaging results obtained by the SCFT−based algorithm [31]. (a) SNR = 0 dB. (b) SNR = −5 dB. (c) SNR = −10 dB.

Figure 14. Maritime ship target imaging results obtained by the proposed algorithm. (a) SNR = 0 dB. (b) SNR = −5 dB. (c) SNR = −10 dB.

Table 1

Signal parameters of CPSs.

	Amplitude ${\mathit{\sigma }}_{\mathit{i}}$	Constant Coefficient ${\mathit{b}}_{\mathit{i}0}$	First-Order Coefficient ${\mathit{b}}_{\mathit{i}1}$	Second-Order Coefficient ${\mathit{b}}_{\mathit{i}2}$	Third-Order Coefficient ${\mathit{b}}_{\mathit{i}3}$
CPS1	1	0.25	6	4	2
CPS2	1	0.6	10	−3	1

Table 2

Simulation parameters.

Parameters	Value	Parameters	Value
Carrier Frequency $f_c$	10 GHz	Platform Velocity V	30 m/s
Signal Bandwidth $B_r$	200 MHz	Velocity $\mathit{v}$	(6, −3, 0) m/s
Pulse Repetition Frequency	1000 Hz	Angular velocity $\mathit{\omega }$	0.04 rad/s
Coordinate of Receiver ${\overrightarrow{P}}_R$	(−8000, 0, 8000) m	Angular acceleration ${\mathit{\omega }}^{\prime }$	0.015 rad/s${}^2$
Coordinate of Transmitter ${\overrightarrow{P}}_T$	(0, 0, 8000) m	Angular acceleration rate ${\mathit{\omega }}^{″}$	0.01 rad/s${}^3$

Table 3

Entropy of the images.

	SNR = 0	SNR = −5	SNR = −10
RD algorithm	9.5216	10.0347	10.6671
The SCFT-based algorithm [31]	8.0105	8.2476	9.3700
The proposed algorithm	7.6766	7.7101	7.8290

Table 4

Sharpness of the images.

	SNR = 0	SNR = −5	SNR = −10
RD algorithm	$-3.5768\times {10}^4$	$-4.5335\times {10}^4$	$-5.6726\times {10}^4$
The SCFT-based algorithm [31]	$-1.8120\times {10}^4$	$-2.0681\times {10}^4$	$-3.5381\times {10}^4$
The proposed algorithm	$-1.5487\times {10}^4$	$-1.5613\times {10}^4$	$-1.6417\times {10}^4$

References

1. Sommer, A.; Ostermann, J. Backprojection Subimage Autofocus of Moving Ships for Synthetic Aperture Radar. IEEE Trans. Geosci. Remote Sens.; 2019; 57, pp. 8383-8393. [DOI: https://dx.doi.org/10.1109/TGRS.2019.2920779]

2. Curlander, J.C.; McDonough, R.N. Synthetic Aperture Radar: Systems and Signal Processing, 1 ed.; Wiley-Interscience: New York, NY, USA, 1991.

3. Cumming, I.G.; Wong, F.H. Digital Processing of Synthetic Aperture Radar Data; Artech House: London, UK, 2005; Volume 1, pp. 108-110.

4. Li, Z.; Wu, J.; Yang, J.; Liu, Z. Bistatic SAR Clutter Suppression: Theory, Method, and Experiment; Springer Nature: Berlin/Heidelberg, Germany, 2022.

5. Walterscheid, I.; Ender, J.H.G.; Brenner, A.R.; Loffeld, O. Bistatic SAR Processing and Experiments. IEEE Trans. Geosci. Remote Sens.; 2006; 44, pp. 2710-2717. [DOI: https://dx.doi.org/10.1109/TGRS.2006.881848]

6. Wu, J.; Li, Z.; Huang, Y.; Yang, J.; Liu, Q.H. An Omega-K Algorithm for Translational Invariant Bistatic SAR Based on Generalized Loffeld’s Bistatic Formula. IEEE Trans. Geosci. Remote Sens.; 2014; 52, pp. 6699-6714. [DOI: https://dx.doi.org/10.1109/TGRS.2014.2301433]

7. Wu, J.; Li, Z.; Huang, Y.; Yang, J.; Liu, Q.H. Omega-K Imaging Algorithm for One-Stationary Bistatic SAR. IEEE Trans. Aerosp. Electron. Syst.; 2014; 50, pp. 33-52. [DOI: https://dx.doi.org/10.1109/TAES.2013.120399]

8. Wu, J.; Sun, Z.; Li, Z.; Huang, Y.; Yang, J.; Liu, Z. Focusing Translational Variant Bistatic Forward-Looking SAR Using Keystone Transform and Extended Nonlinear Chirp Scaling. Remote Sens.; 2016; 8, 840. [DOI: https://dx.doi.org/10.3390/rs8100840]

9. Xi, Z.; Duan, C.; Zuo, W.; Li, C.; Huo, T.; Li, D.; Wen, H. Focus Improvement of Spaceborne-Missile Bistatic SAR Data Using the Modified NLCS Algorithm Based on the Method of Series Reversion. Remote Sens.; 2022; 14, 5770. [DOI: https://dx.doi.org/10.3390/rs14225770]

10. Li, Z.; Huang, C.; Sun, Z.; An, H.; Wu, J.; Yang, J. BeiDou-Based Passive Multistatic Radar Maritime Moving Target Detection Technique via Space–Time Hybrid Integration Processing. IEEE Trans. Geosci. Remote Sens.; 2022; 60, pp. 1-13. [DOI: https://dx.doi.org/10.1109/TGRS.2021.3128650]

11. Li, Z.; Li, S.; Liu, Z.; Yang, H.; Wu, J.; Yang, J. Bistatic Forward-Looking SAR MP-DPCA Method for Space–Time Extension Clutter Suppression. IEEE Trans. Geosci. Remote Sens.; 2020; 58, pp. 6565-6579. [DOI: https://dx.doi.org/10.1109/TGRS.2020.2977982]

12. Zhang, H.; Tang, J.; Wang, R.; Deng, Y.; Wang, W.; Li, N. An Accelerated Backprojection Algorithm for Monostatic and Bistatic SAR Processing. Remote Sens.; 2018; 10, 140. [DOI: https://dx.doi.org/10.3390/rs10010140]

13. Sun, H.; Sun, Z.; Chen, T.; Miao, Y.; Wu, J.; Yang, J. An Efficient Backprojection Algorithm Based on Wavenumber-Domain Spectral Splicing for Monostatic and Bistatic SAR Configurations. Remote Sens.; 2022; 14, 1885. [DOI: https://dx.doi.org/10.3390/rs14081885]

14. Sun, Z.; Chen, T.; Sun, H.; Wu, J.; Lu, Z.; Li, Z.; An, H.; Yang, J. A Novel Frequency-Domain Focusing Method for Geosynchronous Low-Earth-Orbit Bistatic SAR in Sliding-Spotlight Mode. Remote Sens.; 2022; 14, 3178. [DOI: https://dx.doi.org/10.3390/rs14133178]

15. Li, Z.; Wu, J.; Huang, Y.; Yang, H.; Yang, J. An Adaptive Moving Target Imaging Method for Bistatic Forward-Looking SAR Using Keystone Transform and Optimization NLCS. Sensors; 2017; 17, 216. [DOI: https://dx.doi.org/10.3390/s17010216] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28125000]

16. Yang, Q.; Li, Z.; Li, J.; Xiao, Y.; An, H.; Wu, J.; Pi, Y.; Yang, J. An Optimal Polar Format Refocusing Method for Bistatic SAR Moving Target Imaging. IEEE Trans. Geosci. Remote Sens.; 2022; 60, pp. 1-17. [DOI: https://dx.doi.org/10.1109/TGRS.2022.3176402]

17. Li, X.; Jintao, X.; Yulin, H.; Jianyu, Y. Research on airborne bistatic SAR squint imaging mode algorithm and experiment data processing. Proceedings of the 2007 1st Asian and Pacific Conference on Synthetic Aperture Radar; Huangshan, China, 5 November 2007; pp. 618-621. [DOI: https://dx.doi.org/10.1109/APSAR.2007.4418688]

18. Yang, J.; Huang, Y.; Yang, H.; Wu, J.; Li, W.; Li, Z.; Yang, X. A first experiment of airborne bistatic forward-looking SAR-Preliminary results. Proceedings of the 2013 IEEE International Geoscience and Remote Sensing Symposium-IGARSS, IEEE; Melbourne, VIC, Australia, 21 July 2013; pp. 4202-4204. [DOI: https://dx.doi.org/10.1109/IGARSS.2013.6723760]

19. Li, Z.; Ye, H.; Liu, Z.; Sun, Z.; An, H.; Wu, J.; Yang, J. Bistatic SAR Clutter-Ridge Matched STAP Method for Nonstationary Clutter Suppression. IEEE Trans. Geosci. Remote Sens.; 2022; 60, pp. 1-14. [DOI: https://dx.doi.org/10.1109/TGRS.2021.3125043]

20. Li, Z.; Wu, J.; Huang, Y.; Sun, Z.; Yang, J. Ground-Moving Target Imaging and Velocity Estimation Based on Mismatched Compression for Bistatic Forward-Looking SAR. IEEE Trans. Geosci. Remote Sens.; 2016; 54, pp. 3277-3291. [DOI: https://dx.doi.org/10.1109/TGRS.2016.2514494]

21. Li, Z.; Wu, J.; Liu, Z.; Huang, Y.; Yang, H.; Yang, J. An Optimal 2-D Spectrum Matching Method for SAR Ground Moving Target Imaging. IEEE Trans. Geosci. Remote Sens.; 2018; 56, pp. 5961-5974. [DOI: https://dx.doi.org/10.1109/TGRS.2018.2829166]

22. Li, J.; Li, Z.; Yang, Q.; Wang, Y.; Long, J.; Wu, J.; Xia, W.; Yang, J. Joint Clutter Suppression and Moving Target Indication in 2-D Azimuth Rotated Time Domain for Single-Channel Bistatic SAR. IEEE Trans. Geosci. Remote Sens.; 2023; 61, pp. 1-16. [DOI: https://dx.doi.org/10.1109/TGRS.2023.3237553]

23. Chen, V.C.; Qian, S. Joint time-frequency transform for radar range-Doppler imaging. IEEE Trans. Aerosp. Electron. Syst.; 1998; 34, pp. 486-499. [DOI: https://dx.doi.org/10.1109/7.670330]

24. Wood, J.C.; Barry, D.T. Linear signal synthesis using the Radon-Wigner transform. IEEE Trans. Signal Process.; 1994; 42, pp. 2105-2111. [DOI: https://dx.doi.org/10.1109/78.301845]

25. Sun, Y.; Willett, P. Hough transform for long chirp detection. IEEE Trans. Aerosp. Electron. Syst.; 2002; 38, pp. 553-569. [DOI: https://dx.doi.org/10.1109/TAES.2002.1008986]

26. Cohen, L. Time-frequency distributions-a review. Proc. IEEE; 1989; 77, pp. 941-981. [DOI: https://dx.doi.org/10.1109/5.30749]

27. Wang, Y.; Liu, Q. Super-Resolution Sparse Aperture ISAR Imaging of Maneuvering Target via the RELAX Algorithm. IEEE Sensors J.; 2018; 18, pp. 8726-8738. [DOI: https://dx.doi.org/10.1109/JSEN.2018.2868308]

28. Li, Z.; Zhang, X.; Yang, Q.; Xiao, Y.; An, H.; Yang, H.; Wu, J.; Yang, J. Hybrid SAR-ISAR Image Formation via Joint FrFT-WVD Processing for BFSAR Ship Target High-Resolution Imaging. IEEE Trans. Geosci. Remote Sens.; 2022; 60, pp. 1-13. [DOI: https://dx.doi.org/10.1109/TGRS.2021.3117280]

29. Abatzoglou, T.J. Fast Maximum Likelihood Joint Estimation of Frequency and Frequency Rate. IEEE Trans. Aerosp. Electron. Syst.; AES-22, pp. 708-715. [DOI: https://dx.doi.org/10.1109/TAES.1986.310805]

30. 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, pp. 4201-4212. [DOI: https://dx.doi.org/10.1109/TGRS.2012.2189220]

31. 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, pp. 1418-1429. [DOI: https://dx.doi.org/10.1109/TGRS.2013.2251348]

32. Zheng, J.; Su, T.; Zhang, L.; Zhu, W.; Liu, Q.H. ISAR Imaging of Targets With Complex Motion Based on the Chirp Rate–Quadratic Chirp Rate Distribution. IEEE Trans. Geosci. Remote Sens.; 2014; 52, pp. 7276-7289. [DOI: https://dx.doi.org/10.1109/TGRS.2014.2310474]

33. Zheng, J.; Su, T.; Zhu, W.; Zhang, L.; Liu, Z.; Liu, Q.H. ISAR Imaging of Nonuniformly Rotating Target Based on a Fast Parameter Estimation Algorithm of Cubic Phase Signal. IEEE Trans. Geosci. Remote Sens.; 2015; 53, pp. 4727-4740. [DOI: https://dx.doi.org/10.1109/TGRS.2015.2408350]

34. 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. Topics Appl. Earth Observ. Remote Sens.; 2016; 9, pp. 5688-5702. [DOI: https://dx.doi.org/10.1109/JSTARS.2016.2543233]

35. Zhu, D.; Li, Y.; Zhu, Z. A Keystone Transform Without Interpolation for SAR Ground Moving-Target Imaging. IEEE Geosci. Remote Sens. Lett.; 2007; 4, pp. 18-22. [DOI: https://dx.doi.org/10.1109/LGRS.2006.882147]

36. Murray, R.M.; Li, Z.; Sastry, S.S. A Mathematical Introduction to Robotic Manipulation; CRC Press: Boca Raton, FL, USA, 2017.

37. Chen, C.C.; Andrews, H.C. Target-Motion-Induced Radar Imaging. IEEE Trans. Aerosp. Electron. Syst.; 1980; AES-16, pp. 2-14. [DOI: https://dx.doi.org/10.1109/TAES.1980.308873]

38. Xi, L.; Guosui, L.; Ni, J. Autofocusing of ISAR images based on entropy minimization. IEEE Trans. Aerosp. Electron. Syst.; 1999; 35, pp. 1240-1252. [DOI: https://dx.doi.org/10.1109/7.805442]

39. Wahl, D.E.; Eichel, P.H.; Ghiglia, D.C.; Jakowatz, C.V. Phase gradient autofocus-a robust tool for high resolution SAR phase correction. IEEE Trans. Aerosp. Electron. Syst.; 1994; 30, pp. 827-835. [DOI: https://dx.doi.org/10.1109/7.303752]

40. O’Shea, P. A new technique for instantaneous frequency rate estimation. IEEE Signal Process. Lett.; 2002; 9, pp. 251-252. [DOI: https://dx.doi.org/10.1109/LSP.2002.803003]

41. Martorella, M.; Acito, N.; Berizzi, F. Statistical CLEAN Technique for ISAR Imaging. IEEE Trans. Geosci. Remote Sens.; 2007; 45, pp. 3552-3560. [DOI: https://dx.doi.org/10.1109/TGRS.2007.897440]

42. Rabiner, L.R.; Schafer, R.W.; Rader, C.M. The chirp z-transform algorithm and its application. Bell Syst. Tech. J.; 1969; 48, pp. 1249-1292. [DOI: https://dx.doi.org/10.1002/j.1538-7305.1969.tb04268.x]

43. Schulz, T.J. Optimal Sharpness Function for SAR Autofocus. IEEE Signal Process. Lett.; 2007; 14, pp. 27-30. [DOI: https://dx.doi.org/10.1109/LSP.2006.881525]

44. Yang, W.; Guo, J.; Chen, J.; Liu, W.; Deng, J.; Wang, Y.; Zeng, H. A Novel Channel Inconsistency Estimation Method for Azimuth Multichannel SAR Based on Maximum Normalized Image Sharpness. IEEE Trans. Geosci. Remote Sens.; 2022; 60, pp. 1-16. [DOI: https://dx.doi.org/10.1109/TGRS.2022.3219818]