**1. Introduction**

Higher spatial resolution is an important development direction of synthetic aperture radar (SAR). Recent SAR systems are capable of resolutions in the decimeter regime. This requires the usage of large range bandwidth and wide azimuth beamwidth. The high-resolution, together with the all-weather day-and-night imaging capabilities, is turning SAR into an ideal tool for regular mapping and monitoring applications [1,2]. Moreover, microwaves can penetrate into vegetation and even the ground up to a certain depth [3]. The penetration capabilities depend on the carrier frequencies as well as on the complex dielectric constants, densities and conductivities of the observed targets. The high frequencies, like the X-band (8–12 GHz), show typically a high attenuation and are mainly backscattered on the top of the vegetation. Low frequencies, like P- and L-band (0.23∼1 GHz and 1∼2 GHz, respectively) [4], usually penetrate deep into vegetation, snow and ice. A high-resolution low frequency SAR system refers to a SAR system which operates with a low frequency (P- or

L-band) signal with a large fractional bandwidth (>0.2, i.e., the ultra-wideband SAR [2,5–8]) and a wide antenna beamwidth (corresponding to high azimuth resolution in the decimeter regime). The fractional bandwidth is defined by the ratio of the signal bandwidth to the center frequency. The combination of low frequency with large bandwidth and wide beam allows SAR to obtain high-resolution images of concealed targets, with the capability of penetrating the ground or foliage surface, thus it has broad applications for both military and civil purposes in recent years [9–12]. However, the large fractional bandwidth and long azimuth integration time used in high-resolution low frequency SAR bring new challenges to get high-quality images by the conventional image formation.

A crucial problem is the serious coupling between the range and azimuth frequencies in the phase of high-resolution low frequency SAR transfer function [13]. In two-dimension (2D) frequency domain, the phase is range-dependent and can be decomposed into two parts: the range-independent terms and the range-dependent terms. Different algorithms make different approximations of these two parts. At low frequencies, many of the simplifying assumptions made by traditional algorithms are not valid, such as range-Doppler algorithm [14] and chirp scaling algorithm (CSA) [15], resulting in serious image degradations as blurring and resolution loss. The problem stems from high-order range-azimuth phase coupling. Several approaches have been used in the processing of this type of SAR system. The time-domain algorithms and the wavenumber-domain Omega-K (*ω* − *k*) algorithm [16–18] are two common options. The time-domain algorithms, often referred to as backprojection (BP) class algorithms [19], are most accurate and can easily adapt to all SAR configurations. Due to their computational complexity and the poor ability to integrate accurate autofocus algorithm into its imaging process, their use is restricted. The *ω* − *k* algorithm is an ideal solution without approximation in range cell migration (RCM) correction, which can focus data up to very high-resolution values regardless of their azimuth and range bandwidth. However, it is only applicable to spaceborne SAR data with a straight sensor trajectory and can only perform the range-independent motion compensation (MoCo). In addition, the Stolt interpolation makes it to be time-consuming . The extended Omega-K algorithm (EOKA) [20,21] is proposed to integrate the high precise range-dependent MoCo but it is still inefficient due to the Stolt interpolation.

For efficiency reasons, the chirp scaling class algorithms are still attractive. Efforts have been made to modify the CSA to process the low frequency SAR data. Without the interpolation, the chirp scaling class algorithms are effective and phase-preserving. The nonlinear chirp scaling algorithm (NCSA) [13] is proposed to take into account the cubic range-independent coupling phase and the range dependence of secondary range compression (SRC) term, which has better performance than CSA on processing the raw data of highly squint or low frequency SAR cases. Whereas the range dependence of cubic- and higher order terms are neglected. Some modified NCSAs are proposed in References [22–24], which resolve the high-order range-independent coupling terms. However, the cubic and higher order range-dependent coupling terms are still neglected. Besides, the first order approximation of range frequency modulation (FM) rate will introduce a quadratic phase error in the spectrum and degrades the focusing quality of image. Thus, the range focusing depth is restricted, which shows that the NCSA may not be suitable for the high-resolution low frequency SAR processing. A helpful comparison of the BP class algorithm, *ω* − *k* algorithm, EOKA and NCSA can be found in References [22,25]. In References [26,27], a generalized chirp scaling algorithm (GCSA) is developed for the SAR systems operating on wide bandwidths at low frequencies. The GCSA is an efficient arbitrary-order CSA that processes the data using the appropriate number of the approximation terms. Both the higher range-independent terms and range-dependent terms are considered. The GCSA efficiently extends the utility of frequency domain processing for high-resolution low frequency SAR systems. However, the imaging quality of GCSA also decreases as the fractional bandwidth get larger and beamwidth gets wider. In addition, the improvement of focus quality is not significant when the order is greater than 3rd and the edge targets of the range swath still have obvious degradation. Two main reasons lead to this phenomenon. One is that the residual coupling terms are still significant and the other is the error caused by the linear approximation of stationary phase point when solving

the fast Fourier transform (FFT) expression. The linear approximation makes the range-dependent high-order phase terms not effectively compensated, even if a higher-order model is used.

The aim of this study is to overcome these two limits of GCSA and propose a more accurate approach than the GCSA for processing wide-swath, high-resolution low frequency SAR data. To our knowledge, the Lagrange inversion [28–31] gives the power series representation of the inverse of an analytic function, which is quite suitable for calculating the expression of stationary phase point. This paper utilizes Lagrange inversion to calculate a more precise expression of stationary phase point, while compensating for all range-independent coupling terms above 3rd order. In our approach, the high-order chirp scaling technique is extended to achieve the effect of the range-variant filtering required in high order phase terms. The experimental results show that the improved algorithm has a better focusing performance than the original GCSA. The resolution, sidelobe level and range focusing depth are significantly improved.

This paper is organized as follows. In Section 2, the signal model of high-resolution low frequency SAR is analyzed and the limitations of existing GCSA are briefly described. Then, a principle to determine the required order of range-dependent coupling phase is presented. In Section 3, an improved GCSA based on the Lagrange inversion theorem is introduced to focus the high-resolution low frequency SAR data. Focused results obtained by the conventional GCSA and the improved GCSA are presented and analyzed to verify our analysis in Section 4. A discussion is givne in Section 5. Finally, conclusions are drawn in Section 6.

#### **2. Background and Problem Statement**

#### *2.1. Signal Model*

The 2D spectrum is a key for the frequency domain algorithm development. In our analysis, we only consider the phase terms of the SAR signal and ignore the initial phase. The range-dependent SAR transfer function in the 2D frequency domain can be expressed as [14]

$$\Phi\left(f\_{\tau}, f\_{\eta}; R\_0\right) = -\frac{4\pi R\_0 f\_0}{c} \sqrt{D^2(f\_{\eta}) + \frac{2f\_{\tau}}{f\_0} + \frac{f\_{\tau}^2}{f\_0^2}} - \frac{\pi f\_{\tau}^2}{K\_r} \tag{1}$$

where *f<sup>τ</sup>* and *f<sup>η</sup>* represent the range frequency and the azimuth frequency, respectively. *R*<sup>0</sup> is the closest slant range from a point target to the radar trajectory, *f*<sup>0</sup> is the carrier frequency, *c* is the speed of light, *Kr* is the range chirp rate, *<sup>D</sup>*(*fη*) = <sup>1</sup> <sup>−</sup> *<sup>c</sup>*<sup>2</sup> *<sup>f</sup>* <sup>2</sup> *η* 4*V*<sup>2</sup> *<sup>r</sup> f* <sup>2</sup> 0 is the migration factor, *Vr* is the moving velocity of the radar platform.

The first term in Equation (1) represents the coupling relationship between the range frequency and the azimuth frequency, which is called the range-azimuth coupling term. It varies with the slant range. The second term is the range modulation term. The square root term can be expanded into a Taylor series with respect to *fτ* and kept up to the *n*th order,

$$p\left(f\_{\tau}, f\_{\eta}; n\right) = D(f\_{\eta}) + \frac{f\_{\tau}}{f\_0 D(f\_{\eta})} + \frac{D^2(f\_{\eta}) - 1}{2f\_0^2 D^3(f\_{\eta})} f\_{\tau}^2 + \sum\_{i=3}^n \gamma\_i f\_{\tau}^i \tag{2}$$

where *γ<sup>i</sup>* denotes the coefficient of the *i*th term and is given in Appendix A. The first term in Equation (2) corresponds to the azimuth modulation. The second term corresponds to the RCM (first-order coupling term). The third term denotes the SRC (second-order coupling term). The remainder higher order terms donates the high-order range-azimuth coupling. Different frequency algorithms are based on specific order approximations of this equation. For example, the CSA uses a 2nd order model and the NCSA uses a 3rd moder model. For the low frequency SAR with a small bandwidth and narrow beamwidth, quadratic approximation is enough, whereas for large bandwidth the high-order terms become significant.

*Remote Sens.* **2019**, *11*, 1874

Reconsider the phase in Equation (1). Actually, the coupling phase term can be decomposed into two parts: the range-independent terms and the range-dependent terms, namely,

$$\Phi\left(f\_{\tau}, f\_{\eta}; R\_0, n\right) = -\frac{4\pi R\_c f\_0}{c} p\left(f\_{\tau}, f\_{\eta}; n\right) - \frac{4\pi \Lambda R f\_0}{c} p\left(f\_{\tau}, f\_{\eta}; n\right) - \frac{\pi f\_{\tau}^2}{K\_r} \tag{3}$$

where *Rc* represents the slant range of scene center, Δ*R* = *R*<sup>0</sup> − *Rc* represents the difference in slant range between the point target and the scene center point. The first term donates the range-independent coupling phase and the second term donates the range-dependent coupling phase, which varies with the slant range *R*0.

#### *2.2. The Limitations of the Conventional GCSA*

The phase error due to the Taylor approximation can be expressed as

$$\Phi\_{\rm error} \left( f\_{\sf \tau\_{\sf \tau}} f\_{\sf \eta}; R\_{0 \prime} n \right) = -\frac{4 \pi R\_0 f\_0}{c} \left( \sqrt{D^2(f\_{\sf \eta}) + \frac{2f\_{\sf \tau}}{f\_0} + \frac{f\_{\sf \tau}^2}{f\_0^2}} - p(f\_{\sf \tau}, f\_{\sf \eta}; n) \right) \tag{4}$$

Note that this error gets larger when the slant range *R*<sup>0</sup> increases, the maximum range frequency *fτ*,max increases, the maximum azimuth frequency *fη*,max increases, while the center frequency *f*<sup>0</sup> decreases. In high-resolution, low frequency and far range situations, the approximation error is large and a high order model is required.

To illustrate the phase error of Taylor expansion, a numerical analysis is carried out. Assume that the center frequency *f*<sup>0</sup> equals to 600 MHz, bandwidth *Br* is 300 MHz, beamwidth is 29◦ and the target slant range *R*<sup>0</sup> is 12 km. Figure 1 shows the relationship between the phase errors and the range frequency for different order approximations. The maximal phase error of 4th order model is 1025◦ and the maximal phase error of 6th order approximation is about 81.48◦, which indicates that the coupling phase terms above 6th order still have an important influence on the image quality. High-resolution imaging methods should take these terms into account.

**Figure 1.** Phase errors of 2nd to 6th order Taylor series approximation. The point target at a range of 12 km with a center frequency of 600 MHz, a bandwidth of 300 MHz and a beamwidth of 29◦.

In Reference [26], Zaugg et al. gives a guideline to determine required order *n*th from Equation (4) for proper focusing: If less than 30% of the support band has a phase error greater than *π*/10, then one

can predict less than a 20% loss in the azimuth resolution defocusing. The conventional GCSA first compensates the 3rd to *n*th order range-independent coupling terms and then uses the high-order chirp scaling filter to compensate the range-dependent coupling phase. The coupling phase terms above *n*th order are neglected.

In the case of low-frequency high-resolution SAR with large bandwidth and wide beam, the coupling phase above 6th order still has a large proportion. Therefore, a higher order model is needed. However, when we use a model higher than 3rd order, the linear approximation of the stationary phase point when using the principle of stationary phase (POSP) will bring serious errors, which severely degrades the actual performance of high-order models in conventional GCSA. This error makes the high-order range-dependent coupling phase not to be accurately compensated, even though a higher order model is used. This is because the nonlinear FM component of signal cannot be neglected. Besides, the complexity of algorithm design will increase when using a higher order model.

Using the same parameters presented in the previous analysis, Figure 2 shows the residual range migrations of a point target at a range of 12 km after processing by the conventional GCSA. It is assumed here that the reference slant range is 10 km. The range migration crosses several range gates even using a 7th order model. This is because the conventional algorithm does not accurately compensate the range-dependent coupling terms.

**Figure 2.** Range migrations of point target at a range of 12 km after RCMC. (**a**) Algorithm in Reference [26] (6th order model). (**b**) Algorithm in Reference [26] (7th order model).

#### *2.3. New Principle to Determine the Required Order of Range Frequency*

According to the previous analysis, the phase error of 6th order model is still significant. In fact, the range-independent couplings can be compensated by the reference range phase. Figure 3 shows the range-dependent coupling phase errors of different order models. If the range-independent coupling terms are firstly compensated, the maximum range-dependent phase error of 6th order model is 13.58◦, which has a small effect on the imaging results.

Therefore, if the range-independent coupling terms are firstly compensated, we should only consider the range-dependent coupling terms to determine the required order. Thus, a new principle of the proposed method to determine the required order is given here. Assume that the *M*th order range-dependent coupling term should be taken into consideration, the phase error should meet

$$-\frac{4\pi\Lambda Rf\_0}{c}\left(\sqrt{D^2(f\_\eta) + \frac{2f\_\tau}{f\_0} + \frac{f\_\tau^2}{f\_0^2}} - p\left(f\_{\tau\tau}, f\_{\eta'}; M\right)\right)\bigg|\_{f\_\tau = \frac{R\_\tau}{2}, f\_\eta = \frac{2f\_0 V\_r}{c}\sin\theta} \le \delta\tau\tag{5}$$

where *θ* is the beamwidth, *δ* equals to 1/10 in this paper and it can be set to a larger value when the imaging accuracy is not high.

**Figure 3.** Range-dependent coupling phase errors of 2nd to 6th order Taylor series approximation. The point target at the slant range 12 km with a center frequency of 600 MHz, a bandwidth of 300 MHz, a beamwidth of 29◦ and the reference slant range of 10 km.

Inequality (5) is an important basis for the proposed algorithm. It can be seen that the phase error is positively correlated with the range bandwidth *Br*, the beamwidth *θ* and the scene size 2Δ*R* and negatively correlated with the center frequency *f*0. If the center frequency and scene size are given, the required order is determined by the range bandwidth and the azimuth beamwidth.

Figure 4 shows an example of a P-band SAR system. The center frequency is 800 MHz and the scene size equals to 600 m. The required order varies with different beamwidth and fractional bandwidth. As can be seen from Figure 4, under the given parameters, the 2nd order model can only process data with a fractional bandwidth of less than 0.34. As the beamwidth increases, the fractional bandwidth value that can be processed is smaller. For low frequency SAR systems with large bandwidths and wide beams, a higher order range-dependent coupling phase should be considered.

**Figure 4.** An example of the required order of different beamwidth and fractional bandwidth in a P-band synthetic aperture radar (SAR) system.

#### **3. The Improved GCSA Based on Lagrange Inversion Theorem**

#### *3.1. Procedure of Algorithm*

Based on the analysis in Section 2, we proposed an improved GCSA in this section. The flow chart of the proposed algorithm is shown in the Figure 5. There are two main improvements to the algorithm. The first is to perform all the high-order (≥3rd) reference phase compensation in the 2D frequency domain. The second is to use the Lagrange inversion theorem to solve the expression of the range FFT and inverse FFT (IFFT). The derivation of the proposed algorithm is given in Section 3.2. The steps of the proposed algorithm are as follows:

Step 1: Select the parameter *M* according to the principle of Equation (5). Calculate *q*<sup>2</sup> ∼ *qM*, *X*<sup>3</sup> ∼ *XM* and *C*<sup>2</sup> ∼ *CM* according to Equations (25) and (26). Then, 2D FFT is implemented to transfer the data into 2D frequency domain.

Step 2: Multiply the range-independent high-order phase correction (HOPC) and perturbation equations in the 2D frequency domain. And a range IFFT is carried out to transfer the data into the range-Doppler (RD) domain.

Step 3: Multiply the high-order chirp scaling phase function. Then, the data are transformed into 2D frequency domain by range FFT.

Step 4: Multiply the range compression function to perform the RC, SRC and bulk RCMC. Then, the data are transferred into RD domain by IFFT along the range.

Step 5: Multiply the azimuth compression and residual phase correction function. And the images can be obtained by azimuth IFFT.

**Figure 5.** The flow chart of the proposed algorithm.

#### *3.2. Theoretical Formulation*

As mentioned in Section 2, the high-order terms of high-resolution low frequency SAR account for a large proportion. If it can not be eliminated, the image will be deteriorated dramatically. In order to reduce the phase error and improve imaging accuracy, the range-independent high-order (≥3rd) *Remote Sens.* **2019**, *11*, 1874

coupling phases of Equation (3) is firstly compensated in the 2D frequency domain. The 3rd order is chosen because of the need to preserve the chirp information of the range signals. The HOPC function at the reference range is given by

$$H\_{\rm HOPC} = \exp\left[\frac{4\pi R\_c f\_0}{c} \left(\sqrt{D^2(f\_\eta) + \frac{2f\_\tau}{f\_0} + \frac{f\_\tau^2}{f\_0^2}} - D(f\_\eta) - \frac{f\_\tau}{f\_0 D(f\_\eta)} - \frac{D^2(f\_\eta) - 1}{2f\_0^2 D^3(f\_\eta)} f\_\tau^2\right)\right] \tag{6}$$

Suppose that the *M*th range-dependent coupling should be taken into account. The residual signal after HOPC can be expressed as

$$\Phi(f\_{\tau}, f\_{\eta}; R\_0) = -\frac{4\pi R\_0 f\_0}{c} D(f\_{\eta}) - \frac{4\pi R\_0}{c D(f\_{\eta})} f\_{\tau} - \frac{\pi}{K\_m} f\_{\tau}^2 - \frac{4\pi \Delta R f\_0}{c} \sum\_{i=3}^{M} \gamma\_i f\_{\tau}^i \tag{7}$$

where *Km* donates the range FM rate and can be expressed as

$$K\_{\rm ll} = \frac{K\_{\rm r}}{1 - K\_{\rm I} \frac{c \mathcal{R}\_0 f\_{\eta}^2}{2V\_r^2 f\_0^3 D \left(f\_{\eta}\right)^3}}\tag{8}$$

It can be found that the residual high-order coupling becomes zero at the reference range but the residual high-order phase at other ranges increases as the target away from the reference range. Then we filter the data with a turbulent function

$$H\_{Turb} = \exp\left(j\pi \sum\_{i=3}^{M} X\_i f\_{\pi}^i\right) \tag{9}$$

This filtering step provides an accurate accommodation of the range dependence of the high-order terms. After the turbulent compensation, the phase can be expressed as

$$\begin{split} \Phi\_{1}(f\_{\tau}, f\_{\eta}; R\_{0}) &= -\frac{4\pi R\_{0} f\_{0}}{\varepsilon} D(f\_{\eta}) - \frac{4\pi R\_{0}}{\varepsilon D(f\_{\eta})} f\_{\tau} - \frac{\pi}{K\_{\text{w}}} f\_{\tau}^{2} + \pi \sum\_{i=3}^{M} \left( X\_{i} - \frac{4\pi \Delta R\_{0}}{\varepsilon} \gamma\_{i} \right) f\_{\tau}^{i} \\ &= \phi\_{0} + \phi\_{1} f\_{\tau} + \phi\_{2} f\_{\tau}^{2} + \sum\_{i=3}^{M} \phi\_{i} f\_{\tau}^{i} \end{split} \tag{10}$$

with

$$\begin{cases} \phi\_0 = -\frac{4\pi R\_0 f\_0}{\xi} D(f\_\eta) \\ \phi\_1 = -\frac{4\pi R\_0}{\varepsilon D(f\_\eta)} = -2\pi \tau\_d \\ \phi\_2 = -\frac{\pi}{K\_m} \\ \phi\_i = \pi \left[ X\_i - 2f\_0 D(f\_\eta) \gamma\_i \Delta \tau \right] (3 \le i \le M) \end{cases} \tag{11}$$

Then the range IFFT is performed along the range direction. Based on the POSP, the relationship between *τ* and *f<sup>τ</sup>* can be expressed as

$$\tau = -\frac{1}{2\pi} (\phi\_1 + 2\phi\_2 f\_\tau + \dots + M\phi\_M f\_\tau^{M-1}) \tag{12}$$

In Reference [26], only the first order term of *fτ* is retained when solving the stationary phase point, that is, *<sup>f</sup><sup>τ</sup>* <sup>=</sup> <sup>−</sup> *<sup>π</sup> φ*2 *τ* + *<sup>φ</sup>*<sup>1</sup> 2*π* .This approximation is effective in most cases. However, as the center frequency decreases, the beamwidth and bandwidth increase, this approximation introduces severe phase errors, resulting in significant degradation throughout the image. To solve this problem, we use the Lagrange inversion to find a more accurate solution for *fτ*.

In mathematical analysis, the Lagrange inversion theorem gives the Taylor series expansion of the inverse function of an analytic function and it can be expressed as Theorem 1 [28–31].

**Theorem 1.** *Suppose w is defined as a function of z by an equation of the form w* = *h*(*z*)*, where h is analytic at a point z*<sup>0</sup> *and h* (*z*0) = 0*. Then it is possible to invert or solve the equation for z, expressing it in the form z* = *g*(*w*) *given by a power series*

$$\mathcal{g}(w) = z\_0 + \sum\_{n=1}^{\infty} \mathcal{g}\_n(w - h(z\_0))^n \tag{13}$$

*where*

$$\mathbf{g}\_{n} = \frac{1}{n!} \lim\_{z \to z\_{0}} \frac{d^{n-1}}{dz^{n-1}} \left[ \left( \frac{z - z\_{0}}{h(z) - h(z\_{0})} \right)^{n} \right] \tag{14}$$

*g*(*w*) *represents an analytic function of w in a neighbourhood of w* = *h*(*z*0)*.*

As is shown in Equation (12), *τ* is a power series of *fτ*, that is, *τ* = *h*(*fτ*). Let *z*<sup>0</sup> = 0 and *<sup>h</sup>*(*z*0) = <sup>−</sup> *<sup>φ</sup>*<sup>1</sup> <sup>2</sup>*<sup>π</sup>* . According to the Lagrange inversion theorem, the stationary phase point can be given by

$$\begin{split} f\_{7} &= -\frac{\pi}{\Phi\_{2}} \left( \tau + \frac{\Phi\_{1}}{2\pi} \right) - \frac{3\pi^{2}\Phi\_{2}}{2\Phi\_{2}^{3}} \left( \tau + \frac{\Phi\_{1}}{2\pi} \right)^{2} + \frac{\pi^{3} \left( -9\Phi\_{3}^{2} + 49\varphi\_{4} \right)}{2\Phi\_{2}^{3}} \left( \tau + \frac{\Phi\_{1}}{2\pi} \right)^{3} \\ &- \frac{5\pi^{4} \left( 27\Phi\_{3}^{3} - 24\varphi\_{2}\mathfrak{g}\_{3}\mathfrak{g}\_{4} + 4\mathfrak{g}\_{2}^{2}\mathfrak{g}\_{2} \right)}{8\Phi\_{2}^{7}} \left( \tau + \frac{\Phi\_{1}}{2\pi} \right)^{4} + \frac{3\pi^{3} \left( -189\Phi\_{3}^{4} + 252\mathfrak{g}\_{2}\mathfrak{g}\_{3}^{2}\mathfrak{g}\_{4} - 60\mathfrak{g}\_{2}^{2}\mathfrak{g}\_{3}\mathfrak{g}\_{3} + 8\mathfrak{g}\_{2}^{2} \left( -4\mathfrak{g}\_{4}^{2} + \mathfrak{g}\_{2}\mathfrak{g}\_{0} \right) \right)}{8\Phi\_{2}^{9}} \left( \tau + \frac{\Phi\_{1}}{2\pi} \right)^{5} + \dots \end{split} \tag{15}$$

The detailed derivation of Equation (15) is given in Appendix C. Therefore, the IFFT expression of Equation (10) can be expressed as

$$\Phi\_2(\tau, f\_\eta; R\_0) = \phi\_0 + A\_2(\tau - \tau\_d)^2 + A\_3(\tau - \tau\_d)^3 + \dots + A\_M(\tau - \tau\_d)^M \tag{16}$$

with (Here, only *A*<sup>2</sup> ∼ *A*<sup>6</sup> are given due to space restraints.)

$$\begin{cases} A\_{2} = -\frac{\pi^{2}}{\Phi\_{3}^{2}}\\ A\_{3} = -\frac{\pi^{3}\Phi\_{3}}{\Phi\_{2}^{2}}\\ A\_{4} = \frac{\pi^{4}\left(-9\rho\_{3}^{2} + 4\rho\_{2}\phi\_{4}\right)}{4\Phi\_{3}^{4}}\\ A\_{5} = -\frac{\pi^{5}\left(27\rho\_{3}^{2} - 24\rho\_{2}\phi\_{3}\phi\_{4} + 4\rho\_{2}^{2}\phi\_{5}\right)}{4\rho\_{2}^{7}}\\ A\_{6} = \frac{\pi^{6}\left(-189\phi\_{3}^{4} + 252\rho\_{2}\phi\_{3}^{2}\phi\_{6} - 60\rho\_{2}^{2}\phi\_{6}\phi\_{7} + 8\rho\_{2}^{2}\left(-4\rho\_{4}^{2} + \rho\_{2}\phi\_{6}\right)\right)}{8\rho\_{2}^{9}}\\ \dots \end{cases} \tag{17}$$

where *<sup>τ</sup><sup>d</sup>* = *<sup>h</sup>*(*z*0) = 2*R*0/(*cD* - *fη* ) is the time delay in RD domain. The variation of *τ<sup>d</sup>* with *f<sup>η</sup>* is called range migration, which must be removed before azimuth compression. The shape of the range migration trajectory depends on the target slant *R*0.

The high-order CS function can be given by

$$H\_{\rm CS} \left( \pi, f\_{\eta} \right) = \exp \left[ j \pi q\_2 \left( \pi - \pi\_{\rm ref} \right)^2 + j \pi \sum\_{i=3}^{M} q\_i \left( \pi - \pi\_{\rm ref} \right)^i \right] \tag{18}$$

where *<sup>τ</sup>ref* = <sup>2</sup>*Rc*/(*cD* - *fη* ) is the reference trajectory. This step aims to compensate the range-dependent coupling caused by the large fractional bandwidth and wide beamwidth. The phase after high-order CS filter can be expressed as

$$\Phi\_3(\tau, f\_\emptyset; R\_0) = \phi\_0 + A\_2(\tau - \tau\_d)^2 + \sum\_{i=3}^M A\_i(\tau - \tau\_d)^i + \pi q\_2 \left(\tau - \tau\_{\rm ref}\right)^2 + \pi \sum\_{i=3}^M q\_i \left(\tau - \tau\_{\rm ref}\right)^i \tag{19}$$

*Remote Sens.* **2019**, *11*, 1874

According to the chirp scaling principle, the desired trajectory *τ<sup>s</sup>* of the target located at range *R*<sup>0</sup> has the same shape as the reference trajectory and the relationship between *τs*, *τref* and *τ<sup>d</sup>* can be expressed as

$$\begin{cases} \tau\_{r\!\!f} = \tau\_{\!\!s} - \kappa \Delta \tau\\ \tau\_{\!\!d} = \tau\_{\!\!s} - (\kappa - 1) \Delta \tau \end{cases} \tag{20}$$

where Δ*τ* = 2Δ*R*/(*cD* - *fη* ) and *α* = *D* - *fη* /*D* - *fηc* . Using the relationship in Equation (20) and expand Equation (19) at *τs*, we obtain

$$\Phi\_{\mathsf{4}}(\mathsf{r}, f\_{\eta}; R\_{0}) = \phi\_{0} + \pi \mathsf{C}\_{0}(\Delta \mathsf{r}) + \pi \sum\_{i=1}^{M} \mathsf{C}\_{i}(\mathsf{r} - \mathsf{r}\_{\mathsf{s}})^{i} \tag{21}$$

The expressions for coefficients *Ci* are given in Appendix B. As can be seen form Equation (11), parameters *Ai* and *Ci* imply the range FM rate *Km*, which is dependent on azimuth frequency and slant range. Imaging performance is affected by approximations to *Km*. To model the range range-dependence of *Km*, we expand it at the reference slant range with Taylor series and keep up to the second order,

$$K\_{m,app} \approx K\_f + K\_s K\_f^2 \Delta \tau + K\_s^2 K\_f^3 \Delta \tau^2 \tag{22}$$

where *Kf* represents the FM rate at the reference slant range and can be expressed as

$$K\_f = \frac{K\_r}{1 - K\_r \tau\_{ref} K\_s} \tag{23}$$

with

$$K\_s = \frac{c^2 f\_\eta^2}{4v^2 f\_0^3 D \left(f\_\eta\right)^2} \tag{24}$$

According to Appendix B, the expressions of *qi*(*i* ≥ 2) and *Xi*(*i* ≥ 3) can be expressed as

⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ *q*<sup>2</sup> = *Kf* 1−*α α <sup>q</sup>*<sup>3</sup> <sup>=</sup> *<sup>K</sup>*<sup>2</sup> *<sup>f</sup> Ks*(1−*α*) 3*α <sup>X</sup>*<sup>3</sup> <sup>=</sup> *Ks*(*α*−2) 3*Kf*(*α*−1) *<sup>q</sup>*<sup>4</sup> <sup>=</sup> <sup>−</sup>*K*<sup>3</sup> *<sup>f</sup>*(9*Kf Ks*(*α*−1)*X*3+2*Ks*2−6*D*(*f<sup>η</sup>* )*f*0(*α*−1)*γ*3) 12*α <sup>X</sup>*<sup>4</sup> <sup>=</sup> <sup>9</sup>*Kf Ks*(*α*−2)*X*3−27*K*<sup>2</sup> *f*(*α*−1)*X*<sup>2</sup> <sup>3</sup>+2*Ks*2−6*D*(*f<sup>η</sup>* )*f*0(*α*−1)*γ*<sup>3</sup> 12*Kf*(*α*−1) *<sup>q</sup>*<sup>5</sup> <sup>=</sup> <sup>−</sup>45*K*<sup>6</sup> *<sup>f</sup> Ks*(*α*−1)*X*<sup>2</sup> 3+12*K*<sup>5</sup> *<sup>f</sup> <sup>X</sup>*3(*K*<sup>2</sup> *<sup>s</sup>* <sup>−</sup>3*D*(*f<sup>η</sup>* )*f*0(*α*−1)*γ*3) 20*α* −16*K*<sup>5</sup> *<sup>f</sup> Ks*(*α*−1)*X*4−4*K*<sup>4</sup> *<sup>f</sup> D*(*f<sup>η</sup>* )*f*0(3*Ksγ*3+2(*α*−1)*γ*4) 20*α <sup>X</sup>*<sup>5</sup> <sup>=</sup> <sup>−</sup>−45*K*<sup>2</sup> *<sup>f</sup> Ks*(*α*−2)*X*<sup>2</sup> <sup>3</sup>+12*Kf <sup>X</sup>*3(−*K*<sup>2</sup> *<sup>s</sup>* <sup>+</sup>10*Kf*(*α*−1)*X*4+3*D*(*f<sup>η</sup>* )*f*0(*α*−2)*γ*3) 20*Kf*(*α*−1) −4*D*(*f<sup>η</sup>* )*f*0(3*Ksγ*3+2(*α*−2)*γ*4)−16*Kf Ks*(*α*−2)*X*4+135*K*<sup>3</sup> *f*(*α*−1)*X*<sup>3</sup> 3 20*Kf*(*α*−1) ... (25)

And coefficients *Ci* becomes

$$\begin{cases} \mathsf{C}\_{1} = 0\\ \mathsf{C}\_{2} = q\_{2} + K\_{f} \\ \mathsf{C}\_{3} = q\_{3} + K\_{f}^{3}X\_{3} \\ \mathsf{C}\_{4} = q\_{4} + \frac{9}{4}K\_{f}^{5}X\_{3}^{2} + K\_{f}^{4}X\_{4} \\ \mathsf{C}\_{5} = q\_{5} + \frac{27}{4}K\_{f}^{7}X\_{3}^{3} + 6K\_{f}^{6}X\_{3}X\_{4} + K\_{f}^{5}X\_{5} \\ \ldots \end{cases} \tag{26}$$

*Remote Sens.* **2019**, *11*, 1874

After the chirp scaling operation, a range FFT is carried to transform Equation (21) into the 2D frequency domain. Similarly, *f<sup>τ</sup>* is a power series of *τ*,

$$f\_{\tau} = \frac{1}{2} \left( 2\mathbb{C}\_2 \left( \tau - \tau\_\delta \right) + 3\mathbb{C}\_3 \left( \tau - \tau\_\delta \right)^2 + \dots \text{MC}\_M \left( \tau - \tau\_\delta \right)^{M-1} \right) \tag{27}$$

Use the Lagrange inversion theorem again. Let *z*<sup>0</sup> = *τ<sup>s</sup>* and *h*(*z*0) = 0, the stationary phase point is given by

$$\begin{split} \tau &= \tau\_s + \frac{1}{C\_2} f\_\tau - \frac{3C\_3}{2C\_2^3} f\_\tau^2 + \frac{\left(\frac{9C\_3^2}{2} - 2C\_2C\_4\right)}{C\_2^3} f\_\tau^3 - \frac{\left(\frac{19C\_3^3}{8} - 15C\_2C\_3C\_4 + \frac{5}{2}C\_2^2C\_5\right)}{C\_2^3} f\_\tau^4 \\ &+ \frac{\left(\frac{50C\_3^4}{8} - \frac{19}{2}C\_2C\_3^2C\_4 + 12C\_2^2C\_4^2 + \frac{45}{2}C\_2^2C\_3C\_5 - 3C\_2^3C\_6\right)}{C\_2^3} f\_\tau^5 + \dots \end{split} \tag{28}$$

Thus, the phase after range FFT can be expressed as

$$\Phi\_5(f\_\tau, f\_\eta) = \phi\_0 + \pi \mathbb{C}\_0(\Delta \tau) - 2a\tau\_d f\_\tau + \pi \sum\_{i=1}^M E\_i f\_\tau^i \tag{29}$$

with

$$\begin{cases} E\_1 = -2\tau\_{ref}(1-\alpha) \\ E\_2 = -\frac{1}{C\_2} \\ E\_3 = \frac{C\_3}{C\_2^3} \\ E\_4 = \frac{\left(-9\mathcal{C}\_3^2 + 4\mathcal{C}\_2\mathcal{C}\_4\right)}{4C\_2^3} \\ E\_5 = \frac{\left(27\mathcal{C}\_3^3 - 24\mathcal{C}\_2\mathcal{C}\_3\mathcal{C}\_4 + 4\mathcal{C}\_2^2\mathcal{C}\_5\right)}{4C\_2^3} \\ \dots \end{cases} \tag{30}$$

The first term of Equation (29) represents the azimuth compression phase, the second term is the residual phase, the third term is the linear phase corresponding to the target position, the remainder are related to the range compression (RC), SRC and bulk RCM, which are range invariant. Thus, the bulk RCM, SRC and RC can be compensated by a conjugate multiply in the 2D frequency domain. The filtering function can be expressed as

$$H\_{\rm RC}(f\_{\tau}, f\_{\eta}) = \exp(-j\pi \sum\_{i=1}^{M} E\_i f\_{\tau}^i) \tag{31}$$

After the compensation, the phase of signal can be expressed as

$$\Phi\_{\delta}(f\_{\tau}, f\_{\eta}) = \phi\_0 + \pi \mathcal{C}\_0(\Delta \tau) - 2\pi \alpha \tau\_d f\_{\tau} \tag{32}$$

An IFFT is carried out along the range direction and the signal phase becomes

$$\Phi\_{\mathcal{T}}(\mathbf{r}, f\_{\mathcal{I}}) = \phi\_0 + \pi \mathcal{C}\_0(\Delta \mathbf{r}) \tag{33}$$

Therefore, the azimuth compression function is given by

$$H\_{\rm AC}(\tau, f\_{\eta}) = \exp\left[j\pi \left(\frac{4\pi R\_0 f\_0}{c} \left(D(f\_{\eta}) - 1\right) - \mathbb{C}\_0(\Lambda \tau)\right)\right] \tag{34}$$

Finally, the echo data are transformed into 2D time domain by a azimuth IFFT and the focused image is obtained.

#### **4. Experiment Results and Analysis**

In this section, we provide some imaging results to demonstrate the performance of proposed algorithm and the analysis of principle. The system parameters are listed in Table 1. The parameters of platform velocity, pulse duration and center slant range of each SAR system are set to the same value. The reference slant range is selected as the scene center slant range. The theoretical azimuth and range resolutions are evaluated by *ρa*= 0.886*c*/(4 *f*<sup>0</sup> sin(*θ*/2)) and *ρ<sup>r</sup>* = 0.886*c*/(2*Br*). The maximum Doppler frequency can be expressed as

$$f\_D = \frac{2V\_r}{\lambda\_{\rm min}} \sin\left(\frac{\theta\_{\rm max}}{2}\right) \tag{35}$$

where *λ*min and *θ*max are the minimum signal wavelength and the maximum azimuth angle. The PRF must not be chosen less than two times of *fD*. The range oversampling rate is set to 1.2.


**Table 1.** SAR system simulation parameters.

Firstly, to investigate the effects of two improvements in the proposed algorithm, a P-band SAR with a center frequency of 600 MHz was simulated. The fractional bandwidth is set to 50% (corresponding to a range resolution of 0.44 m) and the beamwidth is set to 29◦ (corresponding to an azimuth resolution of 0.44 m). Nine targets are arranged in the illuminated scene along the azimuth center at different distances form the reference range with an interval of 200 m. The conventional GCSA in Reference [26] is used in comparison with the proposed algorithm. According to the principle in Section 2.3, the data is processed with a 6th order model. To highlight the effect of Lagrange inversion, the conventional GCSA with Lagrange inversion (ignoring the range-independent coupling terms above 6th order) is also used in comparison. The 2D focused images of targets at the ranges *Rc*, *Rc* + 800 m and *Rc* + 1600 m are shown in Figure 6. In these figures, contour maps donate the 2D focusing quality. To evaluate the quality of these images quantitatively, the resolution (Res), peak sidelobe ratio (PSLR) and integrated sidelobe ratio (ISLR) along the range and azimuth directions are presented in Table 2. No wighting function or sidelobe control approach is used to obtain a fair comparison. Note that the ideal response is not completely symmetric in the range direction due to the significant range-azimuth coupling.

Figure 6a–c and Table 2 illustrate that the conventional GCSA causes deterioration of the images. The images of three targets are defocused in two directions, even at the scene center targets. The Res, PSLR and ISLR degrade considerably. This issue becomes increasingly serious as the distance increases. This problem is caused by the residual high-order range-independent coupling terms and the phase error introduced by approximate of stationary phase point. As the fractional bandwidth increases, even the scene center point has a severe degradation. Especially, the asymmetric sidelobe in range dimension is obvious. Figure 6d–f shows the imaging results obtained by the GCSA with Lagrange inversion. It can be seen that all three targets are effectively focused, which indicates that the Lagrange inversion can eliminate the range-dependent coupling terms. However, since the range-independent coupling terms above 6th order are neglected, the sidelobes in range direction are severe and the images have a certain quality degradation. The imaging results obtained by the proposed algorithm

are shown in Figure 6g–i, with better quality than Figure 6d–f. From the quality indices presented in Table 2, the imaging coherency is good over the entire swath. These benefits stem from the compensation of coupling terms above 6th order and Lagrange inversion.

**Figure 6.** Focused results by different algorithm for P-band SAR data with a center frequency of 600 MHz. (**a**–**c**) Conventional generalized chirp scaling algorithm (GCSA). (**d**–**f**) Conventional GCSA with Lagrange inversion. (**g**–**i**) Proposed algorithm. The three subgraphs of each row correspond to contours of the targets located at *Rc*, *Rc* + 800 m and *Rc* + 1600 m, respectively. The dynamic range of contour is −35 dB∼0 dB.

To evaluate the accuracy of proposed algorithm better, we measure the Res and differential resolutions (DRES) [32] at different slant range. Figure 7 shows the Res and DRES of nine point targets processed in different algorithms. The references for the DRES measurements are the range and azimuth resolutions obtained by *ω* − *k* algorithm. It can seen that the range and azimuth resolutions loss of GCSA is greater than 13% compared with the *ω* − *k* algorithm but the resolutions loss of the proposed algorithm is less than 1%. We can conclude that the focusing performances of the proposed algorithm are much better than the ones of GCSA and closed to the ones of *ω* − *k* algorithm. The image quality and focusing depth in range dimension are greatly improved.


**Table 2.** Measured parameters of the imaging results for Figure 6.

**Figure 7.** Resolutions (Res) and differential resolutions (DRES) in azimuth and range given by different algorithm, where the resolutions obtained by *ω* − *k* algorithm are references. The DRES presents the loss in spatial resolutions. Nine targets are arranged in the illuminated scene along the azimuth center at different distances form the reference range with an interval of 200 m. (**a**) Azimuth Res. (**b**) Range Res. (**c**) Azimuth DRES. (**d**) Range DRES.

Secondly, in order to better validate the performance of the proposed algorithm, the point scatterer was placed at the edge of the swath, where *R*0−*Rc* = 2 km. A typical L-band SAR system with a center frequency of 1.36 GHz was simulated. The beamwidth is set to 11◦ (corresponding to an azimuth resolution of 0.5 m). The fractional bandwidth is set to 20%, 40%, 60% and 80%, respectively. According to Equation (5), four sets of data are focused by the 3rd order, 4th order, 6th order and 7th order models, respectively.

Figure 8 shows the Res and DRES versus fractional bandwidth for the GCSA and proposed algorithm. The resolutions obtained by *ω* − *k* algorithm are references. When the fractional bandwidth is 20%, both the GCSA and proposed algorithm can achieve good focusing performance. As the fractional bandwidth increases, the performance of GCSA drops dramatically. If the resolution loss threshold is 10%, the GCSA can only process the data where the fractional bandwidth is less than 30%. However, the proposed algorithm achieves nearly the theoretically resolutions for fractional bandwidths up to 80% for L-band. The resolution broadening is less than 1%, which shows the performance is greatly improved over that of the original GCSA.

**Figure 8.** Resolutions (Res) and differential resolutions (DRES) in azimuth and range given by the generalized chirp scaling algorithm (GCSA) and proposed algorithm, where the resolutions obtained by *ω* − *k* algorithm are references. (**a**) Azimuth Res. (**b**) Range Res. (**c**) Azimuth DRES. (**d**) Range DRES.

Next, to illustrate the worst case shape of the proposed image of the point target, Figure 9 illustrates the contour plots, range profiles and azimuth profiles of the processed images at 80% fractional bandwidth. In both figures, the contour maps denote the 2D focusing quality and the profiles represent the focusing quality along the azimuth and range directions. The measured parameters are shown in Table 3. As can be seen, the results of GCSA in this case suffers form severe distortion and broadening. The signal in both range and azimuth are almost defocused. However, the proposed algorithm preserves the focusing performance. It is easy to recognize that the images obtained by the proposed algorithm are well focused, as are shown in Figure 9d–f. The nearly theoretical values of spatial resolution, PSLR and ISLR are obtained, which demonstrates the validity of the proposed algorithm. The good performance is given by the high-order range-independent phase filtering and the Lagrange inversion, which greatly reduces the range-dependent phase error. By comparing the range and azimuth profiles, it is evident that the Lagrange inversion makes the range-dependent coupling terms effectively compensated. The focusing depth in range dimension

is greatly improved. The improved algorithm is consistent with the conventional GCSA in terms of computational complexity but the focusing performance is significantly improved. The improved method provides an attractive solution for processing the low frequency large bandwidth and wide beamwidth SAR data.

**Figure 9.** Focused results of point scatterers for L-band SAR data at 80% fractional bandwidth, using the generalized chirp scaling algorithm (GCSA) and proposed algorithm. (**a**–**c**) Conventional GCSA. (**d**–**f**) Proposed algorithm. The three subgraphs of each row correspond to the contour maps, range profiles and azimuth profiles, respectively. The dynamic range of contour is −35 dB∼0 dB.


**Table 3.** Measured parameters of imaging results for Figure 9.

Finally, to further test the analysis presented in this paper, a SAR real image (Longmen, Henan, China) is used as the input radar cross section to generate SAR echo. The center frequency of the transmitted signal is 400 MHz, the bandwidth is 250 MHz, the beamwidth is 25◦, the velocity is 120 m/s, the PRF is 200 Hz, the pulse duration is 10 us and the center slant range is 10 km. The scene size is 3.0 km in range and 1.5 km in azimuth. The resolutions are 0.53 m in range and 0.76 m in azimuth. In the simulation, the input image is a real complex image. Each cell of the complex image is treated as a point scatterer (this actually contains the target signal and noise). No additional noise was added during the simulation. Figure 10 shows the imaging results processed by the conventional GCSA and the improved algorithm.

As is shown in Figure 10, at the center range region, the focusing quality of the two images is nearly the same. However, at the near and far range region, the texture character of the image obtained by improved GCSA is clearer than that obtained by conventional GCSA. It is obvious that better focusing results are obtained by the proposed method. The image has very high quality. The roads, rivers and farmland can be clearly distinguished. The edge scene along range dimension is well focused. It is evident that the range focusing depth is greatly improved.

(b)

**Figure 10.** Comparison of imaging results of real data processed by different algorithms. (**a**) Conventional generalized chirp scaling algorithm (GCSA). (**b**) Proposed algorithm.

In order to have a distinct contrast, three subregions marked by red solid rectangle are extracted and analyzed in detail. The zooms regions are shown in Figure 11a–f. The entropy [33,34] of images is calculated to compare the focus quality and shown in Table 4. It is generally acknowledged that SAR images with better quality focus have smaller entropy. It is easy to find that images obtained by proposed method have smaller entropy. Therefore, the effectiveness of proposed algorithm is again validated.

*Remote Sens.* **2019**, *11*, 1874

**Figure 11.** Zoom images. (**a**–**c**) Zoom images of subregions extracted from Figure 10a. (**d**–**f**) Zoom images of subregions extracted from Figure 10b.

**Table 4.** Image entropy of the zoom images.


#### **5. Discussion**

The performance of the proposed algorithm is mainly limited by the approximation error of range FM rate *Km* (Equation (22)). This approximation will introduce a quadratic phase error (QPE) in the spectrum and degrades the quality of the image. In the CSA, the range dependence of SRC is neglected and the FM rate is calculated at the reference range *Km*,*app* = *Kf* . The NCSA uses a linear approximation of the FM rate and the GCSA uses a 2nd order approximation model.

The QPE can be expressed as

$$\Phi\_{\rm QPE}(f\_{\rm \tau}, f\_{\rm \bar{\eta}}; \mathcal{R}\_0) = \pi f\_{\rm \tau}^2 \left( \frac{1}{K\_m} - \frac{1}{K\_{m, app}} \right) \tag{36}$$

From Equation (8), we can see that the Taylor expansion is feasible only under the following condition

$$\left| \frac{K\_r c \mathcal{R}\_0 f\_\eta^2}{2V\_r^2 f\_0^3 D \left(f\_\eta\right)^3} \right| \neq 1 \tag{37}$$

Let *<sup>G</sup>*(*Kr*, *<sup>R</sup>*0, *Vr*, *<sup>f</sup>*0, *<sup>f</sup>η*) = *KrcR*<sup>0</sup> *<sup>f</sup>* <sup>2</sup> *η* 2*V*<sup>2</sup> *<sup>r</sup> f* <sup>3</sup> <sup>0</sup> *<sup>D</sup>*(*f<sup>η</sup>* ) <sup>3</sup> , this function is positively related to *Kr*, *R*0, *f<sup>η</sup>* and negatively related to *Vr*, *f*0. And *Kr* = *Br*/*Tr*, where *Tr* is the pulse duration. The azimuth frequency varies within the following range <sup>−</sup>*PRF* <sup>2</sup> <sup>+</sup> *<sup>f</sup>η<sup>c</sup>* <sup>≤</sup> *<sup>f</sup><sup>η</sup>* <sup>≤</sup> *<sup>f</sup>η<sup>c</sup>* <sup>+</sup> *PRF* <sup>2</sup> , where PRF is the pulse repetition frequency. Thus, *G* is an even function about *fη* and inequality (37) actually implied condition *G* < 1. This inequality is satisfied in most L- and P-band SAR systems. However, as the center frequency decreases, the maximum azimuth frequency increases and the slant range increases, this inequality may not be satisfied. At this time, the Taylor approximation of range FM rate will introduce a large error and worsen the performance of algorithm.

It should be noted that *G* < 1 is required for all frequency domain approximation algorithms (such as the chirp scaling class and the range-Doppler class algorithms) [14]. This is satisfied for a linear FM signal with a large time bandwidth product (TBP). The TBP is defined as the product of the pulse duration *Tr* and bandwidth *Br*. For example, assuming the SAR system with a following parameters: center frequency *f*<sup>0</sup> = 400 MHz, bandwidth *Br* = 200 MHz, beamwidth *θ* = 29◦, velocity *Vr* = 100 m/s, pulse duration *Tr* = 2 us (TBP = 400), target slant range *R*<sup>0</sup> = 12 km and center slant range *Rc* = 10 km. Figure 12a shows the variation of *G* with azimuth frequency *fη*. Figure 12c–e show the QPEs in zero-order, 1st-order and 2nd-order approximations, respectively. It can be seen that *G* is not less than 1 at this time. Due to the existence of breakpoints, the QPEs of 1st-order and 2nd-order approximations are larger than zero-order approximation. The maximum QPE of 2nd-order model is about 5000◦, which will seriously deteriorate the image quality.

**Figure 12.** Parameters G and quadratic phase error (QPE). (**a**) Parameter G when *Tr* = 2 us (TBP = 400). (**b**) Parameter G when *Tr* = 10 us (TBP = 2000). (**c**–**e**) QPEs of zero-order, 1st-order and 2nd-order approximations corresponding to (**a**). (**f**–**h**) QPEs of zero-order, 1st-order and 2nd-order approximations corresponding to (**b**).

Change the pulse duration *Tr* = 10 us (TBP = 2000), the variation of *G* with azimuth frequency is show in Figure 12b. The QPEs of different approximations are also shown in Figure 12f–h. *G* is far less than 1, meeting the inequality (37). At this time, the maximum QPE of 1st-order model is 230◦ and the maximum QPE of 2nd-order model is 14◦, which has a limit effect on the imaging results.

By comparison, it can be concluded that the proposed algorithm has the best performance in frequency domain approximation algorithms when Equation (37) is satisfied. If Equation (37) is not satisfied, the performance of the approximation algorithm is degraded and the *ω* − *k* algorithm and time-domain algorithms are a better choice. At the same time, it should be pointed out that most of the existing airborne SAR systems have a large TBP [35,36], which shows that the proposed algorithm still has a broad application prospects. In particular, the proposed algorithm can be applied to high-resolution highly squint SAR imaging, where the phase coupling is also severe.

#### **6. Conclusions**

The high-resolution low frequency SAR has the characteristics of large bandwidth and long integration time. This trait will cause serious range-azimuth phase coupling, which limits the performance of conventional GCSA and results in image defocusing. The longer the integration time or larger bandwidth is, the more serious the deterioration will be. This paper proposes an improved GCSA based on Lagrangian inversion theorem for high-resolution low frequency SAR data processing.

Through the theoretical analysis, we find two main reasons about the defocusing. Firstly, the influence of the residual high-order phase is still significant when the fractional bandwidth is large and/or integration time is long. Secondly, the linear approximation of stationary phase point will make the high-order range-dependent phase coupling not effectively compensated. Aim to solve these two problems, this paper firstly proposes a new criterion for determining the order of Taylor expansion. The range-independent coupling phase terms above 3rd order are first compensated. Moreover, the Lagrange inversion theorem is introduced to obtain a more accurate stationary phase point. The performance and accuracy of the improved GCSA has been demonstrated using the simulated data in P- and L-band. The experimental results show that for P-band SAR systems with resolutions of 0.44 m, the proposed algorithm can focused full-swath targets with a resolution loss of less than 1%. For L-band SAR system with an azimuth resolution of 0.5 m, the edge point scatterer can be focused well even at 80% fractional bandwidth. The proposed algorithm has similar performance to the *ω* − *k* algorithm and is significantly better than the GCSA.The improved method provides the possibility to efficiently process full-swath high-resolution low frequency SAR data.

**Author Contributions:** All authors have made great contributions to the work. X.C. and T.Y. carried out the theoretical framework. X.C., T.Y., Z.H. and F.H. conceived and designed the experiments; X.C. performed the experiments and wrote the manuscript; X.C. and T.Y. analyzed the data; Z.D. gave insightful suggestions for the work and the manuscript.

**Acknowledgments:** This research was supported by the National Natural Science Found of China under Grant 61771478. The authors would like to thank all the anonymous reviewers for their valuable comments and helpful suggestions which lead to substantial improvements of this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

Here, some expressions are given for some characteristics of the SAR signal. Using the Taylor expansion principle, the coefficients *γ<sup>i</sup>* in Equation (2) can be expressed as

$$\begin{aligned} \gamma\_3 &= -\frac{D^2(f\_\eta) - 1}{2D^3(f\_\eta)f\_0^3} \\ \gamma\_4 &= -\frac{(D^2(f\_\eta) - 1)(D^2(f\_\eta) - 5)}{8D^3(f\_\eta)f\_0^4} \\ \gamma\_5 &= \frac{(D^2(f\_\eta) - 1)(3D^2(f\_\eta) - 7)}{8D^6(f\_\eta)f\_0^7} \\ \gamma\_6 &= \frac{(D^2(f\_\eta) - 1)(D^4(f\_\eta) - 14D^2(f\_\eta) + 21)}{16D^{11}(f\_\eta)f\_0^6} \end{aligned} \tag{A1}$$

### **Appendix B**

In this Appendix, we will derive the calculation of *qi*(*i* ≥ 2), *Xi*(*i* ≥ 3) and *Ci*(*i* ≥ 0).

Based on Equations (11), (17) and (20), we expand Equation (19) into a Taylor series with respect to *τs*, the Equation (21) is obtained. The coefficient *C*<sup>0</sup> is a *M*th order polynomial of Δ*τ* and it can be given by

$$\begin{aligned} \mathbb{C}\_{0}(\Delta \tau) &= \left[a^{2}q\_{2} + (a-1)^{2}K\_{\mathfrak{m}}\right] \Delta \tau^{2} + \left[a^{3}q\_{3} + (a-1)^{3}K\_{\mathfrak{m}}^{3}X\_{3}\right] \Delta \tau^{3} \\ &+ \frac{1}{8}\left[8a^{4}q\_{4} + 2(a-1)^{4}K\_{\mathfrak{m}}^{4}\left(9K\_{\mathfrak{m}}X\_{3}^{2} + 4X\_{4}\right) - 16D(f\_{\mathfrak{l}})f\_{0}(a-1)^{3}K\_{\mathfrak{m}}^{3}\gamma\_{3}\right] \Delta \tau^{4} + \dots \end{aligned} \tag{A2}$$

In order to model the rang dependence of coefficients *Ci*(*i* > 0), they are also approximated as a second order series in Δ*τ*,

$$\begin{split} \mathbf{C}\_{1} &= 2\pi \left[ (a-1)K\_{m} + aq\_{2} \right] \Delta \tau + 3\pi \left[ \alpha^{2} q\_{3} + (a-1)^{2} K\_{m}^{3} X\_{3} \right] \Delta \tau^{2} \\ \mathbf{C}\_{2} &= \pi \left( K\_{m} + q\_{2} \right) + 3\pi \left[ aq\_{3} + (a-1) K\_{m}^{3} X\_{3} \right] \Delta \tau \\ &+ \frac{3}{2} \pi \left[ 4a^{2} q\_{4} + (a-1) K\_{m}^{3} \left( 9 \left( a-1 \right) K\_{m}^{2} X\_{3}^{2} + 4 \left( a-1 \right) K\_{m} X\_{4} - 4D(f\_{\eta}) f\_{0} \gamma\_{3} \right) \right] \Delta \tau^{2} \\ \mathbf{C}\_{3} &= \pi \left( q\_{3} + K\_{m}^{3} X\_{3} \right) + \pi \left[ 4a q\_{4} + (a-1) K\_{m}^{4} \left( 9 K\_{m} X\_{3}^{2} + 4 X\_{4} \right) - 2D(f\_{\eta}) f\_{0} K\_{m}^{3} \gamma\_{3} \right] \Delta \tau \\ &+ \frac{1}{8} \pi \left[ 80 a^{2} q\_{5} + 20 (a-1)^{2} K\_{m}^{5} \left( 27 K\_{m}^{2} X\_{3}^{3} + 24 K\_{m} X\_{3} X\_{4} + 4 X\_{5} \right) \right. \\ &- 32D(f\_{\eta}) f\_{0} \left( a-1 \right) K\_{m}^{4} \left( 9 K\_{m} X\_{3} \gamma\_{3} + 2\gamma\_{4} \right) \right] \Delta \tau^{2} \\ &- \dots \end{split}$$

For each *Ci*, the higher order terms of Δ*τ* are very small and can be neglected. Combining Equation (22) and (A3), the coefficients *Ci*(*i* > 0) can be rewritten as

$$\begin{aligned} \mathbf{C}\_{1} &= \mathbf{C}\_{11}\Delta\tau + \mathbf{C}\_{12}\Delta\tau^{2} \\ \mathbf{C}\_{2} &= \mathbf{C}\_{20} + \mathbf{C}\_{21}\Delta\tau + \mathbf{C}\_{22}\Delta\tau^{2} \\ \mathbf{C}\_{3} &= \mathbf{C}\_{30} + \mathbf{C}\_{31}\Delta\tau + \mathbf{C}\_{32}\Delta\tau^{2} \\ &\dots \\ \mathbf{C}\_{M} &= \mathbf{C}\_{M0} + \mathbf{C}\_{M1}\Delta\tau + \mathbf{C}\_{M2}\Delta\tau^{2} \end{aligned} \tag{A4}$$

with

*C*<sup>11</sup> = 2 *αq*<sup>2</sup> + *Kf* (*α* − 1) *C*<sup>12</sup> = 3*α*2*q*<sup>3</sup> + *K*<sup>2</sup> *<sup>f</sup>* (*α* − 1) 2*Ks* + 3*Kf* (*α* − 1) *X*<sup>3</sup> *C*<sup>20</sup> = *Kf* + *q*<sup>2</sup> *C*<sup>21</sup> = 3*αq*<sup>3</sup> + *K*<sup>2</sup> *f Ks* + 3*Kf* (*α* − 1) *X*<sup>3</sup> *C*<sup>22</sup> = <sup>1</sup> 2 12*α*2*q*<sup>4</sup> + *K*<sup>3</sup> *<sup>f</sup>*(18*Kf Ks* (*<sup>α</sup>* − <sup>1</sup>) *<sup>X</sup>*<sup>3</sup> + <sup>27</sup>*K*<sup>2</sup> *<sup>f</sup>*(*α* − 1) 2 *X*2 <sup>3</sup>+ <sup>2</sup>(*K*<sup>2</sup> *<sup>s</sup>* + 6*Kf*(*α* − 1) 2 *X*4) −12*D*(*fη*)*f*<sup>0</sup> (*α* − 1) *γ*3) *C*<sup>30</sup> = *q*<sup>3</sup> + *K*<sup>3</sup> *<sup>f</sup> X*<sup>3</sup> *C*<sup>31</sup> = 4*αq*<sup>4</sup> + 3*K*<sup>4</sup> *<sup>f</sup> KsX*<sup>3</sup> + *<sup>K</sup>*<sup>4</sup> *<sup>f</sup>* (*α* − 1) 9*Kf X*<sup>2</sup> <sup>3</sup> + 4*X*<sup>4</sup> − <sup>2</sup>*D*(*fη*)*f*0*K*<sup>3</sup> *<sup>f</sup> γ*<sup>3</sup> *C*<sup>32</sup> = <sup>1</sup> 8 48*K*<sup>5</sup> *f K*2 *<sup>s</sup> X*<sup>3</sup> + 20*K*<sup>5</sup> *<sup>f</sup>*(*α* − 1) 2 27*K*<sup>2</sup> *<sup>f</sup> <sup>X</sup>*<sup>3</sup> <sup>3</sup> + 24*Kf X*3*X*<sup>4</sup> + 4*X*<sup>5</sup> − <sup>48</sup>*K*<sup>4</sup> *<sup>f</sup> Ks D*(*fη*)*f*0*γ*<sup>3</sup> +80*α*2*q*<sup>5</sup> + 8*K*<sup>5</sup> *<sup>f</sup> Ks* (*α* − 1) 45*Kf X*<sup>2</sup> <sup>3</sup> + 16*X*<sup>4</sup> − <sup>32</sup>*D*(*fη*)*f*0*K*<sup>4</sup> *<sup>f</sup>* (*α* − 1) 9*Kf X*3*γ*<sup>3</sup> + 2*γ*<sup>4</sup> ... (A5)

According to (A4), the coefficients *Ci* are range-dependent. The expression for these coefficients

contains linear and quadratic terms in Δ*τ*. To remove the range dependence of these terms, it is required to set the coefficients of Δ*τ* to zero. Let *C*<sup>11</sup> in (A5) be zero, we can obtain

$$q\_2 = K\_f \frac{1 - \alpha}{\alpha} \tag{A6}$$

Let *C*<sup>12</sup> and *C*<sup>21</sup> be zero, we can obtain

$$\begin{array}{l} q\_3 = \frac{k\_f^2 K\_\circ (1 - a)}{\lambda} \\ X\_3 = \frac{K\_\circ (a - 2)}{3 K\_f (a - 1)} \end{array} \tag{A7}$$

Similarly, *CM*−1,2 and *CM*<sup>1</sup> be zero, the expressions of *qM* and *XM* can be obtained. And the coefficients *Ci* becomes

$$\begin{aligned} \mathbf{C}\_1 &= 0\\ \mathbf{C}\_2 &= \mathbf{C}\_{20} \\ \mathbf{C}\_3 &= \mathbf{C}\_{30} \\ \dots \\ \mathbf{C}\_M &= \mathbf{C}\_{M0} \end{aligned} \tag{A8}$$

Thus, Equations (25) and (26) are obtained.

#### **Appendix C**

This Appendix is used to explain how to use the Lagrange inversion theorem to solve the stationary phase point. According to Equations (12) and (13), *τ* can be expressed as a series form of *fτ*. A Lagrange inversion expression for a general power series is given here. For a function expressed in a series

$$
\omega = h\left(z\right) = a\_0 + a\_1\left(z - z\_0\right) + \dots + a\_{\text{ll}}\left(z - z\_0\right)^{\text{n}}\tag{A9}
$$

$$f\left(z\right) = \frac{z - z\_0}{h\left(z\right) - a\_0} = \frac{1}{a\_1 + a\_2 \left(z - z\_0\right) + \dots + a\_{\text{ll}} \left(z - z\_0\right)^{n-1}}\tag{A10}$$

Thus, we can get

$$g\_1 = f(z)|\_{z = z\_0} = \frac{1}{a\_1} \tag{A11}$$

$$\text{g}\_2 = \frac{1}{2} \frac{df^2(z)}{dz} \Big|\_{z=z\_0} = -\frac{a\_2}{a\_1^3} \tag{A12}$$

$$\log\_3 = \frac{1}{6} \frac{d^2 f^3(z)}{dz^2} \Big|\_{z=z\_0} = \frac{2a\_2^3}{a\_1^5} - \frac{a\_3}{a\_1^4} \tag{A13}$$

$$\log\_4 = \frac{1}{24} \frac{d^3 f^4(z)}{dz^3} \Big|\_{z=z\_0} = -\frac{5a\_2 \,^3 - 5a\_1 a\_2 a\_3 + a\_1 \,^2 a\_4}{a\_1 \overline{\,}} \tag{A14}$$

$$\mathcal{g}\_5 = \frac{1}{120} \frac{d^4 f^5(z)}{dz^4} \Big|\_{z=z\_0} = \frac{14a\_2^4 - 21a\_1a\_2^2a\_3 + 3a\_1^2a\_3^2 + 6a\_1^2a\_2a\_4 - a\_1^3a\_5}{a\_1^9} \tag{A15}$$

$$g\_6 = \frac{1}{720} \frac{d^5 f^6(z)}{dz^5} \Big|\_{z=z\_0} = \frac{1}{a\_1^{17}} \left( -42a\_2^5 + 84a\_1 a\_2^3 a\_3 - 28a\_1^2 a\_2^2 a\_4 + 7a\_1^2 a\_2 \left( -4a\_3^2 + a\_1 a\_5 \right) + a\_1^3 \left( 7a\_3 a\_4 - a\_1 a\_6 \right) \right) \tag{A16}$$

Equations (12) and (27) are special forms of Equation (A9) and the inversion expressions can be easily obtained according to Equation (14).
