Squinted Airborne Synthetic Aperture Radar Imaging with Unknown Curved Trajectory

Abstract

The imagery of airborne highly squinted synthetic aperture radar (SAR) with curved trajectory is a challenging task due to the translational-variant range cell migration (RCM) and azimuth modulation. However, in most cases of practical application, the curved trajectory cannot be accurately known, which brings greater difficulties to the imaging problem. To accommodate these issues, we propose a novel motion modelling and optimisation based imaging algorithm for the highly squinted SAR with unknown curved trajectory. First, to correct the translational-variant RCM, a coarse-to-fine RCM correction scheme as well as a range perturbation approach is applied. Afterwards, an optimisation model of motion information under the criterion of minimum entropy is built during the azimuth processing by nonlinear chirp scaling (NLCS). Correspondingly, a differential evolution (DE) optimisation strategy is proposed to estimate the motion information in an iterative manner. We empirically compare the proposed algorithms with several state-of-the-art highly squinted curved SAR imaging algorithms. Numerical results show the effectiveness of the proposed method in the case without any prior information of the curved trajectory.

1. Introduction

Synthetic aperture radar (SAR) [1,2,3] finds wide applications in both civilian and military domains due to its capability of operating under any weather conditions and any time of the day. SAR is an active radar with its own illumination to image a ground scene of interest with high resolution. Many types of improved SAR systems have emerged in recent years [4,5,6]. According to the angle between the radiation direction of the radar and movement direction, SAR systems can be divided into three categories: (1) side-looking SAR, (2) squint-looking SAR and (3) forward-looking or back-looking SAR, as shown in Figure 1. Among them, due to its manoeuvrability and flexibility, squinted SAR is a research hotpot [7,8,9,10]. In a conventional squinted SAR system, the platform moves along a straight trajectory. However, in many new applications, more and more SAR systems are mounted on the high-speed and manoeuvrable platform—such as unmanned aerial vehicles, air fighters, missile, satellite and so on—and the radar platform moves along a curved trajectory [11,12,13,14]. In this paper, we focus on the airborne squinted SAR system.

Squinted SAR imaging with known curved trajectory—where we have the the prior information of the moving curved trajectory of the radar platform—has been studied in various literature, and can be categorised into two groups:

• Time-domain imaging algorithm: The time-domain imaging algorithm—also known as Back Projection (BP) algorithm—is viewed as the ideal solution to the focusing problem; however, it is also time-consuming. There are many modifications to BP, and among them, the most popular one is fast factorised BP (FFBP) [15,16,17,18,19,20]. In a recursive manner, the computation cost is greatly reduced by FFBP and also could be further accelerated by parallel processing. Another advantage of time-domain imaging algorithms is that they are adaptable, and there is no need for significant modification to applying conventional time-domain imaging algorithms—where platform moves along a linear trajectory—to the squinted SAR imaging with curved trajectory [21,22].
• Frequency-domain imaging algorithm: Some frequency domain algorithms—including chirp scaling (CS) [23], nonlinear chirp scaling (NLCS) [24,25], Omega-K and their extensions [12,26]—are proposed to achieve well-focused image for squinted SAR with curved trajectory. These frequency-domain algorithms are well-designed to cope with the transnational variant range cell migration (RCM) and azimuth phase in a perturbation, equalisation or interpolation manner.

However, in practical application, there is always some motion uncertainty—which means the deviations of the actual trajectory from the pre-designed one—for the moving platform, mainly caused by air turbulence in airborne SAR system. As a consequence, the focusing quality of the existing imaging algorithms mentioned above is strongly impaired. Principally speaking, we can get the motion information of the moving platform through orbit and altitude data provided by an ancillary instrument such as inertial measurement units (IMU), inertial navigation system (INS) and global positioning system (GPS) [27,28,29,30,31]. However, unfortunately, measurement uncertainties on the IMU, INS, GPS and so on would limit the accuracy of the moving information, and the measurement remains unknown in many real applications including some unmanned aerial vehicles without an ancillary instrument. Thus, imaging techniques of squinted SAR without any prior information of the curved moving trajectory are of great value.

Note that this problem is similar with the autofocusing problem [32,33,34] in SAR motion compensation, as both of them focus on eliminating the uncertainty of the platform’s motion that deviates from the expected trajectory. However, these two types of problems are not completely consistent, i.e., they and their corresponding solutions differ in the aims and constraints. The differences lie in the following factors.

• Purpose: The imaging and autofocusing are an organic combination to achieve high-precision imaging together. The imaging algorithms aim to focusing SAR echo under several geometry parameters—motion parameters like initial position, speed and acceleration of the platform. Comparatively, the purpose of the autofocusing algorithms is to compensate for the motion errors—caused by unexpected deviation from the geometry configuration—which are not taken into consideration in the imaging algorithms.
• Constraint: The imaging algorithms cannot deal with the high frequency motion uncertainty as the they are parametrically designed for the geometry parameters. The constraint for the autofocusing algorithms is that they cannot compensate for the motion uncertainty with large amplitude. The corresponding reason is that the complex motion errors—that are always with high-order terms or even non-parameterisable—lead to many assumptions including approximation on residual range cell migration (RCM), neglecting translational variance on azimuth modulation.

Consequently, the existing autofocusing methods cannot solve the imaging problem without any prior information of the curved trajectory.

To accommodate for this problem, a novel squinted SAR imaging algorithm without any prior information of the curved trajectory is proposed. The main contributions are as follows.

• To correct the translational variant RCM, a coarse-to-fine RCM correction scheme integrated with a range perturbation approach is proposed. The coarse-to-fine RCM correction scheme works in an iterative manner, which could guarantee the prerequisite of the following processors on azimuth modulation. A range perturbation approach is utilised to correct the translational variant RCM based on the estimated motion parameters.
• We establish an optimisation model for the motion parameters under the minimum entropy criterion based on NLCS processing. At this stage, the translational variance of the azimuth phase modulation has been taken into consideration and the problem is changed into a multi-variable minimisation problem.
• The minimisation problem is solved by a differential evolution (DE) strategy. The estimated motion parameters would be utilised to guide the RCM correction in the next iteration, which results in a more accurate estimation iteratively.
• We conduct experiments using synthetic SAR data to show the effectiveness of the proposed method.

The rest of this paper is organised as follows. Section 2 introduces squinted SAR with curved trajectory. In Section 3, a motion modelling and optimisation-based imaging algorithm is derived in detail. We make some comments on the proposed method in Section 4, and numerical results are given in Section 5. Section 6 concludes this paper.

2. Squinted Sar with Curved Trajectory

Assuming that the transmitted signal is $p\left(\tau \right)$, the SAR echo corresponding to arbitrary target $P(x,y)$ the imaging scene could be formulated as

$$\begin{array}{c}\hfill s\left(\tau ,t\right)={\omega }_a\left(t\right)exp\left(-j4\pi f_c\frac{r\left(t;x,y\right)}{c}\right)p\left(\tau -\frac{2r\left(t;x,y\right)}{c}\right),\end{array}$$

(1)

where $\tau$ and t denote the fast and slow time, respectively, while x and y represent the azimuth and range directions with respect to imaging scene $\sigma$, respectively. In (1), ${\omega }_a$ is the azimuth envelope, $r(x,y,n)$ the range history, $f_c$ the carrier frequency and c the signal speed. After range compression, SAR echo is changed into

$$\begin{array}{c}\hfill s_r\left(\tau ,t\right)={\omega }_a\left(t\right)exp\left(-j4\pi f_c\frac{r\left(t;x,y\right)}{c}\right)\mathrm{sinc}\left(\tau -\frac{2r\left(t;x,y\right)}{c}\right),\end{array}$$

(2)

where $\mathrm{sinc}\left(\right)$ represents the range envelope after range compression.

Figure 2 presents the geometry configuration of the squinted airborne SAR with motion errors. Assume that the radar platform moves along the direction of y axis with a constant velocity v as well as accelerations $a_x$ and $a_z$ along x- and z-axes, and the cases with more complex motions will be discussed in Section 4. Note that higher-order motion along y-axis is originated from along-track nominal velocity changes, and it is generally compensated via azimuth resampling of SAR raw data [7,35]. Consequently, higher-order motion along y-axis is not taken into consideration in this paper. x is the ground direction orthogonal to y. We choose the centre point of the imaging scene, i.e., target O in Figure 2, as the coordinate origin. t denotes the azimuth time, and it is chosen to be zero at the composite beam centre crossing time of target O.

Suppose target $P(x,y)$ is an arbitrary target located at the imaging area, and the corresponding range history of P can be formulated as

$$\begin{array}{c}\hfill r(t;x,y)=\sqrt{{\left(x_0+a_x{t}^2-x\right)}^2+{\left(y_0+vt-y\right)}^2+{\left(z_0+a_y{t}^2\right)}^2},\end{array}$$

(3)

where $x_0$, $y_0$ and $z_0$ denote the initial position coordinate of the radar platform with respect to x, y and z directions, respectively. $a_x$ and $a_z$ denote the accelerations along x and z directions, respectively.

In order to facilitate subsequent imaging processing, a third-order polynomial range model is utilised to accurately characterise $r(t;x,y)$ as

$$\begin{array}{c}\hfill r(t;x,y)=k_0+k_1(t-\frac{y}{v})+k_2{(t-\frac{y}{v})}^2+k_3{(t-\frac{y}{v})}^3+\dots \end{array}$$

(4)

where

$$\begin{array}{c}\hfill k_0=\sqrt{{\left(x_0+\frac{a_x{y}^2}{v^2}-x\right)}^2+y_0^2+{\left(z_0+\frac{a_x{y}^2}{v^2}\right)}^2},\end{array}$$

(5)

$$\begin{array}{cc}\hfill k_1& ={\left.\frac{\partial r(t;x,y)}{\partial t}\right|}_{t=y/v}\hfill \\ \hfill & =\frac{\left(4a_{x}y(-x+x_0+\frac{a_x{y}^2}{v^2})\right)/v+2v{y}_0+\left(4a_{z}y(\frac{a_z{y}^2}{v^2}+z_0)\right)/v}{2{({(-x+x_0+\frac{a_x{y}^2}{v^2})}^2+y_0^2+{(\frac{a_z{y}^2}{v^2}+z_0)}^2)}^{0.5}},\hfill \end{array}$$

(6)

and

$$\begin{array}{cc}\hfill k_2=& {\left.\frac{{\partial }^{2}r(t;x,y)}{2\partial t^2}\right|}_{t=y/v}\hfill \\ \hfill =& \frac{v^4+(6a_x^2+6a_z^2)y^2+v^2(-2a_{x}x+2a_x{x}_0+2a_z{z}_0)}{2v^2{({(-x+x_0+\frac{a_x{y}^2}{v^2})}^2+y_0^2+{(\frac{a_z{y}^2}{v^2}+z_0)}^2)}^{0.5}}\hfill \\ \hfill & -\frac{{(\left(2a_{x}y(-x+x_0+\frac{a_x{y}^2}{v^2})\right)/v+v{y}_0+\left(2a_{z}y(a_z{y}^2+v^2{z}_0)\right)/v^3)}^2}{4{({(-x+x_0+\frac{a_x{y}^2}{v^2})}^2+y_0^2+{(\frac{a_z{y}^2}{v^2}+z_0)}^2)}^{1.5}}\hfill \end{array}$$

(7)

In (4), $k_3$ is the third-order term of the polynomial; however, we do not show the concrete expression of $k_3$ which would not be used in the following.

3. Motion Modelling and Optimisation

Many literatures are focused on the curved trajectory squinted SAR imaging problem; however, all of them need to know the accurate information of motion parameters. Once there is a slight error occurring in these parameters, the corresponding imaging quality is greatly reduced, which will be illustrated in the Section 5. In this paper, our goal is to develop a imaging algorithm not needing the prior information of the motion parameters, which could achieve a well-focused SAR image as well as estimate the motion parameters at the same time.

In order to solve the squinted SAR imaging problem with unknown curved trajectory, there are three issues to be solved: (1) translational variant RCM correction, (2) translational variant azimuth modulation and (3) unknown motion parameters. Two solutions are proposed to solve these issues: (1) equalise and correct the translational variant RCM by range perturbation processing, and (2) model and estimate the motion parameters through a NLCS processing with azimuth variance of the Doppler parameters taken into consideration, which solve the last two issues at the same time. However, there is a contradiction between the two solutions—the RCM should be corrected to ensure an accurate NLCS-based motion estimation, while an accurate RCM correction need the estimated motion parameters. Here, we apply a coarse-to-fine RCM correction scheme to address this contradiction.

3.1. Coarse-to-Fine Rcm Correction Scheme

The coarse-to-fine RCM correction scheme begins with the observation that if the the residual RCM after RCM correction is smaller than one range resolution, the motion estimation can be implemented correctly. The correction scheme works in an iterative manner.

At first, we use the autofocus RCM correction algorithm in [36] to compensate for the translational-invariant component of the RCM, and the translational variant part remains to be corrected as it is larger than one range resolution in most cases. Then, the echo is changed into a coarser range resolution by adjusting the bandwidth so that the uncompensated RCM is within one range resolution. Once done, the motion parameters can be accurately estimated and azimuth focusing can be implemented from the SAR data with coarser resolution through a NLCS processing. Using estimated and updated motion parameters, the RCM is compensated by a range perturbation processing with translational variance taken into consideration. Then, we adjust range resolution into a finer one, and implement a more accurate motion parameter estimation. The block diagram of the coarse-to-fine RCM correction scheme is shown in Figure 3.

Note that the RCM correction in the first iteration is implemented using the RCM estimation and correction algorithm in [36], and that in the following iterations is performed by range perturbation processing.

The proposed algorithm stops until the change rate of the estimation of motion parameters is smaller than a predetermined threshold. The change rate in the $k+1$th iteration is defined as

$$\begin{array}{c}\hfill c^{k+1}=\frac{|a_x^{k+1}-a_x^k|+|a_z^{k+1}-a_z^k|}{|a_x^k|+|a_z^k|},\end{array}$$

(8)

where $a_x^k$ and $a_z^k$ represent the estimation of $a_x$ and $a_z$ at the kth iteration.

3.2. RCM Correction

In the first iteration, we do not have any estimation on the motion parameters $a_x$ and $a_z$, and therefore the RCM is corrected in an autofocus manner [36,37]. In this algorithm, the RCM of the most dominant target can be estimated and compensated without needing the motion information; however, the translational-variance of RCM is neglected.

In the following iterations, the RCM will be compensated more accurately based on the estimation of $a_x$ and $a_z$ from the last step “Motion estimation and azimuth focusing” shown in Figure 3. In the squinted SAR, the range migration is dominated by linear range walk, and a linear range cell migration correction (LRCMC) should be employed first. However, it introduces the azimuth-variance with respect to both RCM and azimuth modulation. Here, a translational-variant RCM correction processing based on range perturbation [23] is utilised as follows.

• Multiplying LRCMC factor $H_{rcm1}(t,f_{\tau })$ with signal $s_1(t,f_{\tau })$ in the range frequency and azimuth time domain to correct the linear range walk, where $H_{rcm1}(t,f_{\tau })$ is

  $$\begin{array}{c}\hfill H_{rcm1}(t,f_{\tau })=exp\left(j\frac{4\pi (f_c+f_{\tau })t}{c}{k}_{10}\right),\end{array}$$

  (9)
  $f_{\tau }$ denotes the range frequency, $s_1(t,f_{\tau })$ is obtained by applying Fourier transform in $\tau$ dimension for data $s_1(t,\tau )$ in (2) and

  $$\begin{array}{c}\hfill k_{10}=k_1{\mid }_{x=0,y=0}=\frac{v{y}_0}{{(x_0^2+y_0^2+z_0^2)}^{0.5}}.\end{array}$$

  (10)
• Multiplying range perturbation factor $H_1(f_t,f_{\tau })$ with signal obtained by step 1 in 2-D frequency domain to equalise the translational variance of range curvature introduced by LRCMC, where

  $$\begin{array}{c}\hfill H_1(f_t,f_{\tau })=exp\left(j\pi \frac{c^2{f}_{\tau }{f}_t^3}{2{(f_c+f_{\tau })}^3}(p_1+\frac{p_{2}c}{2(f_c+f_{\tau })}{f}_t)\right).\end{array}$$

  (11)
  In (11), $f_t$ denotes the azimuth frequency, and $p_1$ and $p_2$ are

  $$\begin{array}{c}\hfill p_1=-\frac{k_{11}^2{k}_{20}^2+4k_{10}{k}_{11}{k}_{20}{k}_{21}+k_{10}^2{k}_{21}^2}{12k_{10}{k}_{20}^7{k}_{11}^4(3k_{10}{k}_{21}+2k_{11}{k}_{20})},\end{array}$$

  (12)

  $$\begin{array}{c}\hfill p_2=-\frac{2k_{11}{k}_{20}{k}_{21}+k_{10}{k}_{21}^2}{12k_{10}{k}_{20}^7{k}_{11}^4(3k_{10}{k}_{21}+2k_{11}{k}_{20})},\end{array}$$

  (13)

  where

  $$\begin{array}{c}\hfill k_{11}={\left.\frac{\partial k_1}{\partial y}\right|}_{x=0,y=0}=\frac{2(a_x{x}_0+a_z{z}_0)}{v\sqrt{x_0^2+y_0^2+z_0^2}},\end{array}$$

  (14)

  $$\begin{array}{cc}\hfill k_{20}=k_2{|}_{x=0,y=0}=& \frac{v^2+2a_x{x}_0+2a_z{z}_0}{\sqrt{x_0^2+y_0^2+z_0^2}}\hfill \\ \hfill & -\frac{v^2{y}_0^2}{{(x_0^2+y_0^2+z_0^2)}^{1.5}},\hfill \end{array}$$

  (15)

  and

  $$\begin{array}{c}\hfill k_{21}={\left.\frac{\partial k_2}{\partial y}\right|}_{x=0,y=0}=-\frac{4y_0(a_x{x}_0+a_z{z}_0)}{{(x_0^2+y_0^2+z_0^2)}^{1.5}}.\end{array}$$

  (16)
  After the range perturbation processing, the range curvatures for different targets in a given range bin are uniform and can be compensated simultaneously.
• Multiplying factor $H_2(f_t,f_{\tau })$ with signal obtained by step 2 to compensate for the range curvature correction (RCC) and second compression correction (SRC) in 2-D frequency domain, where

  $$\begin{array}{c}\hfill H_2(f_t,f_{\tau })=exp\left(-j\pi ({\mathsf{\Phi }}_{rcc}\left(f_t\right)f_{\tau }+{\mathsf{\Phi }}_{src}\left(f_t\right)f_{\tau }^2)\right).\end{array}$$

  (17)
  ${\mathsf{\Phi }}_{rcc}\left(f_t\right)$ and ${\mathsf{\Phi }}_{src}\left(f_t\right)$ stand for the correction phases for RCC and SRC, respectively, where

  $$\begin{array}{cc}\hfill {\mathsf{\Phi }}_{rcc}\left(f_t\right)=& \frac{{2r(t;x,y)|}_{t=0,x=0,y=0}}{c}-\frac{k_{10}^2}{2c{k}_2}-\frac{k_{10}^3{k}_{30}}{4c{k}_{20}^3}\hfill \\ \hfill & +\frac{c{f}_t^2}{8k_{20}{f}_c^2}+\frac{3\lambda k_{10}{k}_{30}{f}_t^2}{16k_{20}^3{f}_c}+\frac{{\lambda }^2{k}_{30}{f}_t^3}{16k_{20}^3{f}_c},\hfill \end{array}$$

  (18)

  and

  $$\begin{array}{c}\hfill {\mathsf{\Phi }}_{src}\left(f_t\right)=\frac{3{\lambda }^2{k}_{30}{f}_t^3}{16k_{20}^3{f}_c^2}-\frac{c{f}_t^2}{4k_{20}{f}_c^3}+\frac{3\lambda k_{10}{k}_{30}{f}_t^2}{8k_{20}^3{f}_c^2}.\end{array}$$

  (19)
  In (18) and (19),

  $$\begin{array}{c}\hfill {r(t;x,y)|}_{t=0,x=0,y=0}=\sqrt{x_0^2+y_0^2+z_0^2}.\end{array}$$

  (20)
  In addition, $k_{20}$ and $k_{30}$ are given by (21) and (22), respectively.

  $$\begin{array}{cc}\hfill k_{20}=k_2{|}_{x=0,y=0}& ={\left.\frac{{\partial }^{2}r(t;x,y)}{2\partial t^2}\right|}_{t=y/v,x=0,y=0}\hfill \\ \hfill & =\frac{v^2+2a_x{x}_0+2a_z{z}_0}{2\sqrt{x_0^2+y_0^2+z_0^2}}-\frac{v^2{y}_0^2}{4{(x_0^2+y_0^2+z_0^2)}^{1.5}}\hfill \end{array}$$

  (21)

  $$\begin{array}{cc}\hfill k_{30}=k_3{|}_{x=0,y=0}& ={\left.\frac{{\partial }^{3}r(t;x,y)}{6\partial t^3}\right|}_{t=y/v,x=0,y=0}\hfill \\ \hfill & =\frac{v^3{y}_0^3}{2{(x_0^2+y_0^2+z_0^2)}^{2.5}}-\frac{v{y}_0(v^2+2a_x{x}_0+2a_z{z}_0)}{2{(x_0^2+y_0^2+z_0^2)}^{1.5}}\hfill \end{array}$$

  (22)
• The high-order term of RCM is compensated by

  $$\begin{array}{c}\hfill H_{rcm2}(t,f_{\tau })=exp\left(j\frac{4\pi (f_c+f_{\tau })}{c}{r}_h\left(t\right)\right),\end{array}$$

  (23)

  where

  $$\begin{array}{c}\hfill r_h{\left(t\right)=r\left(t\right)|}_{x=0,y=0}-k_{10}t-k_{20}{t}^2\end{array}$$

  (24)

  is the high-order term of RCM.

3.3. Range Resolution Adjustment

In subsection B, most but not all of the RCM has been corrected. To deal with the residual RCM, we manage to adjust the range resolution to a coarse one so that the residual RCM is within one range resolution. With iterative RCM correction scheme proceeding, the estimation of $a_x$ and $a_z$ become more and more accurate, and the coarse resolution can be adjusted to a finer one gradually [38,39]. The coarse-to-fine range resolution can be implemented as

$$\begin{array}{c}\hfill {\rho }_r^{k+1}={\rho }_r\left(1+a{b}^k\right),\end{array}$$

(25)

where $a>1$ is a constant, which stands for the range resolution amplification factor; $b<1$ is the changing rate of range resolution; the range resolution of the kth iteration is denoted by ${\rho }_r^k$; and ${\rho }_r$ is the range resolution of the original SAR data. Here, we choose $a=10$ and $b=0.4$ in this paper. Due to the fact that the range resolution is linear to the inverse of the bandwidth $B_r$, as

$$\begin{array}{c}\hfill {\rho }_r=\frac{c}{2B_r},\end{array}$$

(26)

the range resolution adjustment can be implemented through changing the bandwidth of the range signal. Specifically speaking, the bandwidth $B_r$ can be adjust by filtering the SAR data in the range frequency domain. At the $k+1$th iteration, the frequency filtering width is set to be

$$\begin{array}{c}\hfill B_r^{k+1}=\frac{B_r}{1+a{b}^k}.\end{array}$$

(27)

3.4. Motion Modelling Based on Nlcs

As analysed in [40], the targets in different range bin—with the same ${r(t;x,y)|}_{t=0}$—will be relocated into the same range gate after LRCMC, which introduce azimuth-variance with respect to Doppler parameters. For the purpose of evaluating the azimuth variance of the Doppler parameters, the relationship between coordinates x and y in a certain range bin is given by

$$\begin{array}{c}\hfill r_0=\sqrt{{(x-x_0)}^2+{(y-y_0)}^2+z_0^2},\end{array}$$

(28)

where $r_0$ is the slant range for the given range bin. The corresponding meaningful solution for (28) is

$$\begin{array}{c}\hfill x=f(r_0,y)=x_0-\sqrt{r_0^2-{(y-y_0)}^2-z_0^2},\end{array}$$

(29)

To modulate the azimuth-variant Doppler parameters and focus the curved trajectory squinted SAR echo, we model the Doppler centroid $f_{dc}$ and Doppler FM rate $f_{dr}$ as

$$\begin{array}{c}\hfill f_{dc}=f_{dc0}+f_{dc1}y,\end{array}$$

(30)

and

$$\begin{array}{c}\hfill f_{dr}=f_{dr0}+f_{dr1}y+f_{dr2}{y}^2.\end{array}$$

(31)

According to the range history defined in (3) and its polynomial model, the parameters in (30) and (31) can be obtained by

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{f}_{dc0}=(-2k_1/\lambda ){|}_{x=f(r_0,y),y=0}\\ f_{dr0}=(-2k_2/\lambda ){|}_{x=f(r_0,y),y=0}\\ f_{dc1}=(-2\partial k_1(x,y)/\left(\lambda \partial y\right)){|}_{x=f(r_0,y),y=0}\\ f_{dr1}=(-2\partial k_2(x,y)/\left(\lambda \partial y\right)){|}_{x=f(r_0,y),y=0}\\ f_{dr2}=(-{\partial }^2{k}_2(x,y)/\left(\lambda \partial y^2\right)){|}_{x=f(r_0,y),y=0}\end{array}\right.,\end{array}$$

(32)

where

$$\begin{array}{c}\hfill k_1{\mid }_{x=f(r_0,y),y=0}=\frac{v{y}_0}{r_0},\end{array}$$

(33)

and

$$\begin{array}{cc}\hfill & k_2{|}_{x=f(r_0,y),y=0}\hfill \\ \hfill & =\frac{v^2+2a_x\sqrt{r_0^2-y_0^2-z_0^2}+2a_z{z}_0}{r_0}-\frac{v^2{y}_0^2}{r_0^3}.\hfill \end{array}$$

(34)

In addition, the partial derivatives of $k_1(x,y)$ and $k_2(x,y)$ are

$$\begin{array}{cc}\hfill & {\left.\frac{\partial k_1(x,y)}{\partial y}\right|}_{x=f(r_0,y),y=0}\hfill \\ \hfill & =\frac{2a_x\sqrt{r_0^2-y_0^2-z_0^2}+2a_z{z}_0}{v{r}_0}-\frac{v{y}_0^2}{r_0^3},\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill & {\left.\frac{\partial k_2(x,y)}{\partial y}\right|}_{x=f(r_0,y),y=0}\hfill \\ \hfill & =\frac{a_x{y}_0}{r_0\sqrt{r_0^2-y_0^2-z_0^2}}+\frac{3v^2{y}_0^3}{2r_0^5}\hfill \\ \hfill & -\frac{0.5v^2{y}_0-3y_0(a_x\sqrt{r_0^2-y_0^2-z_0^2}+a_z{z}_0)}{r_0^3},\hfill \end{array}$$

(36)

and

$$\begin{array}{cc}\hfill {\left.\frac{{\partial }^2{k}_2(x,y)}{\partial y^2}\right|}_{x=f(r_0,y),y=0}=& \frac{12(a_x^2+a_z^2)}{v^2{r}_0}-\frac{12a_z^2{z}_0^2}{v^2{r}_0^3}-\frac{15v^2{y}_0^4}{r_0^7}-\frac{12a_x^2(r_0^2-y_0^2-z_0^2)}{v^2{r}_0^3}\hfill \\ \hfill & +\frac{36a_z{y}_0^2{z}_0}{r_0^5}-\frac{v^2}{r_0^3}-\frac{24a_x{a}_z{z}_0{(r_0^2-y_0^2-z_0^2)}^{0.5}}{v^2{r}_0^3}\hfill \\ \hfill & -\frac{2a_x}{{(r_0^2-y_0^2-z_0^2)}^{0.5}{r}_0}+\frac{36a_x{y}_0^2{(r_0^2-y_0^2-z_0^2)}^{0.5}}{r_0^5}\hfill \\ \hfill & -\frac{12a_x{y}_0^2}{{(r_0^2-y_0^2-z_0^2)}^{0.5}{r}_0^3}-\frac{2a_x{y}_0^2}{{(r_0^2-y_0^2-z_0^2)}^{1.5}{r}_0}.\hfill \end{array}$$

(37)

By utilising the Doppler parameter model established above, an azimuth nonlinear chirp scaling (NLCS) [41] is implemented to equalise the translational variant Doppler parameters and get the focused imaging result as follows.

• Perform a fourth-order filtering processing with respect to range-Doppler domain by multiplying the phase

  $$\begin{array}{c}\hfill \varphi \left(f_t\right)=exp\left\{j\pi (Y_3{f}_t^3+Y_4{f}_t^4)\right\}.\end{array}$$

  (38)
• Equalise the Doppler parameters by multiplying equalisation factor $\phi \left(t\right)$ with respect to 2 dimensional time domain

  $$\begin{array}{c}\hfill \phi \left(t\right)=exp\{j\pi q_2{t}^2+j\pi q_3{t}^3+j\pi q_4{t}^4\}.\end{array}$$

  (39)
  The expressions of the parameters in the fourth-order filtering and equalisation factors are

  $$\begin{array}{c}\hfill \left\{\begin{array}{c}{q}_2=-f_{dc1}\\ q_3=\frac{f_{dc1}N}{3{(f_{dr0}-f_{dc1})}^2}\\ Y_3=\frac{f_{dr1}}{3(f_{dc1}-f_{dr0})f_{dr0}{f}_{dc1}}\\ Y_4=-\frac{L/6+M(f_{dr0}-f_{dc1})/4}{{(f_{dr0}-f_{dc1})}^3{f}_{dr0}^2{f}_{dc1}}\\ q_4=\frac{M/4+{(f_{dr0}-f_{dc1})}^3{f}_{dr0}^3{f}_{dc1}{Y}_4}{f_{dc1}-f_{dr0}}\end{array}\right..\end{array}$$

  (40)

  where

  $$\begin{array}{c}\hfill \left\{\begin{array}{c}N=f_{dr1}(f_{dr0}-f_{dc1})\\ M=-3Y_3{f}_{dc1}{f}_{dr1}{f}_{dr0}^2-3q_3{f}_{dr1}\\ L=6q_3{f}_{dr1}(f_{dr0}-f_{dc1})\\ +6Y_3{f}_{dc1}{f}_{dr1}{f}_{dr2}(f_{dr0}{f}_{dc1}-f_{dc1}^2)\\ -[f_{dr2}{(f_{dr0}-f_{dc1})}^2-f_{dr1}^2(f_{dr0}-f_{dc1})]\end{array}\right..\end{array}$$

  (41)
• Compress the azimuth signal by multiplying azimuth compression factor $exp\left({\mathsf{\Phi }}_{az}\left(f_t\right)\right)$ in azimuth frequency domain, where ${\mathsf{\Phi }}_{az}\left(f_t\right)$ is

  $$\begin{array}{cc}\hfill {\mathsf{\Phi }}_{az}\left(f_t\right)=& \pi \frac{1}{q_2+f_{dr0}}{f}_t^2\hfill \\ \hfill & -\pi \frac{Y_3{f}_{dr0}^3+q_3}{{(q_2+f_{dr0})}^3}{f}_t^3\hfill \\ \hfill & -\pi \frac{Y_4{f}_{dr0}^4+q_4}{{(q_2+f_{dr0})}^4}{f}_t^4.\hfill \end{array}$$

  (42)

Afterwards, the imaging result can be obtained, which is discretely represented by $z(m,n)$, where $m=1,2,\dots ,M$ indexes the range time $\tau$ and $n=1,2,\dots ,N$ stands for the samples taken along azimuth time t, M and N are the range and azimuth sample numbers, respectively.

Then, we want to use some metric to evaluate the imaging result $z(m,n)$, and formulate the motion parameter estimation problem into an optimisation problem. Lots of image evaluation metrics, including entropy [42,43,44] and sharpness [45,46], have been applied to SAR imaging application. Here, we choose entropy due to its robustness in quality assessment for SAR images. It is generally viewed that a better focusing quality corresponds to a smaller entropy, and therefore the motion parameter estimation can be implemented by finding the global minimum entropy solution, where entropy of the final image is

$$\begin{array}{c}\hfill E=lnS-\frac{1}{S}\sum _{m=0}^{M-1}\sum _{n=0}^{N-1}{\left|z(m,n)\right|}^{2}ln{\left|z(m,n)\right|}^2\end{array}$$

(43)

In (43), S is the total energy of the imaging result $z(m,n)$,

$$\begin{array}{c}\hfill S=\sum _{m=0}^{M-1}\sum _{n=0}^{N-1}{\left|z(m,n)\right|}^2\end{array}$$

(44)

Then, the optimisation model to estimate the motion parameters $a_x$ and $a_z$ under the criterion of minimum entropy is

$$\begin{array}{c}\hfill \underset{a_x,a_z}{argmin}E(a_x,a_z).\end{array}$$

(45)

3.5. Motion Estimation Based on Differential Evolution

In order to search the global minimisation solution to the problem in (45), we use differential evolution (DE) [47,48] optimisation scheme. DE algorithm is an algorithm suitable for solving the global optimal solution in multidimensional space. The algorithm uses differential information of a certain scale among individuals to mutate, a binomial or exponential method to cross, and finally greedy methods to get good individuals into the next generation. The specific procedures to apply DE to solve problem (45) are shown in Figure 4.

In DE algorithm, the global minimisation solution is searched in a D-dimensional parameter space, where $D=3$ for the proposed method. It starts with a randomly initialised D-dimensional parameter vector, also known as genome/chromosome. Suppose the subsequent generations of DE are denoted by $g=0,1,\dots ,G$, and the ith vector in the population at gth generation is

$$\begin{array}{c}\hfill {\mathbf{A}}_{i,g}=\left[a_{x,i,g},a_{z,i,g}\right].\end{array}$$

(46)

As for the acceleration parameters to be estimated here, the values should be limited to a certain range, usually because these parameters are associated with a physical component or with a natural boundary. The initial population (at $g=0$) should cover this range as much as possible by evenly and randomly allocating individuals within the search space under specified minimum and maximum bounds constraints:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\mathbf{A}}_{max}=\left[a_{xmax},a_{zmax}\right]\\ {\mathbf{A}}_{min}=\left[a_{xmin},a_{zmin}\right]\end{array}\right..\end{array}$$

(47)

Then, the procedures to estimate $a_x$ and $a_z$ based on the model in (45) are given as follows.

• Initialisation of the population: Randomly initialise the population as

  $$\begin{array}{c}\hfill a_{j,i,0}=a_{jmin}+d_{i,j}(a_{jmax}-a_{jmin}),\end{array}$$

  (48)

  where $j=x$ or z, and $d_{i,j}$ is a random valuable between 0 and 1.
• Difference-vector based mutation operator: Three mutually distinct vectors, ${\mathbf{A}}_{r_1^i,g}$, ${\mathbf{A}}_{r_2^i,g}$ and ${\mathbf{A}}_{r_3^i,g}$, are randomly selected from the population and the donor vector for the ith individual is set as

  $$\begin{array}{c}\hfill {\mathbf{V}}_{i,g}={\mathbf{A}}_{r_1^i,g}+F({\mathbf{A}}_{r_2^i,g}-{\mathbf{A}}_{r_3^i,g})\end{array}$$

  (49)

  where F is a scaling factor within $[0.4,1]$.
• Crossover/Recombination operator: This operation uses a recombination operator to exchange donor vector ${\mathbf{V}}_{i,g}$ and the target vector ${\mathbf{A}}_{i,g}$ to generate a test vector ${\mathbf{U}}_{i,g}=\left[u_{x,i,g},u_{y,i,g},u_{z,i,g}\right]$, where

  $$\begin{array}{c}\hfill u_{j,i,g}=\left\{\begin{array}{c}{v}_{j,i,g},\mathrm{if}{d}_{i,j}\le C_r\mathrm{or}j=j_{rand}\\ a_{j,i,g},\mathrm{otherwise}\end{array}\right.\end{array}$$

  (50)

  where $C_r$ is a pre-designed crossover factor with $[0,1]$, and $j_{rand}=x$, y or z is a randomly chosen index ensuring that at least one dimension of test vector ${\mathbf{U}}_{\mathbf{i},\mathbf{g}}$ is chosen from the donor vector ${\mathbf{V}}_{i,g}$.
• NLCS processing and entropy calculation: Carry out the NLCS processing and calculate the corresponding image entropy for the imaging result based on the parameter vectors ${\mathbf{A}}_{\mathbf{i},\mathbf{g}}$ and ${\mathbf{U}}_{\mathbf{i},\mathbf{g}}$, and the entropies corresponding to ${\mathbf{A}}_{\mathbf{i},\mathbf{g}}$ and ${\mathbf{U}}_{\mathbf{i},\mathbf{g}}$ are denoted as $e_A$ and $e_U$, respectively.
• Selection operator: According to the requirements of minimising the objective function (image entropy), the next-generation population is selected from the target vector ${\mathbf{A}}_{i,g}$ and the test vector ${\mathbf{U}}_{i,g}$. The specific selection process is as follows.

  $$\begin{array}{c}\hfill {\mathbf{A}}_{i,g}=\left\{\begin{array}{c}{\mathbf{A}}_{i,g},\mathrm{if}{e}_{\mathbf{A}}<e_{\mathbf{U}}\\ {\mathbf{U}}_{i,g},\mathrm{otherwise}\end{array}\right..\end{array}$$

  (51)
  Go back to the second step until meeting the stop condition.
• Stop condition: Continue the steps from (2) to (5) until the generation number g reaches its maximum value G.

4. Discussion

In this section, we make several comments about our work:

• Motion model: In this paper, for the simplicity and understandability of the formulas, we only consider accelerations of the radar platform in the motion model. However, in principle, the proposed imaging strategy can be formulated based on a more complex motion model with more higher order terms of the motion trajectory. A more complicated motion model will result in a much more complex derivation of Doppler parameters $f_{dc}$ and $f_{dr}$ in (30) and (31) and a higher-dimensional estimation during DE algorithm.
• Imaging processing: To cope with the translational variant RCM and azimuth modulation in curved trajectory squinted SAR imaging, there have been some proposed imaging algorithms. In this paper, we develop range perturbation [23] and NLCS [41] processes here, which have some approximations and assumptions. We would like to make our effort to combine other more accurate imaging processes—such as Omega-k based algorithm [12,26], NLCS jointly combined with 2-D singular value decomposition (SVD) [25] and time-domain algorithm [17,21]—with our proposed motion modelling and estimation strategy in future work.
• Computational cost: Due to the need of iteratively motion estimation using DE and coarse-to-fine RCM correction scheme, the computational cost of the proposed method is much larger than the existing imaging algorithm for squinted SAR with curved trajectory. However, the proposed method deals with a much more complicated but practical issue—imaging without the accurate motion information of the radar platform. It is a trade-off between complexity of the problem and computational cost. However, with the improvement of computing power, we believe that in the future, computational cost will not be a hindrance to the widespread application of the proposed algorithm.

5. Results

In this section, numerical simulations are conducted to evaluate the imaging performance and verify the effectiveness of the proposed method. The imaging performance of the proposed method are further compared with a state-of-the-art method for curved trajectory SAR imaging, i.e., NLCS-based 2-D singular value decomposition (SVD) [25]. Note that the NLCS-based 2-D SVD algorithm is for the case that trajectory is known, and there is no published literature for the case without any information of the movement. Consequently, for the simulation data, the NLCS-based 2-D SVD is implemented with 2 different cases: (1) with accurate movement information and (2) with inaccurate movement information.

It should be stressed that, for the simulated data, our proposed method is implemented without any information of the movement, and we cannot expect it outperforming the NLCS-based 2-D SVD algorithm with accurate movement information. Actually, the similarity of the imaging results between the proposed algorithm and NLCS-based 2-D SVD algorithm with accurate movement information can prove the effectiveness of the proposed method as the proposed method does not use any prior information of the movement. While for the experimental data, the better performance of the proposed method than the NLCS-based 2-D SVD algorithm with INS data can verify the effectiveness of the proposed method.

5.1. Experiment on Point-Like Targets

In the simulation, the imaging algorithms are tested on simulated data of point targets. The spatial distribution of the point targets are illustrated in Figure 5. Thirty-six ideal point targets are evenly distributed in the x-y plane with 500 m distance between adjacent targets both in x and y directions. Four point targets, $P_1$, $P_2$, $P_3$ and $P_4$, are selected for the presentation in the following.

The corresponding system and geometry parameters for the simulation are listed in

Table 1. Simulation parameters.

$\mathbf{Parameter}$	$\mathbf{Value}$
$f_c$: Carrier frequency	$9.6\mathrm{GHz}$
$B_r$: Signal bandwidth	$100\mathrm{MHz}$
$T_r$: Signal time width	$1.5¯\mathrm{s}$
$f_s$: Range time sampling rate	$150\mathrm{MHz}$
$T_s$: Synthetic aperture time	$0.7\mathrm{s}$
v: Radar platform velocity	86 m/s
$PRF$: Pulse repetition frequency	$200\mathrm{Hz}$
$x_0$: Radar position in x direction at initial time	$-1\mathrm{km}$
$y_0$: Radar position in y direction at initial time	$-3\mathrm{km}$
$z_0$: Radar position in z direction at initial time	$2\mathrm{km}$
$a_x$: acceleration in x direction	$2.5\mathrm{m}/{\mathrm{s}}^2$
$a_z$: acceleration in z direction	$1.9\mathrm{m}/{\mathrm{s}}^2$

.

Note that during the implementation of the NLCS-based 2-D SVD algorithm with inaccurate movement information, we set the inaccurate accelerations—which are used in the imaging processing—are $2.45\mathrm{m}/{\mathrm{s}}^2$, $1.37\mathrm{m}/{\mathrm{s}}^2$ and $1.88\mathrm{m}/{\mathrm{s}}^2$ with respect to x, y and z directions, respectively. According to the stop condition in (25), the proposed algorithm stops with 12 iterations.

First, in order to show the effectiveness of the proposed algorithm, we show the RCM corrected data of the NLCS-based 2-D SVD algorithm with inaccurate movement information, NLCS-based 2-D SVD algorithm with accurate movement information and the proposed algorithm as shown in Figure 6a–c, respectively.

It can be observed that the RCM trajectories are beyond one range bin for the correction result of NLCS-based 2-D SVD algorithm with inaccurate movement information. Comparatively, RCM can be accurately corrected by the proposed method as the curved RCM trajectories have been changed into straight lines as shown in Figure 6c. The NLCS-based 2-D SVD algorithm with accurate movement information can correct the RCM as well.

The imaging results are obtained by exploiting the NLCS-based 2-D SVD algorithm with inaccurate movement information, NLCS-based 2-D SVD algorithm with accurate movement information and the proposed algorithm as shown in Figure 7a–c, respectively. In Figure 7a, the imaging result of NLCS-based 2-D SVD algorithm using inaccurate movement information shows that all the point targets are severely defocused in the presented case. Even a small error occurring the motion parameters will result in a serious degradation on the imaging quality. Compared with Figure 7a, the image qualities in Figure 7b is clearly improved as the accurate movement information is utlized. However, as for the proposed algorithm, we don’t need any prior information of the motion of the radar platform, and we can still get the well-focused imaging result as shown in Figure 7c.

To further validate the proposed imaging approach, contour plots of the impulse response function on four different point targets, $P_1$, $P_2$, $P_3$ and $P_4$, are shown in the following. To give the details, the results are interpolated eight times. The contours of point targets $P_1$, $P_2$, $P_3$ and $P_4$ processed by the NLCS-based 2-D SVD algorithm with inaccurate movement information, NLCS-based 2-D SVD algorithm with accurate movement information and the proposed algorithm are shown in Figure 8, Figure 9 and Figure 10, respectively. Apparently, the NLCS-based 2-D SVD algorithm cannot focus the targets without accurate movement information, but we can have a well-focused results by the proposed method without any prior information of the radar movement.

The quantitative analysis of the measured imaging metrics are given in

Table 2. Imaging quality evaluation.

		Azimuth			Range
		PSLR(dB)	ISLR(dB)	IRW(m)	PSLR(dB)	ISLR(dB)	IRW(m)
Proposed method	P1	−13.20	−10.19	1.35	−13.19	−10.22	1.35
	P2	−13.23	−10.23	1.31	−13.22	−10.20	1.34
	P3	−13.22	−10.22	1.32	−13.24	−10.23	1.33
	P4	−13.14	−10.15	1.37	−13.17	−10.17	1.36
NlCS based 2-D SVD algorithm (with accurate movement information)	P1	−13.24	−10.22	1.33	−13.23	−10.22	1.34
	P2	−13.23	−10.23	1.31	−13.24	−10.23	1.33
	P3	−13.24	−10.23	1.30	−13.24	−10.24	1.33
	P4	−13.22	−10.21	1.35	−13.22	−10.22	1.35

to further evaluate the performance of the proposed method. It should be noted that we do not evaluate the imaging results of the NLCS-based 2-D SVD algorithm with inaccurate movement information as they are obviously unfocused in as shown in Figure 7a and Figure 8. The theoretical values of peak sidelobe ratio (PSLR) and integrated sidelobe ratio (ISLR) are −13.24 dB and −10.23 dB, respectively. The 3 dB impulse response width (IRW) is 1.3 m in the azimuth direction and 1.33 m in the range direction. It can be observed from Table 2 that both of the imaging performance metrics of the proposed algorithm and NLCS-based 2-D SVD algorithm with accurate movement information are close to the theoretical values. The azimuth resolution has a maximum broadening of about 5.4% for the proposed algorithm, and comparatively, it is 3.8% for the NLCS-based 2-D SVD algorithm with accurate movement information. This indicates that the proposed method without movement information can achieve a similar imaging quality as the NLCS-based 2-D SVD algorithm with accurate movement information, which verify the effectiveness of the proposed method.

5.2. Experiment on Extended Targets

The results of the proposed technique with state-of-the-art curved trajectory squinted SAR imaging algorithm using SAR data with extended targets are presented in this section. The system and geometry parameters of the experiment is shown in

Table 3. Experimental parameters.

$\mathbf{Parameter}$	$\mathbf{Value}$
$f_c$: Carrier frequency	$9.6\mathrm{GHz}$
$B_r$: Signal bandwidth	$300\mathrm{MHz}$
$T_r$: Signal time width	$1.5¯\mathrm{s}$
$f_s$: Range time sampling rate	$400\mathrm{MHz}$
$T_s$: Synthetic aperture time	$0.7\mathrm{s}$
v: Radar platform velocity	147 m/s
$PRF$: Pulse repetition frequency	$1000\mathrm{Hz}$
$x_0$: Radar position in x direction at initial time	$-0.8\mathrm{km}$
$y_0$: Radar position in y direction at initial time	$-1.6\mathrm{km}$
$z_0$: Radar position in z direction at initial time	$1.3\mathrm{km}$

. Real SAR image is utilised as the extended targets to generate SAR echo.

According to the stop condition in (25), the proposed algorithm stops with 16 iterations. The imaging results of different imaging algorithms are shown in Figure 11, where the image size is 1.8 km × 0.75 km (Range × Azimuth). The horizontal direction is range, and the vertical direction is azimuth. The imaging area is an airport and some neighbouring buildings. The imaging results of NLCS-based 2-D SVD algorithm with inaccurate movement information, NLCS-based 2-D SVD algorithm with accurate movement information and the proposed algorithm are shown in Figure 11a–c, respectively. Note that the in the NLCS-based 2-D SVD algorithm, the motion information of the radar platform is measured from the INS, which is not accurate enough for SAR imaging. In Figure 11a, the imaging result suffers from sever blurring due to the error of the INS data. Compared with Figure 11a, an obvious focusing improvement of the imaging result processed by the proposed algorithm can be observed. The outline of the buildings and roads can be easily recognised. Some strong scatters are perfectly focused with many details being easily identified. The proposed algorithm has a very similar performance and focusing quality as the NLCS-based 2-D SVD algorithm with accurate movement information.

6. Conclusions

In terms of the squinted SAR with curved trajectory, this paper makes an attempt of imaging the SAR echo without any movement information of the radar platform. A novel motion modelling and estimation algorithm is proposed to implement SAR imaging and motion information estimation at the same time. In the proposed method, range perturbation and NLCS processes are utilised to solve the problems stemming from the translational variant RCM and Doppler parameters, while DE algorithm is applied to estimate the motion parameters. We empirically compare the proposed algorithms with state-of-the-art imaging algorithms for SAR with curved trajectory on synthetic SAR data, and the results verify the effectiveness of the proposed method.

Some future challenges on our work lie in the following aspects.

• In order to improve the image quality further, we plan to conduct some more accurate imaging processes with our proposed motion modelling and estimation strategy in future work.
• In order to reduce the computational cost of the proposed algorithm, the authors would like to develop some other efficient global optimisation algorithm to replace DE during the estimation of motion parameters.