A Novel Method for Refocusing Moving Ships in SAR Images via ISAR Technique

Abstract

As an active microwave remote sensing device, synthetic aperture radar (SAR) has been widely used in the field of marine surveillance. However, moving ships appear defocused in SAR images, which seriously affects the classification and identification of ships. Considering the three-dimensional (3-D) rotational motions (roll, pitch, and yaw) of the navigating ship, a novel method for refocusing moving ships in SAR images based on inverse synthetic aperture radar (ISAR) technique is proposed. First, a rectangular window is used to extract the defocused ship subimage. Next, the subimage is transformed into the ISAR equivalent echo domain, and the range migration and phase error caused by the identical movement of all ship scatterers are compensated. Then, the optimal imaging time can be selected by the maximum image contrast search method. Finally, the iterative adaptive approach (IAA) is used to obtain the image with high resolution. This method has satisfactory imaging performance in both azimuth resolution and image focus, and the amount of calculation is small due to the processing of subimages. Simulated data and Gaofen-3 real SAR data are used to verify the effectiveness of the proposed method.

1. Introduction

Synthetic aperture radar (SAR) plays an important role in the field of marine surveillance due to its all-weather and all-time imaging capability. Although stationary scenes can be well focused through traditional SAR imaging algorithms (e.g., Chirp Scaling, $\omega K$), moving ships can be defocused in the SAR image, which seriously affects ship classification and identification. Therefore, the research on the refocusing of the moving ship has become extremely urgent.

The effects of slow moving targets on SAR imaging are studied in [1], where the velocity and acceleration of the moving target are considered, and two airborne moving target indication (AMTI) methods are summarized. In [2], the fractional Fourier transform (FrFT) is used to perform optimal processing and analysis of residual chirp signals, thus the azimuth velocity of the target can be estimated and the along-track defocus can be compensated. In [3], the compensation function of the ground moving target in the two-dimensional (2-D) frequency domain is derived, and a method based on maximum image contrast search is proposed to obtain the refocused image of the ground moving target. On the basis of [3], sparse parameter representation method is proposed in [4] to refocus the moving target, which can suppress asymmetric side-lobes and improve the image quality. Based on the estimation of Doppler parameters (the Doppler centroid and the Doppler rate), the target motion parameters can be obtained and the ship inverse synthetic aperture radar (ISAR) imaging is carried out in [5]. An ISAR refocusing method for SAR ship image based on the fast minimum entropy phase compensation (FMEPC) method is proposed in [6].

However, the above-mentioned imaging methods are proposed for the smoothly moving target model, which are not suitable for the ship with three-dimensional (3-D) rotational motions (roll, pitch, and yaw). Because the Doppler frequency of the ship with 3-D rotational motions is time-varying after motion compensation, it is difficult for the above methods to obtain a focused image for such ships.

The effects of rotational motions (roll, pitch, and yaw) on ISAR imaging of small craft are analyzed in [7]. Moreover, in [8], the effects of ship linear motion and rotational motions on SAR image are researched, and the imaging distortions can be characterized by four main types of transformations, namely translation, rotation, scaling, and shearing according to the analysis and experiments. In [9], the effect of six-degree-of-freedom oscillation of ship target on SAR imaging is studied, the phase compensation method is made by fitting ship attitude angles with multiple sinusoidal functions, however, the ship attitude parameters are generally unknown in practice. In [10], multiple sinusoidal functions are used to match rotation effects of the ship in geosynchronous SAR, however, parameter searching needs a lot of calculation.

For the ISAR imaging with long recorded live ISAR data, a maximum image contrast based automatic time window selection (MC-ATWS) technology is proposed in [11], it can automatic select the position and length of the time window that provides the ISAR image with the highest image contrast, therefore, the relatively stable time of movement can be selected by this method. This method plays an effective role in ISAR imaging of non-cooperative targets, but there is a trade-off between ideal azimuth resolution and image focus. An ISAR motion compensation method via image contrast based technique is proposed in [12]. The applications of [11,12] are further expanded in [13], and a complete ISAR processing flow is proposed to deal with the defocused non-cooperative targets in Cosmo-Skymed spotlight SAR images. Moreover, the contrast maximization based technique is also verified using Radarsat-2 single look complex (SLC) SAR image in [14]. A novel ISAR optimal imaging time interval (OITI) selection algorithm for marine targets based on sea dynamic prior information (SDPI) is proposed in [15], the optimal imaging instant (OII) and optimal imaging duration (OID) are determined by the effective rotational vector (ERV) and maximum contrast criterion, respectively, however, precise SDPI acquisition and the feasibility for real data need to be further investigated. The time window selection algorithm based on measured data is proposed in [16], however, it is only suitable for cooperative ships. Phase nonlinearity is used to select imaging frames where target 2-D motion is dominant enough to obtain focused ISAR images in [17,18], however, those methods need prominent scatterers. The time window selection algorithm based on Doppler frequency information is proposed in [19], but it is not robust due to the instantaneous frequency variation of a ship is different for different scatterers. The time selection method based on the tracking of specific scatterers is proposed in [20], this method searches for the suitable side view imaging moment and top view imaging moment by tracking different ship specific scatterers, however, it is difficult to track the specific scatterers especially in the case of high sea conditions. The slope-based frame selection method proposed in [21] selects the time instants better suited for top or side view image formation, it operates in the image domain by using the slopes of some extracted linear features, however, it is difficult to extract the linear features. The time selection algorithm based on the Doppler center frequency (DCF) is proposed in [22], and the time selection algorithm based on the Doppler spread (DS) technique is proposed in [23], this type of method needs to obtain multiple range Doppler (RD) images, which is time consuming and not suitable for high sea conditions.

A ship ISAR imaging algorithm based on the generalized radon-Fourier transform (GRFT) is proposed in [24] through a 3-D search of range, velocity, and acceleration in the range compression domain, although this method adopts a combination of coarse search and fine search, it still consumes a lot of calculation, in addition, this algorithm only models each point as a combination of velocity and acceleration, ignoring higher order terms. GRFT and gradient-based descent algorithm are used for maneuvering targets imaging in [25], and this method takes a tradeoff between imaging performance and computational efficiency. Joint time-frequency transform is proposed in [26] as an alternative to Fourier transform to solve the image blurring problem caused by time-varying Doppler frequency shift, which does not need to resort to sophisticated motion compensation algorithms, however, joint time-frequency transform suffers from resolution degradation when the short-time Fourier transform (STFT) is used or cross-term interference when the Wigner-Ville distribution (WVD) is used. In [27], squint minimization operation is implemented for the compensation of third range compression in high-squint SAR ship imaging, then WASH-CLEAN method is proposed to automatically focus ship scattering points with different swings, however it is time consuming and it mainly considers the effect of ship motion under low or moderate sea states. Parametric range instantaneous Doppler algorithms express the Doppler phase of the ship scatterers as polynomial functions of observation time and the polynomial model coefficients can be estimated by product high-order matched-phase transform (PHMT) [28], coherently integrated generalized cubic phase function (CIGCPF) [29], coherently integrated cubic phase function (CICPF) [29], bilinear extended fractional Fourier transform (BEFRFT) [30], etc. The Doppler frequency histories of ship scatterers in the observation time can be reconstructed based on the estimated motion parameters, however, such methods need a lot of calculation, meanwhile, error propagation effects exist when estimating the parameters with such methods.

In the field of moving ship imaging, it is considered a very promising idea to process only the subimage which contains the defocused ship. This operation can bring many advantages. First, in the SAR image domain, the ship can be easily separated from the rest of the scene through many ship detection algorithms, in addition, processing subimage greatly saves computer resources, and the amount of calculation is significantly reduced. This paper considers the defocused navigating ship with 3-D rotational motions in the SAR image, the ISAR technique is used to image the ship, and the whole processing flow is based on the subimage, which contains the defocused ship. First, the detected ship is extracted through a rectangular window in the SAR image domain, next, the SAR subimage is transformed into ISAR equivalent echo domain through the azimuth inverse mapping (i.e., inverse fast Fourier transform (IFFT) in the azimuth direction) [31], and the motion compensation including range alignment and phase compensation is carried out, then, the optimal imaging time is selected through the maximum image contrast search method, finally, the iterative adaptive approach (IAA) [32] is used to obtain the focused image. The contradiction between time resolution and frequency resolution is alleviated, thus the proposed method has satisfactory imaging performance in both azimuth resolution and image focus for navigating ships with 3-D rotational motions.

Since the defocusing effect of the moving ship is a common phenomenon in SAR images, which severely limits the applications of SAR images, this study aims to propose a signal model of the navigating ship with 3-D rotational motions and refocus the moving ship with SAR subimage via ISAR technique, so as to provide a basis for ship classification and identification.

2. Geometry and Signal Model

The motion model used in this paper includes velocity and 3-D rotational motions including roll, pitch, and yaw. The geometric relationship between the radar and the ship target is shown in Figure 1a, the 3-D rotational motion model of the ship target is shown in Figure 1b, the velocity of the radar relative to the ground is v, the velocity of the ship relative to the ground is $v_s$, the angle between the orientation of the ship and the X axis is $\theta$.

For a ship scatterer $(x_0,y_0,z_0)$ in the ship coordinate system, its time-varying coordinates is $(x_1,y_1,z_1)$.

$$\begin{array}{c}\hfill \left[\begin{array}{c}{x}_1\\ y_1\\ z_1\end{array}\right]=M_y{M}_p{M}_r\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]\\ \hfill M_y=\left[\begin{array}{ccc}cos\left({\theta }_y\left(t\right)\right)& -sin\left({\theta }_y\left(t\right)\right)& 0\\ sin\left({\theta }_y\left(t\right)\right)& cos\left({\theta }_y\left(t\right)\right)& 0\\ 0& 0& 1\end{array}\right]\\ \hfill M_p=\left[\begin{array}{ccc}cos\left({\theta }_p\left(t\right)\right)& 0& sin\left({\theta }_p\left(t\right)\right)\\ 0& 1& 0\\ -sin\left({\theta }_p\left(t\right)\right)& 0& cos\left({\theta }_p\left(t\right)\right)\end{array}\right]\\ \hfill M_r=\left[\begin{array}{ccc}1& 0& 0\\ 0& cos\left({\theta }_r\left(t\right)\right)& -sin\left({\theta }_r\left(t\right)\right)\\ 0& sin\left({\theta }_r\left(t\right)\right)& cos\left({\theta }_r\left(t\right)\right)\end{array}\right]\end{array}$$

(1)

where t is the azimuth time, ${\theta }_r\left(t\right)$, ${\theta }_p\left(t\right)$, and ${\theta }_y\left(t\right)$ denote the angle of roll, pitch, and yaw, respectively, which can be expressed as

$$\begin{array}{cc}\hfill {\theta }_r\left(t\right)& =A_{r}sin(\frac{2\pi }{T_r}t+{\varphi }_r)\hfill \\ \hfill {\theta }_p\left(t\right)& =A_{p}sin(\frac{2\pi }{T_p}t+{\varphi }_p)\hfill \\ \hfill {\theta }_y\left(t\right)& =A_{y}sin(\frac{2\pi }{T_y}t+{\varphi }_y)\hfill \end{array}$$

(2)

where $A_r$, $T_r$, and ${\varphi }_r$ are the amplitude, period, and central time phase of the roll angle, $A_p$, $T_p$, and ${\varphi }_p$ are the amplitude, period, and central time phase of the pitch angle, $A_y$, $T_y$, and ${\varphi }_y$ are the amplitude, period, and central time phase of the yaw angle.

The coordinates of the scatterer in the global, earth-fixed coordinate system are $(x_2,y_2,z_2)$.

$$\left[\begin{array}{c}{x}_2\\ y_2\\ z_2\end{array}\right]=\left[\begin{array}{ccc}cos\left(\theta \right)& -sin\left(\theta \right)& 0\\ sin\left(\theta \right)& cos\left(\theta \right)& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_1+v_{s}t\\ y_1\\ z_1\end{array}\right]+\left[\begin{array}{c}\sqrt{R_0^2-H^2}\\ 0\\ 0\end{array}\right]$$

(3)

where $R_0$ is the distance between the radar and the ship at central time, H is the flight altitude of the radar.

Then, the distance between the scatterer $(x_2,y_2,z_2)$ and the radar can be written as

$$R\left(t\right)=\sqrt{{\left(x_2\right)}^2+{(y_2-vt)}^2+{(z_2-H)}^2}$$

(4)

As for ISAR imaging, the geometry is shown in Figure 2, the total angular rotation vector is ${\mathbf{\Omega }}_T\left(t\right)$, which is caused by the target rotation due to both the translational and angular motions, The projection of ${\mathbf{\Omega }}_T\left(t\right)$ on the plane orthogonal to the radar line of sight (LOS) is called effective rotation vector ${\mathbf{\Omega }}_{eff}\left(t\right)$, which contributes to the target aspect angle variation. The imaging plane is orthogonal to the effective rotation vector [13].

After range compression, the echo signal can be written as

$$S(t,f_r)=\sum _{i=1}^I{\sigma }_{i}rect\left(\frac{t}{T}\right)rect\left(\frac{f_r}{B}\right)exp(-j\frac{4\pi (f_0+f_r)R_i\left(t\right)}{c})$$

(5)

$$R_i\left(t\right)\approx R_0\left(t\right)+{\mathit{x}}^T{\mathit{i}}_{LOS}$$

(6)

where c is the velocity of light, I is the number of ship scatterers, ${\sigma }_i$ is a complex constant belonging to scatterer i, t is the azimuth time, T is the observation time, $f_r$ is the range frequency, B is the range bandwidth, $f_0$ is the carrier frequency, ${\mathit{i}}_{LOS}$ is the unit vector that identifies the LOS, $\mathit{x}$ is the vector that identifies a scatterer on the target. $R_0\left(t\right)$ represents the distance between the radar and the rotational center of the target, which is the main factor causing the target range migration and azimuth phase error, so $R_0\left(t\right)$ needs to be compensated before azimuth imaging, the operation to compensate the $R_0\left(t\right)$ named motion compensation.

Then, the Doppler frequency can be calculated from

$$f_d\left(t\right)=-\frac{2}{\lambda }\frac{\mathrm{d}{R}_i\left(t\right)}{\mathrm{d}t}$$

(7)

where $\lambda$ is the wavelength. After the phases introduced by the $R_0\left(t\right)$ of the rotational center are compensated, the Doppler frequencies caused by the roll, pitch, and yaw rotational motions are analyzed in this paper, considering the complex expression of 3-D rotational motions of the ship, the rotational motions are analyzed separately.

Assume the ith scatterer coordinates of the ship is $(x_i,y_i,z_i)$ in the ship coordinate system, its time-varying coordinates caused by roll is $(x_{it},y_{it},z_{it})$, then

$$\left[\begin{array}{c}{x}_{it}\\ y_{it}\\ z_{it}\end{array}\right]=M_r\left[\begin{array}{c}{x}_i\\ y_i\\ z_i\end{array}\right]$$

(8)

The Doppler frequency caused by roll is

$$f_{ri}\left(t\right)=-\frac{2}{\lambda }\frac{\mathrm{d}\left({\mathit{r}}_{it}^T{\mathit{i}}_{LOS}\right)}{\mathrm{d}t}$$

(9)

where ${\mathit{r}}_{it}={[x_{it},y_{it},z_{it}]}^T$, ${\mathit{i}}_{LOS}={[a_1,a_2,a_3]}^T$, then

$$f_{ri}\left(t\right)=-\frac{2{\omega }_r\left(t\right)}{\lambda }((y_{i}cos\left({\theta }_r\left(t\right)\right)-z_{i}sin\left({\theta }_r\left(t\right)\right))a_3-(y_{i}sin\left({\theta }_r\left(t\right)\right)+z_{i}cos\left({\theta }_r\left(t\right)\right))a_2)$$

(10)

where ${\omega }_r\left(t\right)=A_r\frac{2\pi }{T_r}cos(\frac{2\pi }{T_r}t+{\varphi }_r)$

Similarly, the Doppler frequencies caused by pitch and yaw are

$$\begin{array}{cc}\hfill f_{pi}\left(t\right)& =-\frac{2{\omega }_p\left(t\right)}{\lambda }((z_{i}cos\left({\theta }_p\left(t\right)\right)-x_{i}sin\left({\theta }_p\left(t\right)\right))a_1-(z_{i}sin\left({\theta }_p\left(t\right)\right)+x_{i}cos\left({\theta }_p\left(t\right)\right))a_3)\hfill \\ \hfill f_{yi}\left(t\right)& =-\frac{2{\omega }_y\left(t\right)}{\lambda }((x_{i}cos\left({\theta }_y\left(t\right)\right)-y_{i}sin\left({\theta }_y\left(t\right)\right))a_2-(x_{i}sin\left({\theta }_y\left(t\right)\right)+y_{i}cos\left({\theta }_y\left(t\right)\right))a_1)\hfill \end{array}$$

(11)

where ${\omega }_p\left(t\right)=A_p\frac{2\pi }{T_p}cos(\frac{2\pi }{T_p}t+{\varphi }_p)$ and ${\omega }_y\left(t\right)=A_y\frac{2\pi }{T_y}cos(\frac{2\pi }{T_y}t+{\varphi }_y)$, when the ${\theta }_r\left(t\right)$, ${\theta }_p\left(t\right)$, and ${\theta }_y\left(t\right)$ are small enough, then

$$\begin{array}{cc}\hfill sin\left({\theta }_r\left(t\right)\right)& \approx 0,sin\left({\theta }_p\left(t\right)\right)\approx 0,sin\left({\theta }_y\left(t\right)\right)\approx 0\hfill \\ \hfill cos\left({\theta }_r\left(t\right)\right)& \approx 1,cos\left({\theta }_p\left(t\right)\right)\approx 1,cos\left({\theta }_y\left(t\right)\right)\approx 1\hfill \end{array}$$

(12)

Combining the Doppler frequencies of roll, pitch, and yaw, then

$$\begin{array}{cc}\hfill f_{roti}\left(t\right)& =f_{ri}\left(t\right)+f_{pi}\left(t\right)+f_{yi}\left(t\right)\hfill \\ \hfill & \approx -\frac{2}{\lambda }((y_i{a}_3-z_i{a}_2){\omega }_r\left(t\right)+(z_i{a}_1-x_i{a}_3){\omega }_p\left(t\right)+(x_i{a}_2-y_i{a}_1){\omega }_y\left(t\right))\hfill \end{array}$$

(13)

Thus, the Doppler frequency variation introduced by 3-D rotational motions of each scatterer is the superposition of multiple sinusoidal signals, which allows us to find the optimal imaging time when the Doppler frequency is approximately constant.

3. Method

In this section, the proposed method for refocusing moving ships is described, it is based on ISAR technique with SAR subimage and the flowchart of the proposed method is shown in Figure 3.

3.1. Motion Compensation

The motion compensation should be done before azimuth imaging, which contains range alignment and phase compensation [6,31].

For the rotational center of the ship, its relative distance to the radar can be expressed as

$$R_0\left(t\right)\approx R_0+{\beta }_{1}t+{\beta }_2{t}^2$$

(14)

Since the azimuth velocity of the ship is always far smaller than the velocity of the radar (especially for the satellite platform), it can be neglected when analyzing the range migration (assume the range migration introduced by the azimuth velocity is compensated through SAR imaging). So the range migration of the moving ship relative to the radar is mainly caused by the slant range velocity, which is ${\beta }_1$, however, the phase introduced by ${\beta }_2$ cannot be neglected, which named Doppler rate mismatch in SAR imaging. So in the SAR SLC image, moving target will be defocused in the azimuth direction, and the energy will spread to different range cells at the same time.

The subimage can be transformed into azimuth time range frequency domain through inverse $\omega K$ algorithm, since the ship always occupies a small range cells, the center range of the subimage can be used as the reference range of the inverse $\omega K$ algorithm, and the inverse Stolt interpolation is neglected to decrease the computation time. Note that zero padding is required in this step. In order to achieve range alignment, the compensation function can be given by

$$H_c(t,f_r)=exp\left(j\frac{4\pi }{c}(f_r+f_0)\alpha t\right)$$

(15)

Then, the SAR image can be obtained by $\omega K$ algorithm and the Stolt interpolation is also neglected here, then, the entropy of the sum range profile (SRP) [33] defined in the ISAR equivalent echo domain can be used as the evaluation function of range alignment.

Then, $\alpha$ can be obtained through search, it corresponds to the value at which the entropy function is minimized. Finally, after coarse alignment, minimizing the entropy of the average range profile (ARP) [34] in the ISAR equivalent echo domain can be used for fine adjustment.

After the range alignment, the phase compensation should be carried out [6], which is an important step and second-order and higher-order phase errors should be compensated in this step, since they seriously affect the quality of focus. There are many phase compensation algorithms for SAR and ISAR, such as Doppler centroid tracking (DCT) [35], phase gradient autofocus (PGA) [36], fast minimum entropy phase compensation (FMEPC) [6,37], and maximum-contrast phase adjustment [37]. Since image entropy is an important indicator for the image overall focusing effect, the FMEPC algorithm is used in this paper to compensate the phase error. FMEPC estimates the phase error by minimizing image entropy iteratively and the specific principle is as follows.

Assume that the ISAR image is obtained by

$$g(k,n)=\sum _{u=0}^{U-1}y(u,n)exp\left(j\varphi \left(u\right)\right)exp(-j\frac{2\pi }{U}ku)$$

(16)

where U and N are the number of azimuth cells and range cells, respectively, $u=0,1,\dots ,U-1$ and $n=0,1,\dots ,N-1$, $y(u,n)$ denotes the equivalent echo after range alignment, k denotes the index of azimuth cells, $\varphi$ denotes the compensation phase, the entropy of the ISAR image [37] is

$$\begin{array}{cc}\hfill E& =\sum _{k=0}^{U-1}\sum _{n=0}^{N-1}\frac{{\left|g(k,n)\right|}^2}{S_g}ln\frac{S_g}{{\left|g(k,n)\right|}^2}\hfill \\ \hfill S_g& =\sum _{k=0}^{U-1}\sum _{n=0}^{N-1}{\left|g(k,n)\right|}^2\hfill \end{array}$$

(17)

Since $S_g$ is a constant, entropy can be redefined as [37]

$$E^{\prime }=-\sum _{k=0}^{U-1}\sum _{n=0}^{N-1}{\left|g(k,n)\right|}^{2}ln{\left|g(k,n)\right|}^2$$

(18)

To minimise $E^{\prime }$, the derivative of $E^{\prime }$ with respect to $\varphi \left(u\right)$ needs to be 0, which can be written as

$$\frac{\partial E^{\prime }}{\partial \varphi \left(u\right)}=0$$

(19)

According to [37], the derivative of $E^{\prime }$ with respect to $\varphi \left(u\right)$ can be obtained as

$$\begin{array}{cc}\hfill \frac{\partial E^{\prime }}{\partial \varphi \left(u\right)}& =2UIm\left\{exp\left(j\varphi \left(u\right)\right)a^{*}\left(u\right)\right\}\hfill \\ \hfill a\left(u\right)& =\sum _{n=0}^{N-1}{y}^{*}(u,n)\frac{1}{U}\sum _{k=0}^{U-1}{[1+ln|g(k,n)|}^2]g(k,n)exp\left(j\frac{2\pi }{U}ku\right)\hfill \end{array}$$

(20)

Thus, $\varphi \left(u\right)=angle\left(a\right(u\left)\right)$, and $angle(•)$ denotes phase calculation operation.

The main steps of FMEPC are as follows.

1.

Initialize the compensation phase $\varphi \left(u\right)$, $u=0,1,\dots ,U-1$;

2.

According to the Equation (16), The echo after range alignment $y(u,n)$ can be used to obtain ISAR image $g(k,n)$;

3.

Calculate the entropy of ISAR image, and set the termination condition (e.g., if the entropy difference between this iteration and the previous iteration is small, then, break; otherwise, go to next step);

4.

Calculate the value of $a\left(u\right)$ in this iteration process;

5.

Update the compensation phase $\varphi \left(u\right)$, $u=0,1,\dots ,U-1$, go to step 2.

Note that the compensation phase $\varphi \left(u\right)$ can be initialized through DCT method. After motion compensation, it can be considered that the translation of the ship is compensated, some ships with smooth motion can be focused well, nevertheless, the ship with 3-D rotational motions cannot be focused due to the time-varying Doppler frequency, which needs further processing, and the final focusing performance will depend on the ISAR azimuth imaging method.

3.2. Image Contrast Based Optimal Imaging Time Selection

After the motion compensation of the ISAR echo, it can be considered that the phase error introduced by the identical motion of the ship scatterers relative to the radar is compensated. The rest 3-D rotational motions of the various ship scatterers around the center of rotation will cause the time-varying Doppler, and the time-varying Doppler will make the ship defocused. However, in some short time intervals, the Doppler frequency variation of the ship scatterers is nearly stable. It is completely possible to form the image by selecting the time interval during which the ship is moving smoothly based on maximum image contrast search method, and the rotation vector is assumed to be approximately constant during the selected imaging time. However, the azimuth resolution will decrease when the time window is used with RD imaging. This paper selects the time interval of smooth motion for imaging with a short constant length. The image contrast [38] can be defined as

$$C=\frac{std\left(\right|I(m,n){|}^2)}{mean\left(\right|I(m,n){|}^2)}$$

(21)

where $std(•)$ represents the standard deviation operation, $mean(•)$ represents the average operation, m and n denote the index of azimuth cells and range cells, respectively, $I(m,n)$ denotes the image obtained by RD imaging.

3.3. IAA

After time selection, the azimuth resolution of ISAR echo is decreased, which is harmful for subsequent classification and identification. In order to optimize the imaging results, the method of IAA is used for azimuth imaging after the optimal imaging time selection, IAA is a nonparametric, weighted least squares (WLS)–based spectral estimation algorithm that can increase the resolution and reduce sidelobes. It has been used in source localization and sensing [39], multiple-input multiple-output (MIMO) radar imaging [40], super-resolution Doppler beam sharpening imaging [32], etc. In this subsection, the IAA is used to achieve azimuth super-resolution imaging after optimal imaging time selection.

Assume that ${\mathit{y}}_n={[y_n\left(1\right),y_n\left(2\right),\dots ,y_n\left(M\right)]}^T$ and ${\mathit{e}}_n={[e_n\left(1\right),e_n\left(2\right),\dots ,e_n\left(M\right)]}^T$ denote the azimuth sample and azimuth noise of the nth range cell after optimal imaging time selection, ${\mathit{s}}_n={[s_n\left(1\right),s_n\left(2\right),\dots ,s_n\left(K\right)]}^T$ denotes the reflectivity of the scatterers at each single frequency in the nth range cell. Then ${\mathit{y}}_n$ can be expressed as

$${\mathit{y}}_n=\mathit{A}{\mathit{s}}_n+{\mathit{e}}_n$$

(22)

where $\mathit{A}$ is the redundant time–frequency dictionary and $\mathit{A}=[{\mathit{a}}_1,{\mathit{a}}_2,\dots ,{\mathit{a}}_K]$, ${\mathit{a}}_k={[1,exp\left(j2\pi f_k\right),\dots ,exp\left(j2\pi f_k(M-1)\right)]}^T$, $f_k=\frac{k-1}{K}-\frac{1}{2}$. Then ${\mathit{s}}_n$ can be obtained using the weighted least squares method [32], given by

$$s_n\left(k\right)=argmin\left|\right|{\mathit{y}}_n-{\mathit{a}}_k{s}_n\left(k\right){\left|\right|}_{{\mathit{Q}}_k^{-1}}^2\triangleq argmin{({\mathit{y}}_n-{\mathit{a}}_k{s}_n\left(k\right))}^H{\mathit{Q}}_k^{-1}({\mathit{y}}_n-{\mathit{a}}_k{s}_n\left(k\right))$$

(23)

where ≜ means “be defined as”, ${\mathit{Q}}_k$ denotes the interference covariance matrix which is composed of all the Doppler components except the Doppler component corresponding to $f_k$, given by

$${\mathit{Q}}_k=\mathit{R}-p_n\left(k\right){\mathit{a}}_k{\mathit{a}}_k^H$$

(24)

where $\mathit{R}$ is the covariance matrix and $\mathit{R}={\sum }_{k=1}^K{p}_n\left(k\right){\mathit{a}}_k{\mathit{a}}_k^H$, $p_n\left(k\right)$ is the power of the target located at $f_k$ and $p_n\left(k\right)={\left|s_n\left(k\right)\right|}^2$.

Assume that $\mathit{P}=diag\left([p_n\left(1\right),p_n\left(2\right),\dots ,p_n\left(K\right)]\right)$, then $\mathit{R}=\mathit{A}\mathit{P}{\mathit{A}}^H$, where $diag(•)$ means the transformation from vector to diagonal matrix. We can get the weighted least squares estimation of $s_n\left(k\right)$ in an iterative way [32], given by

$${\hat{s}}_n\left(k\right)=\frac{{\mathit{a}}_k^H{\mathit{Q}}_k^{-1}{\mathit{y}}_n}{{\mathit{a}}_k^H{\mathit{Q}}_k^{-1}{\mathit{a}}_k}$$

(25)

where $\widehat{}$ means the estimated value of a variable.

According to the matrix inversion lemma, the ${\mathit{Q}}_k$ in Equation (25) can be replaced by $\mathit{R}$, thus the estimation of $s_n\left(k\right)$ [32] can be written as

$${\widehat{s}}_n\left(k\right)=\frac{{\mathit{a}}_k^H{\mathit{R}}^{-1}{\mathit{y}}_n}{{\mathit{a}}_k^H{\mathit{R}}^{-1}{\mathit{a}}_k}$$

(26)

$\mathit{R}$ is initialized to be an identity matrix, and Algorithm 1 illustrates the specific implementation of IAA. Note that for the case that $\mathit{R}$ is ill conditioned, the problem may be handled by regularized IAA with minor modification to the IAA [40].

Algorithm 1. The specific implementation of IAA.

Initialization:   ${\widehat{p}}_n\left(k\right)=\frac{|{\mathit{a}}_k^H{\mathit{y}}_n{|}^2}{{\left({\mathit{a}}_k^H{\mathit{a}}_k\right)}^2},k=1,2,\dots K$;

Repeat:

$\hat{\mathit{P}}=diag\left([{\widehat{p}}_n\left(1\right),{\widehat{p}}_n\left(2\right),\dots ,{\widehat{p}}_n\left(K\right)]\right)$;

$\mathit{R}=\mathit{A}\widehat{\mathit{P}}{\mathit{A}}^H$;

for $k=1,2,\dots K$

${\widehat{s}}_n\left(k\right)=\frac{{\mathit{a}}_k^H{\mathit{R}}^{-1}{\mathit{y}}_n}{{\mathit{a}}_k^H{\mathit{R}}^{-1}{\mathit{a}}_k};$

${\widehat{p}}_n\left(k\right)={\left|{\widehat{s}}_n\left(k\right)\right|}^2$;

end for

Until a certain number of iterations is reached.

The number of iterations for IAA is set to 15 in this paper according to [32], which will make IAA very efficient and fast. Perform IAA at each range cell, and the ISAR image will be obtained with better image quality than Fourier transform.

4. Experimental Results

4.1. IAA Super-Resolution Experiment

Since the time selection in the azimuth direction leads to a reduction in azimuth resolution, the spectrum estimation is only carried out in the azimuth direction range cell by range cell, which allows us to simulate a one-dimensional (1-D) signal to illustrate the problem. Assuming that the sampling rate of the signal is 1000 Hz and the number of samples is 1000, thus the frequency resolution is 1 Hz. Six typical targets numbered 1–6 are considered in the simulation, and the corresponding frequencies are $-100$, $-98$, $-31$, $-20$, 21, 30, the corresponding amplitudes are 1, 1, 1, 0.4, 0.2, 0.2. Different targets are for different purposes, target 1 and target 2 represent targets with the same amplitude and adjacent frequency, target 3 and target 4 represent targets with adjacent frequency and large amplitude difference, target 5 and target 6 represent targets with low power, the signal noise obeys complex Gaussian distribution, and the signal-to-noise ratio (SNR) is set to 20 dB. Because of the high frequency resolution, the fast Fourier transform (FFT) can be used for spectrum estimation as shown in Figure 4a. Then the signal is truncated, and only 250 samples in the middle of the signal are selected, thus the frequency resolution is 4 Hz by theoretical calculation. As shown in Figure 4b, the frequencies of target 1 and target 2 are difficult to distinguish by FFT, meanwhile, the amplitude estimation error is large. With the zero padding FFT, the target 1 and target 2 are still difficult to distinguish in Figure 4c. IAA is used for spectrum estimation in Figure 4d, the frequencies of target 1 and target 2 are separated, and the amplitude estimation error is acceptable.

4.2. Moving Ship Simulation Experiment

In order to verify the focusing performance of the proposed method on the defocused ship in the SAR image, two simulation experiments are carried out in the subsection. The simulation model is shown in Figure 5 and the system parameters of the radar platform are shown in

Table 1. The system parameters of the simulation experiment.

Parameter	Value
Carrier frequency (GHz)	5.4
Platform velocity (m/s)	150
Slant range of scene center (Km)	10
A/D sampling rate (MHz)	240
Squint angle ()	0
Bandwidth (MHz)	200
Pulse width ($\mathsf{\mu }$s)	1
Pulse repetition frequency (Hz)	750

.

In the first simulation experiment, a navigating ship with 2-D (range and azimuth) velocities is studied. The angle between the ship orientation and the X axis is 45, the velocity of the ship is set to 5 m/s, and the 3-D (roll, pitch, and yaw) rotational velocities of the ship are set to 0.

Figure 6a shows the defocusing phenomenon caused by the sailing velocity of the ship. Then, SAR subimage is first transformed into the ISAR echo domain through the inverse mapping, and the motion compensation method is carried out. Figure 6b shows the RD imaging result after motion compensation. The ship target is focused after motion compensation, and the image quality is improved, which indicates the effectiveness of the motion compensation method.

However, some ships on the sea will produce 3-D rotational motions due to the sea waves and other factors. In this case, the Doppler frequency of the scatterer on the ship after motion compensation is not constant, and it varies with the change of azimuth time. Since each scatterer on the ship has different time-varying Doppler frequency, it is difficult to image all the ship scatterers.

After motion compensation, the echo data in a short imaging time is selected for ISAR imaging based on the maximum image contrast search method. Since the Doppler frequency in the selected imaging time is assumed to be approximately constant, the IAA is used to achieve ISAR super-resolution imaging.

In the second simulation experiment, a navigating ship with the same velocity and orientation as the first simulation experiment is studied, meanwhile 3-D rotational motions are added, and the parameters of rotational motions are shown in

Table 2. The parameters of rotational motions.

Rotational Mode	Amplitude (°)	Period (s)	Central Time Phase (°)
Roll	5	12.2	0
Pitch	1.7	6.7	0
Yaw	1.9	14.2	0

.

Figure 7a shows a SAR subimage that contains a defocused ship. Due to the translational motion and 3-D rotational motions of the ship, the SAR image of the ship is severely defocused, and the scatterers on the ship are difficult to identify. The image after motion compensation is shown in Figure 7b, the imaging performance is improved, but the ship scatterers still have severe defocusing phenomenon after motion compensation. The reason for this phenomenon is that even if translational motion compensation is performed, the Doppler frequency of the scatterer caused by the 3-D rotational motions is still changing. After optimal imaging time selection, the Doppler frequency of the ship is approximately constant within the selected imaging time and the RD imaging result is shown in Figure 7c, the azimuth resolution is low. IAA imaging result is shown in Figure 7d, the image quality is improved, and the azimuth resolution is improved.

4.3. Real Data of Gaofen-3

Because the stripmap mode of Gaofen-3 has short illuminated time and low resolution, in this subsection, Gaofen-3 sliding spotlight SAR data are used to verify the effectiveness of the proposed method. Some system parameters of Gaofen-3 sliding spotlight mode are shown in

Table 3. The system parameters of Gaofen-3.

Parameter	Value
Carrier frequency (GHz)	5.4
Platform velocity (m/s)	7571
Slant range of scene center (Km)	856.33
A/D sampling rate (MHz)	266.67
Bandwidth (MHz)	240
Pulse width ($\mathsf{\mu }$s)	45
Pulse repetition frequency (Hz)	3591

. The resolution of this mode can reach to 1 m [41], and the defocusing phenomenon of the moving ship is apparent in the SAR image, which seriously affects the applications of SAR image, such as ship classification and identification. Because the general SAR processing flow for the sliding spotlight mode of Gaofen-3 is complicated, and the energy spread across the range cells is not obvious for Gaofen-3 data. Therefore, this paper directly performs range alignment in the ISAR equivalent echo domain by minimizing the entropy of the ARP [34].

In the first SAR subimage, a navigating ship is used to verify the motion compensation method. The 3-D rotational motions of the ship can be neglected in this subimage. The SAR ship subimage is shown in Figure 8a, the ship appears defocused in the subimage, which is mainly caused by the 2-D sailing velocities. The motion compensation result is shown in Figure 8b, the ship is well focused after motion compensation through RD imaging method, this result verifies the effectiveness of the motion compensation method.

In order to verify the effectiveness of the proposed method, two navigating ships with 3-D rotational motions named rotating ship 1 and rotating ship 2 are used. For the rotating ship 1, the defocused SAR subimage is shown in Figure 9a. The motion compensation result is shown in Figure 9b, the ship still appears defocused due to the rotational motions introduced by sea waves and other factors, and the bow and stern cannot be focused at the same time. Figure 9c shows the ship image after optimal imaging time selection through the maximum image contrast search method, the resolution is reduced after time selection, which is the same as the STFT in [26], however, the Doppler frequency is almost constant in the imaging time period, which allows us to obtain ship image through super-resolution spectral estimation. Finally, IAA is used to obtain the ship image range cell by range cell, and the final ship image is shown in Figure 9d.

As shown in Figure 10a, the rotating ship 2 suffers from image distortion due to the velocity and 3-D rotational motions, the processing flow is the same as the rotating ship 1, the image after motion compensation is shown in Figure 10b, the ship still appears defocused, especially the stern, the optimal imaging time selection result is shown in Figure 10c, the result of IAA is shown in Figure 10d, the image quality is improved after the processing, which shows the effectiveness of the proposed imaging method.

5. Discussion

In this paper, the refocusing method of the navigating ship with 3-D rotational motions is studied. Due to the advantages of ISAR technique in non-cooperative target imaging, the ISAR equivalent echo is obtained from SAR subimage which contains a single defocused ship. Then motion compensation including range alignment and phase compensation is carried out to compensate the identical movement of all ship scatterers, the moving ship whose 3-D rotational motions can be neglected will form a well focused image through RD imaging method, as shown in Figure 6, the simulation ship is well focused after motion compensation via RD imaging, moreover, the motion compensation performance is also verified by the real SAR data of the Gaofen-3 in Figure 8. Both simulated data and real SAR data have a good focus effect and high resolution.

However, many ships with 3-D rotational motions will have time-varying Doppler frequency after motion compensation, which will affect the RD imaging result. For simulation ship in Figure 7, the ship still suffers from defocusing after motion compensation, a short imaging time is selected by the maximum image contrast search method, and the RD imaging result is shown in Figure 7c, the motion during the selected imaging time is considered as smooth, thus the Doppler frequency variation is not obvious, the time selection method has the similar principle as STFT, and shorter imaging time will bring worse azimuth resolution since the look angle variation is small. To improve the resolution, the IAA is used to achieve the azimuth imaging, as shown in Figure 7d, the azimuth resolution has been improved after IAA imaging. The contradiction between time resolution and frequency resolution is alleviated, thus the proposed method has satisfactory imaging performance in both azimuth resolution and image focus. The high-resolution SAR data of Gaofen-3 is also used to verify the proposed method in Figure 9 and Figure 10, the proposed method shows good imaging performance for ships with 3-D rotational motions. In the future, optimal imaging time selection algorithm will be further researched to select the smooth motion time accurately and fast with ISAR equivalent echo.

6. Conclusions

A novel method for refocusing moving ships in SAR images is proposed in this paper. Both the 2-D velocities and 3-D rotational motions of the moving ship are considered. Through the motion compensation in the ISAR equivalent echo domain obtained from the SAR subimage, the identical motion of all the scatterers on the ship can be compensated. Furthermore, the optimal imaging time is selected through the maximum image contrast search method. Finally the image is obtained through IAA with high resolution and low side lobes. The proposed method has good focusing performance for moving ships in both azimuth resolution and image focus, and the whole processing flow is based on SAR subimages, which has a small amount of calculation and great convenience.