Squint Mode GEO SAR Imaging Using Bulk Range Walk Correction on Received Signals

Abstract

Geosynchronous synthetic aperture radar (GEO SAR) has the potential for conducting long-term observation of target zones, which is essential for remote sensing applications such as disaster monitoring and vegetation measurements. The squint imaging mode is crucial for long-term observation using GEO SAR. However, this type of SAR imaging is problematic because the squint mode introduces a nonzero range cell walk, which increases the prevalence of invalid data in echoes and intensifies the coupling between the azimuth and range. Therefore, this paper proposes a novel squint mode GEO SAR imaging method based on the correction of the bulk range walk of received signals. Adjusting the starting time of the receiving window significantly reduces the redundancy in echoes. Then, first-order filtering, range cell migration correction, range compression, partial dechirp, and azimuth compression are used to obtain the imaging result. Simulation results for the GEO SAR imaging of Wenchuan County in China demonstrate that the proposed algorithm can achieve a resolution of 5 m within a 30 × 30 km swath over 48% of the orbital period.

1. Introduction

Geosynchronous synthetic aperture radar (GEO SAR) operates on an inclined orbit at a speed of nearly 3.075 km/s and at a height of approximately 35,786 km [1,2]. Compared with low earth orbit (LEO) SAR, GEO SAR moves more slowly and has a wider field of view that is nearly one-third of the earth’s surface [3]. These features enable GEO SAR to observe target zones for a much longer period than LEO SAR, which is beneficial for conducting disaster monitoring [4]. Figure 1 illustrates an observation mode suitable for long-term observation via GEO SAR. This mode is composed of several strip-map observations acquired with different squint angles, and the target zone is imaged according to the overlapping area of the strip-map data received. Here, not only broadside SAR imaging, but also squint mode SAR imaging, must be used to process echoes acquired by this observation mode. The present study focuses on squint mode strip-map GEO SAR imaging.

The relative trajectories between GEO SAR platforms and targets are curved, which varies the range histories spatially in both the azimuth and range directions [5]. A number of GEO SAR imaging algorithms have been developed to correct this two-dimensional (2D) spatial variance. For example, the time-domain back projection (BP) algorithm has been developed for image focusing [6]. However, this algorithm suffers from its tremendous computational burden, even though the fast BP algorithm reduces the computational burden to O(N^5/2) [7]. In contrast, the computational complexity of frequency-domain algorithms is much lower, and numerous algorithms that operate in this domain have been developed. The classic frequency-domain SAR imaging algorithms include the omega-K methods presented by Claudio Prati and Fabio Rocca [8], the methods based on the chirp scaling [9], and the frequency domain algorithms proposed by Franceschetti et al. [10]. As aforementioned, the severe spatial variance makes these classical algorithms fail in GEO SAR. However, inspired by these methods, some scholars developed algorithms for GEO SAR. Bao et al. [11] modified the chirp scaling (CS) algorithm based on a fourth-order polynomial model, and achieved an 8-m resolution within a 40 × 40 km swath. Hu et al. [12] improved the conventional nonlinear CS algorithm by accounting for a linear spatial variance in the range direction in the polynomial model. Tao et al. [13] developed the nonlinear azimuth CS algorithm by introducing spatially variant components in the classic hyperbolic model. Sun et al. [14] considered spatial variation in both the range and azimuth directions in a fifth-order polynomial model, and a 5-m resolution within an 82 × 95 km swath was achieved using range cell migration (RCM) equalization and sub-band synthesis technologies. To avoid the requirement for conducting sub-band synthesis, Chen et al. [15] combined singular value decomposition and azimuth nonlinear scaling, and obtained a 20-m resolution within a 150 × 130 km swath. Yu et al. [16] proposed a time-frequency scaling algorithm to correct linear and quadratic variances in the azimuth direction and obtained well-focused results. Ding et al. [17] modified the quadratic term in the azimuth direction and obtained a 20-m resolution within a 400 × 200 km swath under conditions where the Doppler rate was zero or approximately zero.

Notably, all of the aforementioned frequency-domain algorithms were designed for broadside GEO SAR imaging where the range cell walk is zero. However, the squint mode introduces a nonzero range walk, which increases the prevalence of invalid data in echoes, and intensifies the coupling between the azimuth and the range. As a result, squint mode SAR imaging is considerably more complicated than broadside imaging. Therefore, correcting the range walk is the key to developing high-resolution squint mode GEO SAR imaging. Whereas past efforts have addressed this issue using range migration removal and 2D nonlinear CS [18], the researchers assumed that the range walk was spatially invariant over the entire swath, which restricted the performance of the developed algorithm. Therefore, we propose a novel squint mode GEO SAR imaging method based on the correction of the bulk range walk of echoes. The key component of the method is to exploit the characteristics of the low pulse repetition frequency (PRF) employed in GEO SAR, and accordingly adjust the starting time of the receiving window. Then, first-order filtering, RCM correction (RCMC), range compression, partial dechirp, and azimuth compression are used to obtain the imaging result. Although the bulk range walk can be corrected in the signal processing stage, the proposed method conducts the correction while receiving the echoes (i.e., during receiving), which significantly reduces the echoes redundancy and improve the efficiency of on-board storage.

The remainder of this paper is structured as follows. Section 2 analyzes the range walk of squint mode GEO SAR based on a slant range model, and illustrates the advantages of conducting range walk removal upon receiving. Section 3 discusses the feasibility of conducting range walk removal on receiving by comparing the Doppler bandwidths of GEO SAR and LEO SAR, and presents the proposed approach for adjusting the starting time of the receive window. An imaging algorithm based on echoes with the bulk range walk corrected is then proposed in Section 4. Simulation results for the squint mode GEO SAR imaging of Wenchuan County in China are presented and analyzed in Section 5, and the results are discussed in Section 6. Finally, Section 7 concludes the paper.

2. Slant Range History of Squint Mode GEO SAR

Figure 2 illustrates the observation geometry between the strip-map GEO SAR imaging and a target zone as the GEO SAR platform moves along its orbit. At azimuth time t_ref, the beam center points to the reference target P_ref, which is usually the center of the target zone, and the slant range between the beam center and P_ref is denoted by R_ref. At the azimuth time t_c, the beam center points to P_c, and the corresponding slant distance is R_c. The slant range history of GEO SAR imaging can be analyzed according to Figure 2 as follows.

Ignoring the antenna weighting, the GEO SAR echo from P_c can be expressed as:

$$S_{echo}\left(t,\tau \right)=\sigma \cdot rect\left(\frac{t-t_c}{T_s}\right)\cdot rect\left[\frac{\tau -2R\left(t\right)/c}{T_p}\right]\cdot \exp \left[j\pi K_r{\left(\tau -\frac{2R\left(t\right)}{c}\right)}^2\right]\cdot \exp \left[-j\frac{4\pi }{\lambda }R\left(t\right)\right],$$

(1)

where t and $\tau$ denote the slow time in the azimuth direction and the fast time in the range direction, respectively; $\sigma$ is the radar cross-section of P_c; rect(·) is the rectangle function; T_s is the synthetic aperture time; R(t) represents the slant range history between the beam center and P_c; c is the speed of light; T_p is the pulse width of the transmitted chirp signal; K_r is the frequency modulation rate; and $\lambda$ is the wavelength of the beam. As previously reported [16], R(t) can be described by a polynomial expression of order N, as follows:

$$R\left(t\right)\approx {\displaystyle \sum _{n=0}^N\frac{r_n}{n!}{\left(t-t_c\right)}^n},$$

(2)

where r_n is the nth polynomial coefficient that can be obtained by fitting based on ephemeris data and geographic information specific to the target zone. A value of N = 5 has been demonstrated to be sufficiently large for GEO SAR imaging [16], and that value is assumed hereafter. However, because R(t) varies for different target zones at different locations, the coefficients r_n are spatially variant, and can be expressed as functions of ${\Delta R=R}_c-R_{ref}$ and $\Delta {t=t}_c-t_{ref}$, as follows:

$$\{\begin{array}{l}{r}_0=r_{0,ref}+\Delta R\\ r_1=r_{1,ref}+k_{1,1,r}\cdot \Delta R+{k_{1,1,a}|}_{\Delta R}\cdot \Delta t\\ r_2=r_{2,ref}+k_{2,1,r}\cdot \Delta R+k_{2,2,r}\cdot \Delta R^2+{k_{2,1,a}|}_{\Delta R}\cdot \Delta t+{k_{2,2,a}|}_{\Delta R}\cdot \Delta t^2\\ r_3=r_{3,ref}+k_{3,1,r}\cdot \Delta R+{k_{3,1,a}|}_{\Delta R}\cdot \Delta t+{k_{3,2,a}|}_{\Delta R}\cdot \Delta t^2\\ r_4=r_{4,ref}+k_{4,1,r}\cdot \Delta R+{k_{4,1,a}|}_{\Delta R}\cdot \Delta t\\ r_5=r_{5,ref}\end{array}$$

(3)

where the auxiliary coefficients r_n,ref, $k_{n,m,r}$, and ${k_{n,m,a}|}_{\Delta R}$, where ${k_{n,m,a}|}_{\Delta R}$ denotes the spatially variant coefficient $k_{n,m,r}$ at $\Delta R$, can be calculated as presented previously in Yu et al. [16].

As Figure 3 shows, echoes usually contain range walk ${\Delta R}_W$ and range curvature ${\Delta R}_B$. According to Equation (2), the first-order coefficient r₁ mainly causes the range walk. Therefore, the range walk and the range curvature are given as follows:

$$\{\begin{array}{l}\Delta R_W=r_1{T}_s\approx r_{1,ref}{T}_s\\ \Delta R_B=\underset{t\in T_s}{\max }\left\{R\left(t\right)-\Delta R_W-R_c\right\}\approx r_{2,ref}{T}_s^2/8\end{array}$$

(4)

In frequency domain squint mode GEO SAR imaging, both $\Delta R_W$ and $\Delta R_B$ should be corrected to make the migration curve limited in one single range bin.

To evaluate the proportion of $\Delta R_W$ requiring correction, we define a ratio ${\chi }_W$ as:

$${\chi }_W=\frac{\Delta R_W}{\Delta R_W+\Delta R_B}\approx \frac{8r_{1,ref}}{8r_{1,ref}+r_{2,ref}{T}_s}.$$

(5)

Figure 4 illustrates the value of ${\chi }_W$ obtained for different squint angles ${\theta }_s$, where the GEO SAR platform operates on an orbit with an inclination angle of 60°. Here, we note that ${\chi }_W$ generally increases with increasing ${\theta }_s$. ${\chi }_W$ is greater than 90% for the most of the orbital period. This indicates that $\Delta R_W$ accounts for the vast majority of the correction required. We note that for a swath of width $\Delta W$, the volume of recorded echo data can be given as:

$$V_{echo}\approx \frac{2}{c}\left(\Delta W\cdot \sin {\theta }_{inc}+\frac{c{T}_p}{2}+\max \left\{\Delta R_B,\Delta R_W\right\}\right)\cdot F_s\cdot N_a,$$

(6)

where ${\theta }_{inc}$ is the incidence angle corresponding to the center of the target zone, $F_s$ is the range sampling rate, and $N_a$ is the azimuth sampling rate. From Equation (5), the ratio $\Delta R_W∕\Delta R_B$ is given by ${\chi }_W∕\left(1-{\chi }_W\right)$, which means $\Delta R_W$ is much higher than 9 times of $\Delta R_B$. Equation (6) indicates that V_echo is significantly reduced if $\Delta R_W$ is corrected on receiving rather than storing the data for later image processing. Consequently, the efficiency of on-board GEO SAR data storage increases.

3. Bulk Range Walk Correction on Receive

Compared with LEO SAR, the Doppler band of GEO SAR is much narrower. This feature makes it feasible to conduct bulk range walk correction by adjusting the starting time of the receiving window. The specific features of GEO SAR and the proposed adjustment principle are presented in the following discussion.

3.1. Doppler Bandwidth of GEO SAR

Since the relative trajectory between a GEO SAR platform and a target zone is curved, the azimuth resolution is no longer inversely proportional to the Doppler bandwidth. The relationship between the Doppler bandwidth and azimuth resolution is established in the wavenumber domain in this paper.

As illustrated in Figure 5, a GEO SAR platform moves from S₁ to S₂ while observing P_c. S₀ is the midpoint of the arc $\widehat{S_1{S}_2}$. The corresponding Doppler bandwidth is:

$$B_a=|f_{d2}-f_{d1}|,$$

(7)

where f_d1 and f_d2 are the instantaneous Doppler frequencies at S₁ and S₂, respectively, which are given as follows:

$$f_{d1}=\left(2/\lambda \right)\cdot {\overrightarrow{V}}_{st1}\cdot \left({\overrightarrow{R}}_t-{\overrightarrow{R}}_{s1}\right)/|{\overrightarrow{R}}_t-{\overrightarrow{R}}_{s1}|$$

(7a)

$$f_{d2}=\left(2/\lambda \right)\cdot {\overrightarrow{V}}_{st2}\cdot \left({\overrightarrow{R}}_t-{\overrightarrow{R}}_{s2}\right)/|{\overrightarrow{R}}_t-{\overrightarrow{R}}_{s2}|$$

(7b)

where $\overrightarrow{V_{st1}}$ and $\overrightarrow{V_{st2}}$ are the respective relative velocity vectors between the GEO SAR platform and P_c at S₁ and S₂, $\overrightarrow{R_{s1}}$ and $\overrightarrow{R_{s2}}$ are the respective platform position vectors at S₁ and S₂, and $\overrightarrow{R_t}$ is the position vector of P_c. In Figure 5, $\overrightarrow{i_{st0}}$, $\overrightarrow{i_{st1}}$, and $\overrightarrow{i_{st2}}$ are the unit vectors from P_c. to S₀, S₁, and S₂, respectively. The wavenumber vectors corresponding to $\overrightarrow{i_{st0}}$, $\overrightarrow{i_{st1}}$, and $\overrightarrow{i_{st2}}$ are given as follows, respectively:

$$\overrightarrow{k_{r0}}=\frac{4\pi }{\lambda }\overrightarrow{i_{st0}},\overrightarrow{k_{r1}}=\frac{4\pi }{\lambda }\overrightarrow{i_{st1}},\mathrm{and}\overrightarrow{k_{r2}}=\frac{4\pi }{\lambda }\overrightarrow{i_{st2}}.$$

(8)

In the wavenumber domain, the azimuth resolution ${\rho }_a$ is determined by the azimuth wavenumber vector $\overrightarrow{k_{gcr}}$, which is represented by the red arrow in Figure 5. Here, $\overrightarrow{k_{gcr}}$ is the projection of the vector $\overrightarrow{k_{r12}}=\overrightarrow{k_{r2}}-\overrightarrow{k_{r1}}$ in the azimuth direction. Therefore, $\overrightarrow{k_{gcr}}$ can be given as:

$$\overrightarrow{k_{gcr}}=\overrightarrow{k_{r12}}\cdot \sin \gamma ,$$

(9)

where $\gamma$ is the angle between $\overrightarrow{k_{r0}}$ and $\overrightarrow{k_{r12}}$ and is defined as:

$$\gamma =\arcsin \left\{\sqrt{1-{\left[\overrightarrow{i_{st0}}\cdot \left(\overrightarrow{i_{st2}}-\overrightarrow{i_{st1}}\right)\right]}^2}\right\}.$$

Accordingly, the azimuth resolution can be calculated as:

$${\rho }_a=\frac{2\pi }{k_{gcr}}=\frac{\lambda }{2|\overrightarrow{i_{st2}}-\overrightarrow{i_{st1}}|\cdot \sin \gamma }.$$

(10)

The relationship between azimuth resolution and Doppler bandwidths varies with different targets at different orbit positions. Equations (7) and (10) provide a solution. Using the iterative method to determine the starting observation position (e.g., S₁ in Figure 5) and ending observing position (e.g., S₂ in Figure 5) under the predetermined azimuth resolution, the Doppler bandwidths can be obtained by Equation (7). Figure 6 presents the Doppler bandwidths of GEO SAR and LEO SAR with an azimuth resolution of 5 m. The figure indicates that the Doppler bandwidth of GEO SAR is much smaller than that of LEO SAR, which is caused by the slower relative speed between GEO SAR platforms and targets.

3.2. Adjusting the Starting Time of the Receive Window

Spaceborne SAR transmits signals at the PRF, which is usually 1.1–1.2 times the Doppler bandwidth [20]. The echo is recorded between two adjacent transmission pulses in the receive window. The starting time of the receive window determines the nearest range at which the beginning of the echo can be recorded. Therefore, adopting a variable starting time can alter the form of the recorded data. This is illustrated in Figure 7 for the echo data obtained from a single isolated point target when employing a constant starting time given by the blue line in Figure 7a, and when employing the variable starting time, given by the red line in Figure 7b. We note that adjusting the starting time greatly reduces the RCM corresponding to the target, which indicates that the echoes can be recorded using less on-board storage and the range-azimuth coupling decreased. To avoid overlapping between the transmitted pulses and the received signals, the scope for adjusting the starting time must satisfy

$$\Delta t_a<PRT-T_p,$$

(11)

where the pulse repetitive time is PRT = 1/PRF. For LEO SAR with a resolution of 5 m, $\Delta t_a$ is less than 0.7 ms [21], which provides no scope for adjusting the starting time. However, the much lower PRF employed with GEO SAR relative to that employed with LEO SAR provides a value of $\Delta t_a$ as long as 12 ms, which provides ample scope for adjusting the starting time of the receive window.

In this paper, we adjusted the starting time to correct the bulk range walk of the reference target, which is more than 90% of the sum of $\Delta R_W$ and $\Delta R_B$ for an entire swath. Since the range walk is directly defined by r₁ in Equation (4), the range walk corrected means this item is removed. Therefore, the slant range of a reference target after correction becomes:

$$R_{ref\_RWC}\left(t\right)=r_{0\_ref}+{\displaystyle \sum _{n=2}^5\frac{r_n^n}{n!}{t}^n}.$$

(12)

By comparing Equations (2) and (12), we note that the linear component representing the range migration disappears. Therefore, the adjusted starting time $T_{S\_RWC}$ should satisfy:

$$T_{s\_RWC}=T_{s\_general}+2\left[R_{ref\_RWC}\left(t\right)-R_{ref\_act}\left(t\right)\right]/c,$$

(13)

where T_{s_general} is the starting time before adjustment and R_{ref_RWC}(t) is the actual range history of the reference target. Furthermore, correcting the bulk range walk on receive ensures that r_1,ref in Equation (3) is equal to 0, which indicates that:

$$r_1=k_{1,1,r}\cdot \Delta R+{k_{1,1,a}|}_{\Delta R}\cdot \Delta t.$$

(14)

Equation (13) can be further explained based on the observation geometry illustrated in Figure 8. Here, GEO SAR observes the reference target P_ref with a squint angle ${\theta }_s$. The echoes are sampled in fast and slow time. Here, fast time occurs in the beam pointing direction, which is also the direction of transmission and reception, and slow time occurs in the flight direction, i.e., the azimuth direction, prior to adjusting the starting time of the receive window. By applying Equation (13), the receive time is postponed from A to A′ or advanced from B to B′. As a result, slow time is changed from the flight direction to the new direction denoted by the red dotted line. As such, the deviation between the red dotted line and the flight trajectory complies with Equation (13). The operation makes the new slow time and fast time directions orthogonal, and changes the squint mode to an equivalent broadside mode. However, this change from squint mode to broadside mode is only applicable to P_ref, and the other targets in the swath are still subject to a slight degree of squint mode observation.

4. Squint Imaging Algorithm for GEO SAR after Correcting Bulk Range Walk

The proposed imaging process comprising operations in the range and azimuth directions after correcting the bulk range walk on receive is illustrated in Figure 9. The echo data of four isolated point targets with bulk range walk removed is shown in Figure 9a. First-order filtering was applied in the range direction to compensate for the range-dependent term $k_{1,1,r}\cdot \Delta R$ in Equation (14), as shown in Figure 9b. Then, RCMC and range compression were implemented to obtain the signal illustrated in Figure 9c. Next, partial dechirp was used to shorten the span of the aperture time for every target to facilitate the azimuthal segmentation of echoes. The results are shown in Figure 9d. Then, azimuth compensation was conducted, which led to the imaging result illustrated in Figure 9e. The proposed squint mode imaging algorithm is described by the flowchart given in Figure 10. The detailed steps in the algorithm are outlined in the following subsections.

4.1. Operations in the Range Direction

After transforming Equation (1) into the range frequency domain, the return signal is given as:

$$S_{rngF}\left(t,f_{\tau };t_c,\Delta R\right)=\sigma \cdot rect\left(\frac{t-t_c}{T_s}\right)\cdot rect\left(\frac{f_{\tau }}{B_r}\right)\cdot \exp \left(-j\pi \frac{f_{\tau }^2}{K_r}\right)\cdot \exp \left[-j\frac{4\pi \left(f_{\tau }+f_0\right)}{c}R\left(t\right)\right],$$

(15)

where $f_{\tau }$ is the range frequency, $B_r=K_r{T}_p$ is the bandwidth of the transmitted signal, and $f_0$ is the carrier frequency. To increase the accuracy of RCMC, the term $k_{1,1,r}\cdot \Delta R$ in Equation (14) should first be removed. Because $k_{1,1,r}\cdot \Delta R$ depends on the position of a target, the return signal must be segmented in the range direction. The length of each segment $\Delta R_{seg}$ can be 1.5–2 times $\Delta R_B$. Then, the range-dependent term $k_{1,1,r}\cdot \Delta R$ can be compensated in each segment by the following first-order filter:

$$H_{Lrcmc}=\exp \left\{j\frac{4\pi \left(f_0+f_{\tau }\right)}{c}\cdot k_{1,1,r}\cdot \Delta R_{seg}\cdot t\right\}.$$

(16)

Subsequently, RCMC and range compression are implemented using the fourth-order nonlinear chirp scaling imaging algorithm [22], which extends the RCMC ability in the range direction, and the echo becomes (Appendix A):

$$S\left(t,\tau ;t_c,\Delta R\right)=\sigma \cdot rect\left(\frac{t-t_c}{T_s}\right)\cdot \sin c\left(\tau -{\tau }_{mig}\right)\cdot \exp \left\{-j\frac{4\pi }{\lambda }R\left(t\right)\right\},$$

(17)

where ${\tau }_{mig}=2R_c/c-T_{s\_RWC}$. In addition, R(t) has the same form as that given in Equation (2), where $r_1=k_{1,1,a}\cdot \Delta t_c$, and all other parameters, including r₀, r₂, r₃, r₄, and r₅, are the same as those given in Equation (3).

4.2. Operation in the Azimuth Direction

Because the signal has been well compressed in the range direction, the azimuth operations can be implemented in each range gate. The first time-scaling given as:

$$H_{AS1}\left(t;\Delta R\right)=\exp \left\{j\frac{2\pi }{\lambda }{k_{1,1,a}|}_{\Delta R}{t}^2\right\}$$

(18)

is used to remove the azimuth-dependent term ${k_{1,1,a}|}_{\Delta R}\cdot \Delta t_c$. Then, the frequency-scaling

$$H_{Y_3{Y}_4}(f_t;\Delta R)=\exp \left\{j\frac{2\pi {Y_3|}_{\Delta R}}{3\lambda }{\left(\frac{\lambda }{2}\right)}^3{f}_t^3-j\frac{\pi {Y_4|}_{\Delta R}}{\lambda }{\left(\frac{\lambda }{2}\right)}^4{f}_t^4\right\}$$

(19)

and the second time-scaling

$$H_{AS2}(t)=\exp \left\{-j\frac{4\pi }{\lambda }\cdot {\displaystyle \sum _{n=3}^5\frac{p_n}{n!}}{t}^n\right\}$$

(20)

are subsequently applied to correct the quadratic and linear spatial variances in r_n (n = 3, 4, 5), respectively. Expressions for the undefined parameters in Equations (19) and (20), including Y₃, Y₄, p₃, p₄, and p₅, are given in Appendix B. After removing the azimuth spatial variance in r_n (n = 3, 4, 5), the higher-order items in the frequency domain can be compensated by applying

$$H_{apc}\left(f_t\right)=\exp \left\{-j\frac{4\pi }{\lambda }{\displaystyle \sum _{m=3}^{10}{P}_m}{\left(\frac{\lambda f_t}{2}\right)}^m\right\},$$

(21)

where the values of P_m (m = 3–10) are defined in Appendix C. Implementing Equations (18)–(21) transforms the return signal as follows:

$$\begin{array}{c}{S}_{apc}\left(t,\tau ;t_c,\Delta R\right)\approx \sigma \cdot rect\left(\frac{t-t_c}{T_s}\right)\cdot \sin c\left(\tau -{\tau }_{mig}\right)\hfill \\ \cdot \exp \left\{-j\frac{4\pi }{\lambda }\left[r_0-\frac{{k_{1,1,a}|}_{\Delta R}}{2}{t}_c^2+\frac{1}{2}\left(r_2-{k_{1,1,a}|}_{\Delta R}\right){\left(t-t_c\right)}^2+{\displaystyle \sum _{n=3}^5\frac{{r_n|}_{\Delta t_c=0}}{n}{\left(t-t_c\right)}^n}\right]\right\}\hfill \end{array}$$

(22)

The time-frequency relationships corresponding to Equations (17) and (22) are illustrated in Figure 11a,b, respectively. Here, the Doppler bandwidth is transformed from $B_{a\_origin}\approx 2|r_2|T_s/\lambda$ to $B_{a\_AS}\approx 2|r_2-{k_{1,1,a}|}_{\Delta R}|T_s/\lambda$. However, the Doppler bandwidth must be recovered to preserve the azimuth resolution. This is implemented using partial dechirp, which consists of applying the following two steps:

$$\{\begin{array}{l}{H}_{PD\_1}\left(f_t\right)=\exp \left\{-j2\pi \cdot \chi \cdot {\varphi }_0\left(f_t\right)\right\}\\ H_{PD\_2}\left(t\right)=\exp \left\{-j\frac{4\pi }{\lambda }\cdot \frac{{k_{1,1,a}|}_{\Delta R}}{2\left(1-\chi \right)}{t}^2\right\}\end{array}$$

(23)

Here, ${\varphi }_0\left(f_t\right)$ can be calculated according to Equation (A8) in Appendix A, and $\chi$ is the partial dechirp coefficient, where $0<\chi <1$. The application of Equation (23) transforms the echo as follows:

$$\begin{array}{l}{S}_{PD}\left(t,\tau ;t_c,\Delta R\right)=\sigma \cdot rect\left[\frac{t-t_c}{\left(1-\chi \right)T_s}\right]\cdot \sin c\left(\tau -{\tau }_{mig}\right)\\ \cdot \exp \left\{-j\frac{4\pi }{\lambda }\left[{r^{\prime }}_0+{r^{\prime }}_1\cdot \left(t-t_c\right)+\frac{{r^{\prime }}_2}{2}\cdot {\left(t-t_c\right)}^2+{\displaystyle \sum _{n=3}^5\frac{{r_{n,ref}|}_{\Delta R}}{n}{\left(t-t_c\right)}^n}\right]\right\}\end{array}$$

(24)

where $r_0^{\prime }=r_0+{k_{1,1,a}|}_{\Delta R}\cdot \left[1/\left(1-\chi \right)-1/2\right]$, $r_1^{\prime }=2{k_{1,1,a}|}_{\Delta R}\cdot t_c/\left(1-\chi \right)$, and $r_2^{\prime }=r_2/\left(1-\chi \right)$. The azimuth time-frequency relationship corresponding to Equation (24) is shown in Figure 11c. The Doppler bandwidth is recovered and the span of the aperture time is reduced to $T_{s\_pd}=\left(1-\chi \right)T_s$.

After conducting partial dechirp, azimuthal segmentation is employed, where the width of each segment is $\left(1-\chi \right)T_s$. We can compensate for the azimuth-dependent term in Equation (24) by applying

$$H_{subScene}\left(t;t_{c,sub}\right)=\exp \left\{j\frac{4\pi }{\lambda }\frac{{k_{1,1,a}|}_{\Delta R}{t}_{c,sub}}{\left(1-\chi \right)}\left(t-t_{c,sub}\right)\right\}$$

(25)

to each azimuth segment, where $t_{c,sub}$ denotes the center time of each azimuth segment. Once the spatial variance in the azimuth direction has been removed, azimuth compression is implemented by applying:

$$H_{aziCompress}\left(f_t\right)=\exp \left\{-j\frac{\pi \lambda \left(1-\chi \right)}{{2r_2|}_{\Delta t_c=0}}{f}_t^2\right\}.$$

(26)

Then, the final imaging result is obtained by transforming the return signal into the 2D time domain.

5. Simulation and Analysis

5.1. Simulation Method and Parameters

Since the flight speed of GEO SAR is nearly equivalent to the earth’s rotation, the relative speed between targets and the satellite drastically vary with the satellite’s position, termed the argument of latitude (AOL). This varies the slant range history in different azimuth time, i.e., the Doppler characteristic is spatially variant. Only one tailored frequency-domain imaging algorithm cannot focus well in the full orbit time. In order to evaluate the feasibility over the full orbit time, the spatial variance condition of the echoes should be quantized. Measuring the variance involves describing the relationship between the linear and nonlinear component of the Doppler signal in a single synthetic aperture time. The higher linear proportion means that the echo is less spatial, and vice versa. According to the relationship between the slant range and Doppler, the maximum linear component in one synthetic aperture time is $|r_2|T_s/2$ according to Equation (2). The maximum nonlinear part is $|r_3|T_s^2/8$. Due to r₂ potentially being zero in some AOL positions, we mapped the ratio between these two into a negative exponential function, and introduced the spatial variance index:

$$g_{\mathrm{var}}=\exp \left\{-\frac{|r_3|}{|r_2|}\cdot \frac{T_s}{4}\right\}.$$

(27)

This means that when g_var is closer to 1, the linear component accounts for a larger proportion, and then the spatial variance weakens, and vice versa.

Another key factor dominating the performance of an imaging algorithm is the squint angle. As mentioned before, the squint angle is much smaller in GEO SAR, limited to ±8°, than in LEO SAR. This means that the traditional squint angle, or together with the elevation angle, cannot describe the squint characteristic in GEO SAR. Therefore, the ground squint angle was introduced to evaluate the squint condition in GEO SAR [18], which is the projection angle on ground of the traditional squint angle. The relationship of these angles is demonstrated in Figure 12.

The simulation parameters are listed in

Table 1. Simulation parameters.

Parameter	Value
Semi-major Axis (km)	42,164.3
Inclination (°)	60
Eccentricity	0
Longitude of Ascend Node (°)	100
Wavelength (m)	0.24
Longitude of Wenchuan County, China (°)	103.4 E
Latitude of Wenchuan County, China (°)	31.0 N

. The orbital parameters refer to the Global Earthquake Satellite System (GESS) projection [3]. The swath center is selected as Wenchuan County in China, which is a highly seismic zone requiring frequent observation.

5.2. Observation Duration Comparison

Based on the simulation parameters listed in Table 1, the full orbit observations by broadside and squint modes were simulated, and the results are shown in Figure 13. Here, the blue lines represent the nadir trajectory, the red lines represent the feasible observation durations corresponding to the respective observation modes, and the green point represents the target location.

For the broadside mode, there are only three positions that can be used to observe Wenchuan, as shown in Figure 13a. The AOLs of three positions are nearly 30°, 87°, and 148°. At these positions, the squint angle and elevation angle, denoted by $\left({\theta }_L,{\theta }_s\right)$, should be (1.92°,2.10°), (4.82°,0.51°), and (0.99°,−1.56°). The ground squint angles of these three positions are equal to zero. Applying the broadside mode yields an observation duration of about 0.59 h per orbit, which is 2.47% of the orbital period. In contrast, through adjusting the observing angles in different AOLs, the proposed algorithm can image Wenchuan for several strip areas, as shown in Figure 13b. Applying the proposed method, the observation duration increased to about 11.55 h per orbit, which is 48.15% of the orbital period. The results indicate that the proposed method significantly improves the observation ability of GEO SAR.

5.3. Imaging Results and Analysis

Figure 13 shows the feasible observation period in different modes. There is no need to clarify the feasible period in Figure 12a because so many broadside imaging algorithms have that ability. However, the feasible periods in Figure 12b were created by the proposed algorithm, which should be verified. In Figure 13b, three sections, FA, BC, and DE, are stated as the feasible periods. According to the simulation parameters listed in Table 1, we drew the spatial variance index g_var and ground squint angle ${\theta }_g$ curves in the full orbit, which are shown in Figure 14.

According to Section 5.1, the most difficult position for imaging are among A~F. Therefore, the main characteristics of these points are listed in

Table 2. Characteristics of candidate positions.

Index	AOL (°)	Longitute-Latitude of Nadir (°)	gvar	${\mathit{\theta }}_{\mathit{g}}$ (°)	${\mathit{\theta }}_{\mathit{L}}$ (°)	${\mathit{\theta }}_{\mathit{s}}$ (°)	${\mathit{T}}_{\mathit{s}}$ (s)
A	47.10	(81.18E, 39.38N)	0.955	26.95	3.18	1.30	365.08
B	67.29	(82.78E, 53.02N)	0.965	24.48	4.43	0.67	465.53
C	107.66	(114.83E, 55.61N)	0.965	−31.66	4.35	0.37	529.99
D	131.22	(119.07E, 40.65N)	0.959	−39.51	2.74	−0.43	428.13
E	185.05	(97.48E, 4.37S)	0.981	62.53	−3.61	−4.48	714.52
F	343.18	(108.23E, 14.51S)	0.964	−58.63	−4.15	5.52	722.26

.

According to Table 2, position A has the most severe spatial variance (g_var = 0.955), and position D has the largest ground squint angle (${\theta }_g$ = 62.53°). Hence, if the imaging qualification in A and D can be guaranteed, then all the other positions in the feasible period in Figure 13b can be focused well.

We employed the simulated scene illustrated in Figure 15 covering an area of 30 × 30 km to verify the imaging performance of the proposed method. The scene includes three point targets, T₁, T₂, and T₃, arranged at the lower left corner, the center, and the upper right corner, respectively. Two group experiments were carried out. The first was at position A. The algorithm designed by Zhanget al. [18] and the proposed algorithm were compared. The other experiment was at position D, and only the results by the proposed algorithm are shown. All other parameters were the same as those listed in Table 1 and Table 2.

At position A, the imaging results for the simulated scene given in Figure 15 obtained by the 2D nonlinear CS algorithm [18] and the proposed algorithm are presented in Figure 16 and Figure 17, respectively. Here, the first, second, and third rows correspond to points T₁, T₂, and T₃ illustrated in Figure 15, respectively. The first, second, and third columns represent contour maps, azimuth profiles, and range profiles, respectively. In both figures, contour maps denote the 2D focusing quality, and the profiles represent the focusing quality along the range and azimuth directions. The data given in Figure 16 and Figure 17 were evaluated according to resolution (res), peak to side lobe ratio (PSLR), and integrated side lobe ratio (ISLR) along the azimuth and range directions, and the results are presented in

Table 3. Evaluation of imaging quality (in position A).

		Azimuth Direction			Range Direction
Method	Target	Res (m)	PSLR (dB)	ISLR (dB)	Res (m)	PSLR (dB)	ISLR (dB)
Algorithm proposed by Zhang et al. [18]	T1	55.622	−15.129	−19.633	1.578	−12.963	−10.624
	T2	4.810	−12.525	−10.279	1.522	−13.279	−10.050
	T3	33.486	−0.809	−1.828	1.522	−9.213	−7.721
Proposed Algorithm	T1	4.775	−13.256	−10.295	1.534	−13.351	−10.138
	T2	4.800	−13.168	−10.153	1.534	−13.264	−10.039
	T3	4.973	−13.323	−10.091	1.534	−13.310	−10.119

. Here, PSLR is defined as the ratio of the peak level of sidelobes to the peak level of the main lobe, and the ISLR is the ratio of the total power in all sidelobes to the power in the main lobe.

The algorithm proposed by Zhang et al. [18] neglects the residual walk after removing the bulk range walk, and this negatively affects the azimuth spectrum of the imaging results, as indicated by Figure 16 and Figure 17 and Table 3. This issue becomes increasingly serious as the spatial variance increases. As shown in Figure 16, the swath center point T₂ is well focused. However, for points T₁ and T₃, which are at the edges of the swath, the signal in azimuth is almost defocused. Since the spatial variance is so severe that the traditional nonlinear chirp scaling does not work well, the imaging quality in range direction is degraded. However, the proposed algorithm preserves the focusing performance of the 2D nonlinear CS algorithm over the full scene. From the quality indices presented in Table 3, the imaging coherency is good over the entire swath. Since the fourth-order nonlinear chirp scaling is used in range direction, the focusing ability in range direction is improved compared to the algorithm proposed by Zhang et al. [18].

Position D has the largest ground squint angle, which should be simulated to verify the performance of proposed imaging algorithm. Since the performance comparison with the start-of-art algorithm has been presented in position A above, only the results of the proposed algorithm are shown here, in Figure 18 and

Table 4. Evaluation of imaging quality (in position D).

		Azimuth Direction			Range Direction
Method	Target	Res (m)	PSLR (dB)	ISLR (dB)	Res (m)	PSLR (dB)	ISLR (dB)
Proposed Algorithm	T1	4.656	−13.233	−9.695	1.318	−13.298	−10.097
	T2	4.866	−13.310	−10.253	1.309	−13.233	−10.031
	T3	4.952	−13.158	−10.250	1.318	−13.324	−10.139

. The simulation parameters are listed in Table 1 and Table 2, and the simulation layout is the same as shown in Figure 15.

The simulation results shown in Figure 18 and Table 4 together verify the imaging performance in the most severe squint looking conditions. It is obvious that our method has good imaging ability over the full scene. However, the focusing effect in the azimuth direction is little degraded compared to position A, as shown in the azimuth profiles in Figure 18. However, this performance degradation is acceptable for most applications.

6. Discussion

6.1. Accuracy

In order to correct the squint effect, the linear residual terms $k_{1,1,r}\cdot \Delta R$ and ${k_{1,1,a}|}_{\Delta R}\cdot \Delta t_c$ in Equation (3) were removed since these two terms depend on the locations of targets. Segmentation was adopted by applying Equations (16) and (25). After segmentation, each sub-block corresponds to an area of $\Delta R_{seg}\times \left(1-\chi \right)T_s$, as shown in Figure 19. When processing each sub-block of data, the follow-up filters were designed according to the center point P_C. Therefore, for all points except P_C, residual phase errors occur, which have to be controlled to guarantee the imaging accuracy.

Suppose a target located at $\left(\Delta t,\Delta R\right)$ has an accurate 2D frequency domain phase ${\mathsf{\Phi }}^a\left(\Delta t,\Delta R\right)$ and the polynomial expansion phase ${\mathsf{\Phi }}^e\left(\Delta t,\Delta R\right)$. In some positions where the spatial variance is minimal, the range walk removed echo data are like broadside data, which means $r_1\approx 0$ and ${\mathsf{\Phi }}^a\approx {\mathsf{\Phi }}^e$. However, for some severe spatial variant areas, $k_{1,1,r}\cdot \Delta R$ and ${k_{1,1,a}|}_{\Delta R}\cdot \Delta t_c$ change drastically. The expansion ${\mathsf{\Phi }}^e\left(\Delta t,\Delta R\right)$ is no longer accurate. Segmentation can be used to correct it. Therefore, the correct approach to weigh the segmentation is to improve the quality of the phase error of the corner points. That is,

$$\Delta \mathsf{\Phi }\left(\chi ,\Delta R_{seg}/\Delta R_B\right)=\underset{P\in P_M}{\max }\left|{\mathsf{\Phi }}_P^e\left(\Delta t\pm \frac{T_s}{2}\left(1-\chi \right),\Delta R\pm \frac{1}{2}\Delta R_{seg}\right)-{\mathsf{\Phi }}_P^a\left(\Delta t\pm \frac{T_s}{2}\left(1-\chi \right),\Delta R\pm \frac{1}{2}\Delta R_{seg}\right)\right|,$$

(28)

where

$$\{\begin{array}{l}{\mathsf{\Phi }}^a\left(\Delta t,\Delta R\right)=-\frac{4\pi (f_{\tau }+f_0)}{c}\left[r_0-{\displaystyle \sum _{n=1}^9\frac{A_n}{n+1}{\left(-\frac{c{f}_t}{2(f_0+f_{\tau })}-r_1\right)}^{n+1}}\right]\\ {\mathsf{\Phi }}^e\left(\Delta t,\Delta R\right)=\mathsf{\Phi }\left(f_t,f_{\tau }\right)\end{array},$$

(29)

and the $\Delta t$ and $\Delta R$ are inherent in r_n and A_n, which are defined in Equations (3) and (A9). $\mathsf{\Phi }\left(f_t,f_{\tau }\right)$ is shown in Equation (A7).

The relative phase $\Delta \mathsf{\Phi }$ should be controlled at a very low level, which means segmentation contributes to an accurate expression without introducing phase error. To this end, Equation (28) should obey

$$\Delta \mathsf{\Phi }\left(\chi ,\Delta R_{seg}/\Delta R_B\right)<\frac{\pi }{4M},$$

(30)

where M > 1 is an adjustable coefficient.

Based on this criterion, parameters $\Delta R_{seg}$ and $\chi$ can be selected to ensure the residual parts after segmenting are small enough for subsequent operations. It also shows that segmentation does not introduce error by this means. Using the simulation examples demonstrated in Section 5, the simulated relationship among $\Delta R_{seg}$, $\chi$, and $\Delta \mathsf{\Phi }$ is demonstrated in Figure 20.

The final imaging quality does not depend on Equation (30); instead, quality depends on the fourth-order NCS in range direction and time-frequency scaling in the azimuth direction. Equation (30) only ensures that segmentation has no effect on the subsequent image quality. Since the fourth-order NCS and time-frequency are designed for imaging a broad scene, the imaging quality of targets in small 30 × 30 km areas can be ensured in most circumstances. This means that the marginal points of sub-blocks can be focused. The error between sub-blocks can only be distorted, which can be resolved by geometric calibration before stitching.

6.2. Computational Load

The computational load is one of the main factors affecting the practical application of the algorithm. To multiply a pair of complex numbers requires six floating point operations per second (FLOPs), whereas only one FLOP is needed for complex addition. Specifically, creating an N points FFT (or IFFT) operation requires 5Nlog₂(N) FLOPs [20]. Suppose the echo data has N_azi and N_rng sampling points in the azimuth and slant range directions, respectively, after range walk removal. For the BP algorithm, with eight-times interpolation in the range direction, and imaging an area of N_rng × N_rng size, the computational load of BP is:

$$C_{BP}=45N_{azi}{N}_{rng}{\log }_2\left(N_{rng}\right)+7N_{azi}{N}_{rng}^2.$$

(31)

However, for the algorithm proposed, there are 4N_azi times N_rng points FFT operations, 10N_rng times N_azi points FFT operations, and 14N_aziN_rng times multiplications. Therefore, the computational load is:

$$C_p\approx 20N_{azi}{N}_{rng}{\log }_2\left(N_{rng}\right)+50N_{azi}{N}_{rng}{\log }_2\left(N_{azi}\right)+84N_{azi}{N}_{rng}.$$

(32)

For specific simulation circumstances, the computational load can be numerally determined. Take the simulation circumstance shown in Figure 17 as an example. The azimuth sampling number is 6000 and the slant range sampling number after range removal is 2000. Therefore, the computational load of BP and the proposed algorithm are 173.9 and 11.17 GFLOPs, respectively. The computational burden of the proposed algorithm is only one-fifth that of the BP algorithm, which is suitable for practical use. In addition, many segmenting operations exist that can be easily used by parallel computing. This further shortens the time required for the method when implementing in modern computers.

7. Conclusions

We proposed a novel squint mode GEO SAR imaging method, where the bulk range walk is first corrected on receive by adjusting the starting time of the receive window and then applying first-order filtering, RCMC, range compression, partial dechirp, and azimuth compression to obtain the focused single looking complex (SLC) image. Simulation results demonstrated that the proposed algorithm achieves a resolution less than 5 m in the azimuth direction over a 30 × 30 km swath.

In addition to the performance, the feasibility of the proposed algorithm was investigated. Firstly, up-to-date technology can make the antenna have the ability to scan in a range of ±8° in the azimuth direction, which indicates that the staring mode can be realized. Secondly, the low PRF employed with GEO SAR provides ample scope for adjusting the starting time of the receive window and removing the range walk on receive. Lastly, the ephemeris data and geographic information of the target zone can be acquired to calculate the coefficients of the range model and each filter in the proposed algorithm, and focused imaging results can be produced.

The staring mode of GEO SAR enables the continuous observation of the region of interest. Compared with the broadside mode, the proposed method increases the observation duration for Wenchuan County in China by nearly 20 times. As a result, the proposed method enables GEO SAR to conduct long-term observation of target zones, which is essential for the management of global disasters and related affairs. For earthquake-prone areas, like Wenchuan, the starting mode can continuously provide images over a long time period after an earthquake. Some disaster information, including dammed lakes, landslides, and mudslides, can be updated continuously, which is central to disaster relief and rescue.

Another distinctive characteristic of GEO SAR is the short temporal resolution, which is vital to the application of interferometry, such as vegetation measurement and crust deformation detection. However, interferometry cannot be realized by applying several images obtained at different squint angles within one orbit period and coherent change detection needs to use the same observation geometry. For this reason, for interferometric applications, it is necessary to wait at least one orbit period (one day), similarly to the polar orbit repeat pass interferometry case.