Correcting Spatial Variance of RCM for GEO SAR Imaging Based on Time-Frequency Scaling

Compared with low-Earth orbit synthetic aperture radar (SAR), a geosynchronous (GEO) SAR can have a shorter revisit period and vaster coverage. However, relative motion between this SAR and targets is more complicated, which makes range cell migration (RCM) spatially variant along both range and azimuth. As a result, efficient and precise imaging becomes difficult. This paper analyzes and models spatial variance for GEO SAR in the time and frequency domains. A novel algorithm for GEO SAR imaging with a resolution of 2 m in both the ground cross-range and range directions is proposed, which is composed of five steps. The first is to eliminate linear azimuth variance through the first azimuth time scaling. The second is to achieve RCM correction and range compression. The third is to correct residual azimuth variance by the second azimuth time-frequency scaling. The fourth and final steps are to accomplish azimuth focusing and correct geometric distortion. The most important innovation of this algorithm is implementation of the time-frequency scaling to correct high-order azimuth variance. As demonstrated by simulation results, this algorithm can accomplish GEO SAR imaging with good and uniform imaging quality over the entire swath.


Introduction
Geosynchronous synthetic aperture radar (GEO SAR) operates at an altitude~36,000 km [1]. Compared with a low-Earth orbit (LEO) SAR, greater coverage can be achieved by GEO SAR because of its much higher orbit [2]. Furthermore, GEO SAR can guarantee observation of the same location every 24 h with the same incidence angle, which cannot be realized by LEO SAR [3]. Given its special characteristics, GEO SAR has attracted much attention [4]. It has become part of a global earthquake satellite system to monitor the global seismic state, as proposed by NASA and JPL in 2003 [5]. Another GEO SAR system called Geosynchronous Earth Monitoring by Interferometry and Imaging (GEMINI) was put forward in 2012 to acquire Earth surface data through GEO SAR interferometry [6]. For GEO SAR systems, imaging is always a major problem, the key to which is correction of the spatial variance of range cell migration (RCM).
RCM determines the distribution of echoes in the time and frequency domains. The spatial variance of RCM causes the spectrums corresponding to targets at different locations to be different. Therefore, in order to accomplish precise focusing in the frequency domain, the spatial variance must be corrected to eliminate the difference of RCM between any target in the swath and the reference point, usually the swath center. For LEO SAR with zero-Doppler steering, RCM is spatially variant along the range direction only, which can be processed by algorithms, such as range Doppler (RD) [7],

Echo Model and Spatial Variance Analysis
For SAR, after demodulation, the echo corresponding to an isolated point target can be represented by: where η and τ denote the slow time along the azimuth and the fast time along the range, respectively. The constant σ is the backscattering coefficient of the target, c 0 is light speed, and λ is wavelength. a pτq represents the transmitted signal. Here, the chirp signal is the transmitted signal, which implies a pτq " rect " τ{T p ‰ exp " jπK r τ 2 ‰ . K r is the linear frequency modulation rate, T p is pulse duration, and rect r¨s denotes the rectangular envelope. η c is the beam crossing time and T s is the synthetic aperture time.
In Equation (1), R actual pηq is the equivalent slant range between the spaceborne SAR and the target, which is the average of the one-way slant range when transmitting a pulse and that when receiving the echo. Usually R actual pηq can be approximated by a polynomial: where r n is the nth-order coefficient, and N r denotes the order. All SAR imaging algorithms are derived based on an appropriate slant range model, like Equation (2). The higher the order, the better the fit, which leads to better imaging performance. However, a high order increases the complexity of imaging algorithms. Therefore, determining an N r that balances imaging quality and complexity of the algorithm is a major challenge. Usually the aperture time of LEO SAR is so short that Equation (2) with N r " 2 is adequate to approximate the slant range, which indicates that the phase error induced by the approximation is too small to affect imaging quality [18]. For GEO SAR, orbital altitude increases the aperture time over one hundred times. For example, if the ground resolution is 2 m, the aperture times for LEO and GEO SARs are about 6 s and 750 s, respectively. Thus, to describe the much more complicated relative motion with the longer aperture time, the order of Equation (2) must be redetermined by evaluating the impact of different N r on imaging quality, which can be represented by the resolution, peak side-lobe level ratio (PSLR), and integral side-lobe level ratio (ISLR) [18].
The following matching filtering is adopted: where b represents convolution and N r may equal 4, 5 or 8. When N r approaches 8, Equation (3) indicates ideal filtering, which can achieve optimal imaging quality per the theory of matching filtering [18]. Evaluation results of applying Equation (3) and parameters in Table 1 are shown in Figure 1. Azimuthal resolution, PSLR, and ISLR achieved with N r " 5 are as nearly identical to results of the ideal filtering during the entire orbital period. However, when N r " 4, results are much worse. Therefore, the fifth-order approximation to R actual pηq is adequate to acquire optimum imaging quality.

Spatial Variance in the Time Domain
Given the invalidity of the "stop-and-go" approximation, n r ( n  0 ) in Equation (2)  (a-c) show results of azimuthal resolution, PSLR, and ISLR during the entire orbital period, respectively. Results achieved with N r " 5 are as nearly identical to those attained by ideal filtering. With N r " 4 , results are much worse.  (1) and (2), the echo model can be expressed as:

Spatial Variance in the Time Domain
Given the invalidity of the "stop-and-go" approximation, r n (n ě 0) in Equation (2) can be expressed as follows (see Appendix A): where: and x¨y denotes the dot product. C m n " n!{ rpn´mq!m!s and r n,sin denotes the nth-order one-way slant range coefficient.
represents nth-order derivatives of the position vectors of the SAR satellite at the beam crossing time η c and Ñ R pnq g_tar represents the position vector of the target. Usually the attitude steering is applied for GEO SAR [19], and makes the Doppler centroid zero, which leads to r 1 " 0. Figure 2 illustrates the GEO SAR observation geometry. The beam crossing time corresponding to a target is the moment when the zero Doppler plane crosses the target. At the beam crossing time corresponding to the swath center, the distance from GEO SAR to the swath center is defined as the reference range, i.e., r 0,re f , as shown in Figure 2a. When the zero Doppler plane crosses any target in the swath at η c , the distance from SAR to the target is denoted by r 0 , as shown in Figure 2b. As a result, the location of an isolated target can be uniquely represented by the beam crossing time η c and the corresponding distance r 0 .
Equation (5) demonstrates that r n depends on the target position, indicating that r n is spatially variant. Therefore, using polynomial fitting, r n can be expressed as a function of ∆R and η c , which are adopted as spatial variables along the range and azimuth directions, respectively. That is: r 0 " r 0,re f`∆ R r 1 " 0 r 2 " r 2,re f`k2,1,r¨∆ R`k 2,2,r¨p ∆Rq 2`k 2,1,a¨ηc`k2,2,a¨p η c q 2 r 3 " r 3,re f`k3,1,r¨∆ R`k 3,1,a¨ηc`k3,2,a¨p η c q 2 r 4 " r 4,re f`k4,1,r¨∆ R`k 4,1,a¨ηc r 5 " r 5,re f (6) where ∆R " r 0´r0,re f . r n,re f , k n,m,a and k n,m,r can be calculated by polynomial fitting based on ephemeris data and geographic information of the swath. All coefficients are spatially invariant except Sensors 2016, 16, 1091 5 of 30 k n,m,a , which varies with ∆R. Equation (6) demonstrates that there is nonlinear spatial variance in the slant range coefficients, which leads to the same condition in the phase spectrum.

Spatial Variance in the Frequency Domain
Equation (6) demonstrates that coefficient r n is two-dimensionally spatially variant, and so does RCM in the time domain. By implementing the Fourier transform (FT) on s original pη, τq in range, the signal in the range-frequency domain can be expressed as follows: where f τ is range frequency, B r is bandwidth of the transmitted signal, and f 0 is the carrier frequency. Then by performing FT along the azimuth, the two-dimensional spectrum is as follows (see Appendix B): where A n is defined as follows: "´r 3 2 r 6`1 5r 5 r 2 2 r 3`1 0r 2 2 r 2 4´1 05r 2 r 2 3 r 4`1 05r 4 3 ‰ A 6 "¨¨¨(

9)
Sensors 2016, 16, 1091 6 of 30 W a " f η ‰ is the azimuth spectrum amplitude. Since it is only concerned with azimuth focusing, its representation will be given in Section 4. By applying series expansion [20], Equation (8) can be organized in the form of a series of f τ , which is: where: ϕ 0 , φ 1 , φ 2 and φ k denote the azimuth modulation phase, RCM, range linear frequency-modulated phase and high-order range frequency-modulated phase, respectively. Equations (9), (11) and (12) demonstrate that ϕ 0 and φ k depend on r n . Therefore, ϕ 0 and φ k are spatially variant, which can also be explicitly expressed as a function of ∆R and η c : All coefficients in Equation (13) depend on f η , and some of them also vary with ∆R. Details are given in Appendix B. The phase difference between Equations (11) and (13) is <0.012π, indicating that these approximations in Equation (13) will not affect imaging quality.

Basic Methodology of Correcting Spatial Variance
According to Equation (2), RCM is determined by slant range coefficients (i.e., r n ). In order to correct the spatial variance of RCM, this paper adopts time-frequency scaling to modify r n .
To introduce this idea, a one-dimensional signal s ptq is assumed to be: where T is the signal duration time. t c is the center time and d n is the nth-order time-domain phase coefficient. According to Equation (6), d n can be assumed quadratically variant with t c . The aim of the time-frequency scaling method is to remove the spatial variance in d n , which indicates that d n doesn't depend on t c after scaling. By applying series reversion method [21], the corresponding frequency-domain spectrum of Equation (14) can be attained as: D n is the frequency-domain spectrum phase coefficient and can be obtained by applying the method in Appendix B. In order to avoid time-domain aliasing, D 1 will not be modified by scaling. Modification of frequency-domain spectrum phase coefficients leads to modification of time-domain phase coefficients. As a result, the following frequency scaling function can be used: After multiplication by Equation (16), Equation (15) becomes: By applying series reversion method to Equation (17), the following expression can be obtained: where: d 1 n is the reconstructed nth-order time-domain phase coefficient. Like d n , d 1 n is also spatially variant with t c . According to Equation (6), spatial variance is up to the second order along both the range and the azimuth. As a result, d 1 n can be assumed a second-order function of t c , i.e., where d 1 n,re f corresponds to the swath center. Then, the spatial difference between d 1 n and d 1 n,re f is: In order to eliminate the first-order spatial variance of d 1 n , the time scaling function to be multiplied with Equation (18) is designed as: By multiplying Equations (18) and (22), the linear component of ∆d 1 n will be removed. Therefore, for correcting the linear spatial variance of all d 1 n (n ě 2), the complete time scaling function is: After multiplying Equation (18) with Equation (23), the signal becomes: where: The new time-domain phase coefficient is: Similar to d 1 n , d 2 n also varies with t c . By applying the series reversion method to Equation (24), the spectrum is: where B m can be acquired by the method in Appendix B.
Because the time scaling has removed the linear spatial variance, B m satisfies: In order to remove the quadratic spatial variance in B m (m " 2, 3), the following equation should be satisfied by assigning the appropriate value of Z n : According to Equation (13), φ k (k ą 3) is spatially invariant in GEO SAR imaging. And after time-frequency scaling, the linear and quadratic spatial variance has been removed from B m (m " 2, 3), as shown in Equations (28) and (29). Therefore, focusing for the whole swath can be accomplished in the frequency domain. Although the time-frequency scaling is only applied to correct the azimuth variance in this section, it can also be applied for the range variance correction, as shown in Section 4.

Spatial Variance Correction and GEO SAR Imaging
Based on the basic idea in Section 3, a GEO SAR imaging algorithm composed of five steps is proposed. The first is to eliminate linear azimuth variance of φ k through the first azimuth time scaling. The second is to achieve RCM correction and range compression. The third is to correct residual azimuth variance by the second azimuth time-frequency scaling. The fourth and final steps are to accomplish azimuth focusing and correct geometric distortion.

First Azimuth Time Scaling
In Equation (10), φ k determines the coupling between range and azimuth. The first step applies the time scaling to remove linear azimuth variance of φ k and guarantees the quality of RCM correction and range compression. In the range-frequency and azimuth-time domain, the echo can be expressed as Equation (7). The first azimuth scaling function is designed as: where: After multiplying Equations (7) and (30), the signal becomes: r 1 n represents the modified nth-order slant range coefficient, which is: By performing the azimuth FT on Equation (32), the signal in the two-dimensional frequency domain is: where ϕ 1 0 and φ 1 k (1 ď k ď 9) can be attained by replacing r n with r 1 n in Equation (11) and polynomial fitting, i.e., In Figure 3a, three curves show RCMs corresponding to the three targets T A , T B and T C . And, they are in the same range cell. These curves do not coincide because of azimuthal variance in the RD domain, as demonstrated in Figure 3b. After the first azimuth time scaling, RCMs are corrected to be the same and shifted by various offsets. As a result, they are approximately parallel in the RD domain and cross different range cells, as shown in Figure 3c. The range offset causes imaging distortion and is corrected in the final step. Equation (35), and the residual quadratic azimuth variance makes RCM variation be less than one quarter of a range cell in azimuth. Therefore, RCMs corresponding to targets in the same range cell can be regarded as identical [22], indicating that the first azimuth time scaling is beneficial to RCMC and range compression in the next step.
However, for targets in the same range cell a range offset in the RD domain is induced by the first azimuth scaling, which is: In Figure 3a, three curves show RCMs corresponding to the three targets TA, TB and TC. And, they are in the same range cell. These curves do not coincide because of azimuthal variance in the RD domain, as demonstrated in Figure 3b. After the first azimuth time scaling, RCMs are corrected to be the same and shifted by various offsets. As a result, they are approximately parallel in the RD domain and cross different range cells, as shown in Figure 3c. The range offset causes imaging distortion and is corrected in the final step.

RCM Correction and Range Compression
In order to correct the range variance of 1   , 2   , and 3   , the following compensation function can be applied based on the thought of time-frequency scaling which is demonstrated in Section 3: where:

RCM Correction and Range Compression
In order to correct the range variance of φ 1 1 , φ 1 2 , and φ 1 3 , the following compensation function can be applied based on the thought of time-frequency scaling which is demonstrated in Section 3: where: and f ηre f is the Doppler centroid corresponding to the swath center. By applying range IFT on the multiplication of Equations (34) and (37), the signal in the RD domain can be formulated as follows (see Appendix C): where: Coupling of τ 1 and τ 1 denotes range variance. Per Equation (24), the function for correcting the range variance is: where τ 1`τ 1 " τ`φ 1 1,re f , and: The multiplication of Equations (38) and (39) can be approximated as: The last term of Equation (41) will result in asymmetric range sidelobe. To further eliminate the impact of range variance in the last term of Equation (41), data segmentation can be adopted in the RD domain. Every segment corresponds to a sub-swath whose width in the slant range direction is <30 km. The phase variation induced by the term with τ 1 pτ 1 q 3 is <0.04π in the sub-swath, which indicates that the last term in Equation (41) can be ignored in the processing of every segment.
By applying the range FT and the compensation function to Equation (41) successively, the range variance in every segment can be corrected. The compensation function is: where ∆R sub is distance from the swath center to sub-swath center in the slant range direction. Then, compensation results of data segments are stitched together. A uniform filter can be applied to accomplish RCM correction and range compression for the entire swath, which is: By applying Equation (43), the spectrum in the RD domain becomes: The last three terms in Equation (44) can be further compensated for by the following function: Compensated for by Equation (45), the signal is:

Second Azimuth Time-Frequency Scaling
After RCM correction and range compression, the third step is to further correct the residual azimuth variance in ϕ 1 0 . To compensate for the additional phase induced by the first azimuth time scaling, the following function is adopted, which is: where the last phase is spatially variant and leads to azimuth defocusing. Based on the basic methodology in Section 3, to correct the azimuthal variance in Equation (48), the frequency scaling is adopted, which is: where: and: After multiplying Equation (48) with Equation (49), the signal in the time domain is: where: Then, the azimuth time scaling function is applied to Equation (52), which is: where: where: By applying azimuth FT to Equation (57), the signal becomes: where P 3 m can be obtained by replacing r n with r 3 n in Equation (12). The spectrum envelope can be expressed as: In Equation (59), the phase variation induced by azimuthal variance is <0.024π in the entire swath, and does not affect azimuth focusing. As a result, the azimuth focusing can be implemented in the frequency domain.

Azimuth Compression
Equation (59) is spatially invariant and can be focused by applying a uniform azimuth matching filter for the entire swath, i.e., where: P r c ,m can be attained by replacing r n with r c,n in Equation (12). The result of matching filtering is: By implementing azimuth IFT, the focusing result can be expressed as:

Geometric Correction
Equation (65) indicates that the echo has been completely focused. However, the target location is shifted in both azimuth and range directions. As illustrated in Figure 4, two targets, T A and T B , are originally in the same range cell, and η c for T B is zero. In the focusing result, T B is still at its original position. However, T A has been shifted to another position T 1 A . The offsets along the range and azimuth directions are t con and´P 3 1 , respectively, revealing geometric distortion in the focused imaging that should be corrected.
In Figure 4, the dashed line represents the offset trajectory, which can be described by: where:

67)
Q 1,2,a and Q 1,3,a are the second-order and the third-order spatial variance coefficients of P 3 1 along the azimuth, respectively. They can be obtained by replacing r n with r 3 n in Equations (9) and (12). Equation (65) can be transformed into the range-frequency and azimuth-time domain, multiplied by the geometric correction function, i.e., and transformed back into the time domain to achieve the final imaging result, which is: A flowchart of the proposed algorithm is shown in Figure 5. In Figure 4, the dashed line represents the offset trajectory, which can be described by: and 1,3,a Q are the second-order and the third-order spatial variance coefficients of 1 P along the azimuth, respectively. They can be obtained by replacing n r with n r  in Equations (9) and (12). Equation (65) can be transformed into the range-frequency and azimuth-time domain, multiplied by the geometric correction function, i.e., and transformed back into the time domain to achieve the final imaging result, which is: Figure 4. Illustration of geometric distortion. Two targets (T A and T B ) are orignally located in the same range cell, and η c for T B is zero. After imaging by proposed algorithm, T B is still at its original location. However, T A is focused at T 1 A . Offsets in the range and azimuth directions are t con and´P 3 1 , respectively. All targets originally on red line are on dashed line after focusing. Therefore, the focusing results must be corrected from the dashed to red line by geometric correction, from T A flowchart of the proposed algorithm is shown in Figure 5.

Ephemeris data
Geographic information Figure 5. Flowchart of the proposed algorithm.

Simulation Parameters
Simulation parameters are listed in Table 2, which refer to the global earthquake satellite system [5]. The center of the swath is set to be 108.5° E and 35.3° N. Besides these parameters, the spatial variance of RCM also depends on the swath size. As the swath becomes wider, higher-order spatial

Simulation Parameters
Simulation parameters are listed in Table 2, which refer to the global earthquake satellite system [5]. The center of the swath is set to be 108.5˝E and 35.3˝N. Besides these parameters, the spatial variance of RCM also depends on the swath size. As the swath becomes wider, higher-order spatial variance may occur along both range and azimuth. As discussed in Section 4, the proposed algorithm can correct linear and quadric spatial variance. Therefore, in order to avoid the cubic and higher-order spatial variance, a simulated swath covering an area of 83 km (azimuth)ˆ86 km (range) is adopted to verify the algorithm with parameters in Table 2. It is certain that the swath size for the proposed algorithm varies with observation parameters. For example, if the center time is 0 and the swath center is at the equator, the swath size can reach 150 km (azimuth)ˆ150 km (range) (see Appendix D). The simulated swath is portrayed in Figure 6. Nine point targets, from T 1 to T 9 , are deployed. T 5 is at the swath center and other points are at the swath edge.  The simulated swath is portrayed in Figure 6. Nine point targets, from 1 T to 9 T , are deployed.

5
T is at the swath center and other points are at the swath edge.  Figure 6. T5 is at swath center. T1 is 41.5 km and 43 km away from T5 along azimuth and range, respectively.

Imaging Results
Range and azimuth profiles corresponding to every target are illustrated in Figure 7, and evaluation results are listed in Table 3. The broadening coefficient is the ratio of the achieved resolution to the ideal resolution. Suppose that a  and ,  Table 3, the broadening coefficients are almost 1, which signifies no resolution loss in imaging. The ideal PSLR should be −13.26 dB. The loss of PSLR equals the absolute difference

Imaging Results
Range and azimuth profiles corresponding to every target are illustrated in Figure 7, and evaluation results are listed in Table 3. The broadening coefficient is the ratio of the achieved resolution to the ideal resolution. Suppose that ρ a and ρ a,ideal are the azimuth resolutions achieved by the proposed algorithm and the ideal imaging, respectively. The azimuth broadening coefficient is defined as ρ a {ρ a,ideal . The range broadening coefficient can be achieved by the same way. As shown in Table 3, the broadening coefficients are almost 1, which signifies no resolution loss in imaging. The ideal PSLR should be´13.26 dB. The loss of PSLR equals the absolute difference between the actual and ideal PSLRs. Table 3 shows that the maximal loss of range and azimuth PSLRs is ď0.2 dB, and ď0.1 dB, indicating good focusing quality. The difference of range PSLRs, azimuth PSLRs, range ISLRs, and azimuth ISLRs over the whole swath is ď0.22 dB, ď0.17 dB, ď0.29 dB, and ď0.28 dB respectively, indicating uniform imaging quality.  In order to further demonstrate the advantage of the proposed algorithm, it is useful to compare the imaging results obtained by different techniques. Here the imaging results of T3 and T5 are compared by applying the proposed algorithm, and other two algorithms developed by Hu [15] and Li [17].
Theoretically, any algorithm can achieve good focusing quality for T5, because it is at the swath center. For the same target at the swath edge the imaging performances of different algorithms may be different. Imaging profiles corresponding to T5 and T3 are shown in Figures 9 and 10. Evaluation results are listed in Tables 4 and 5. It is shown that for T5 three algorithms have basically the same imaging performance. However, for T3, azimuth defocusing occurs by applying the algorithm in [15]. The algorithm in [17] induces defocusing along azimuth and range directions, because the range variance is not corrected adequately and furthermore the azimuth focusing quality is affected. By comparison, the proposed algorithm has the best performance between these three algorithms.  Usually for LEO SAR azimuth ISLRs and range ISLRs are almost equal. The situation is different for GEO SAR with a resolution of 2 m. Because of the wide-angle observation, the two-dimensional amplitude spectrum is not rectangular, and bifurcation exists along azimuth, as shown in Figure 8. The bifurcation causes the energy of sidelobes to disperse in two directions, while the mainlobe is not affected. As a result, azimuth ISLRs are better than range ISLRs in Table 3.   In order to further demonstrate the advantage of the proposed algorithm, it is useful to compare the imaging results obtained by different techniques. Here the imaging results of T 3 and T 5 are compared by applying the proposed algorithm, and other two algorithms developed by Hu [15] and Li [17].
Theoretically, any algorithm can achieve good focusing quality for T 5 , because it is at the swath center. For the same target at the swath edge the imaging performances of different algorithms may be different. Imaging profiles corresponding to T 5 and T 3 are shown in Figures 9 and 10. Evaluation results are listed in Tables 4 and 5. It is shown that for T 5 three algorithms have basically the same imaging performance. However, for T 3 , azimuth defocusing occurs by applying the algorithm in [15]. The algorithm in [17] induces defocusing along azimuth and range directions, because the range variance is not corrected adequately and furthermore the azimuth focusing quality is affected. By comparison, the proposed algorithm has the best performance between these three algorithms.   (c) Figure 9. Imaging profiles corresponding to T5. The first column represents the two-dimensional point spread function. The second and third columns represent azimuth and range profiles. (a-c) are achieved by applying the proposed algorithm, the algorithm in [15] and the algorithm in [17] respectively. (a) (b) Figure 9. Imaging profiles corresponding to T 5 . The first column represents the two-dimensional point spread function. The second and third columns represent azimuth and range profiles. (a-c) are achieved by applying the proposed algorithm, the algorithm in [15] and the algorithm in [17] respectively.
(c) Figure 9. Imaging profiles corresponding to T5. The first column represents the two-dimensional point spread function. The second and third columns represent azimuth and range profiles. (a-c) are achieved by applying the proposed algorithm, the algorithm in [15] and the algorithm in [17] respectively. (c) Figure 10. Imaging profiles corresponding to T3. The first column represents the two-dimensional point spread function. The second and third columns represent azimuth and range profiles. (a-c) are achieved by applying the proposed algorithm, the algorithm in [15] and the algorithm in [17] respectively.

Computational Load
Computational load is a key element to restrict the application of an algorithm. Although the chirp scaling algorithm (CSA) [12] cannot achieve the 2 m resolution for GEO SAR, it is worth comparing CSA and the proposed algorithm from the aspect of the computational load, because CSA is recognized as an efficient algorithm and has been widely applied. The back projection algorithm (BPA) [23] is also compared here, because BPA can achieve same focusing quality in the time domain.
Computational load is evaluated according to the complex multiplication and addition in the Figure 10. Imaging profiles corresponding to T 3 . The first column represents the two-dimensional point spread function. The second and third columns represent azimuth and range profiles. (a-c) are achieved by applying the proposed algorithm, the algorithm in [15] and the algorithm in [17] respectively.

Computational Load
Computational load is a key element to restrict the application of an algorithm. Although the chirp scaling algorithm (CSA) [12] cannot achieve the 2 m resolution for GEO SAR, it is worth comparing CSA and the proposed algorithm from the aspect of the computational load, because CSA is recognized as an efficient algorithm and has been widely applied. The back projection algorithm (BPA) [23] is also compared here, because BPA can achieve same focusing quality in the time domain.
Computational load is evaluated according to the complex multiplication and addition in the algorithm. Multiplication of two complex numbers and addition of two real numbers need 6 FLOPs and 1 FLOP, respectively. A FFT or IFFT with a length of N points needs 5Nlog 2 pNq FLOPs [18].
Suppose sampling numbers along azimuth and range are N azi and N rng , respectively. The computational loads of the proposed algorithm, CSA, and BPA with 8-fold interpolation, are respectively: In order to achieve 2 m resolution and a swath of 80 kmˆ80 km, N azi and N rng should be 140,000 and 50,000 at least. According to Equation (70), analysis results are listed in Table 6. Although the computational load of the proposed algorithm is more than twice that of CSA, the increasement is acceptable because of the significant improvement of imaging quality. And the computational load is about 1/1000 of that of BPA, indicating that the proposed algorithm is efficient.

Conclusions
This work models the spatial variance in the time and frequency domains based on a fifth-order polynomial slant range model. And a GEO SAR imaging algorithm is proposed, whose basic method is to correct the linear and quadratic spatial variance of RCM in the range and azimuth directions based on time-frequency scaling. As demonstrated by simulation results, this algorithm can accomplish GEO SAR imaging with good and uniform imaging quality over the entire swath.
The algorithm can be applied in the squint mode for more flexible observation, although it is developed under the condition of zero Doppler centroid. For the squint mode, the Doppler centroid is spatially variant. As a result, data segmentation has to be used to divide the echo into blocks. Every block is processed by linear RCM correction [24] and the proposed algorithm successively. Imaging results of blocks are mosaicked to form the final image.

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

Appendix A
Let R t pηq denote the one-way slant range between the satellite and target when transmitting a pulse, whose representation is: where Ñ R g_sat pηq and Ñ R g_tar are position vectors of the satellite and the target in the Earth Centered Rotating (ECR) coordinate system [18], respectively. By applying series expansion, Equation (A1) can be expressed as: R t pηq " r 0,sin`8 ÿ n"1 r n,sin n! pη´η c q n (A2) r n,sin denotes the nth-order coefficient. Furthermore, R t pηq can also be expressed as: Let R r pη`∆ηq denote the one-way slant range when receiving the echo, where ∆η denotes the round-trip delay time. For GEO SAR, the farthest slant range is about 41,680 km, which makes the round-trip delay time less than 0.28 s. As a result, the radius of the GEO SAR trajectory during the round-trip delay time is approximately invariable, which indicatesˇˇˇˇÑ R g_sat pη`∆ηqˇˇˇˇ«ˇˇˇˇÑ R g_sat pηqˇˇˇˇ. By combining Equations (A1) and (A2), R r pη`∆ηq equals: R t pηq and R r pη`∆ηq also satisfies the following propagation equation: Combining Equations (A4) and (A5), the following equation can be achieved: The satellite position at η`∆η can be approximated as: where Ñ V g_sat pηq and Ñ A g_sat pηq are the velocity and acceleration vectors of the satellite, respectively.
By replacing the component of Ñ R g_sat pη`∆ηq´Ñ R g_sat pηq in Equation (A6) with Equation (A7), ∆η can be solved: , ∆η can be approximated as: By substituting Equations (A2) and (A11) into Equation (A9), the equivalent slant range R pηq equals: R pηq " c 0 ∆η{2 " r 0,sin`8 ř where k n denotes From Equation (B4), a set of coefficient equations is acquired as follows: By solving Equation (B5), A i satisfies: Substituting Equation (B3) into Equation (B1), the spectrum phase in the two dimensional frequency domain is acquired: By applying series expansion, Equation (B7) can be expressed as a polynomial of f τ : and: A n p´r 1 q n`1 n`1 P 1 " Because P m depends on r n and the highest orders along range and azimuth are both 2 shown in Equation (6), P m can be modeled as: P m " P m,re f`Qm,1,r ∆R`Q m,2,r p∆Rq 2`Q m,1,a | ∆R η c`Qm,2,a | ∆R pη c q 2 , 0 ď m ď 10 (B11) According to Equations (B6) and (B10), every coefficient in Equation (B11) can be obtained. Then coefficients in Equation (13) can be acquired: and: φ k,re f " p´1q k   Figure C1. Illustration of important curves and variables in RD domain. τ 1 mig,re f and τ 1 mig represent migration curves corresponding to swath center and any target in the swath, respectively. τ 1 mig,c and τ 1 mig,s represent migration curves corresponding to any target whose η c " 0 before and after RCMC, respectively. τ 1 is the difference between τ 1 mig,s and τ 1 mig,re f . 2r pt pη c q{c 0 is the range offset induced by first azimuth scaling.
Considering the spatial variance model in Equation (35), we can represent τ 1 as: Combining Equations (35) and (C4), the differential RCM of the range cell center defined as where τ 1 " τ´τ 1 mig,s . K m and J m can also be represented as functions of τ 1 : # K m " K mre f`Ks τ 1 J m " J mre f`Js τ 1 (C8) Based on Equations (C7) and (C8), the signal in Equation (C2) can be represented as Equation (38).

Appendix D
This appendix will show the swath size achieved by the proposed algorithm can reach 150 km (azimuth)ˆ150 km (range). Simulation parameters are listed in Table D1. The simulated swath is portrayed in Figure D1, and the swath center is at the equator.

Appendix D
This appendix will show the swath size achieved by the proposed algorithm can reach 150 km (azimuth) × 150 km (range). Simulation parameters are listed in Table D1. The simulated swath is portrayed in Figure D1, and the swath center is at the equator.   Figure D1. T5 is at swath center. T1 is 75 km away from T5 along the azimuth and range directions.
Range and azimuth profiles corresponding to every target are illustrated in Figure D2, and evaluation results are listed in Table D2. As shown in Table D2, the difference of range PSLRs, Figure D1. T 5 is at swath center. T 1 is 75 km away from T 5 along the azimuth and range directions.
Range and azimuth profiles corresponding to every target are illustrated in Figure D2, and evaluation results are listed in Table D2. As shown in Table D2, the difference of range PSLRs, azimuth PSLRs, range ISLRs, and azimuth ISLRs over the whole swath is ď0.02 dB, ď0.06 dB, ď0.03 dB, and ď0.09 dB respectively. Some broadening coefficients are less than 1, because the modulation rates corresponding to these points are changed in imaging, and the resolutions become better. The simulation results show that this algorithm can achieve good imaging quality over the swath of 150 kmˆ150 km. azimuth PSLRs, range ISLRs, and azimuth ISLRs over the whole swath is ≤0.02 dB, ≤0.06 dB, ≤0.03 dB, and ≤0.09 dB respectively. Some broadening coefficients are less than 1, because the modulation rates corresponding to these points are changed in imaging, and the resolutions become better. The simulation results show that this algorithm can achieve good imaging quality over the swath of 150 km × 150 km.