An Accurate GEO SAR Range Model for Ultralong Integration Time Based on m th-Order Taylor Expansion

: As the Geosynchronous Earth Orbital Synthetic Aperture Radar (GEO SAR) allows a wide area viewing combined with a short revisit cycle, it is suitable for many applications that require high timeliness, such as natural disaster monitoring, weather supervision, and military reconnaissance. However, the ultralong integration time and the invalidation of “stop-and-go” assumption caused by the raise of orbital height also greatly increase the difﬁculty of signal processing. In this paper, a generalized method for calculating the accurate propagation distance between a GEO satellite and a target with ultralong integration time is proposed. This range model is mainly composed of an accurate pulse transmitting distance and an error compensation term for “stop-and-go” assumption failure. The transmitting distance is obtained by Taylor expansion, and the speciﬁc derivation process of the general formula of the m th-order expansion is given, in this paper. As for the compensation term, this is achieved by approximately calculating the pulse receiving distance based on twice Taylor expansion, the ﬁrst expansion is for fast-time and the other is for slow-time. Finally, a series of simulation experiments were conducted to verify the effectiveness and superiority of this new range model for an ultralong integration time.


Introduction
As an advanced active aerospace microwave remote sensing method, Synthetic Aperture Radar (SAR) can provide high resolution images in all-weather and all-time, which is widely used in military reconnaissance and many civil fields [1]. However, with the demand for shorter revisit cycles and larger coverage of imaging, the existing Low Earth Orbital SAR (LEO SAR) can hardly satisfy people's requirements. In addition, the rapid development of space ballistic missiles and laser weapons in recent years also put a severe test on the battlefield survivability of LEO satellites. Therefore, increasing the orbital altitude of SAR satellites becomes a future development trend.
The concept of Geosynchronous Earth Orbital SAR (GEO SAR) was first proposed by Kiyo Tomiyasu in 1978 [2,3], which indicates that GEO SAR satellite is in the orbit at an altitude of approximately 36,500 km and the period is almost one Earth day. Compared to the conventional SAR satellite located in lower earth orbit, it has many advantages, such as short revisit cycle, large observe swath, and instantaneous accessibility of targets on the Earth surface. However, limited by the antenna technology, rocket carrying capacity, and satellite performance at that time, research work on GEO SAR remained stagnant since it was proposed. Until the end of the 20th century, with the improvement of many key technologies, GEO SAR research turned active again. an accurate range model, named the Doppler parameter range model (DRM), to perform signal processing for GEO SAR, and the simulation results validated that DRM has a satisfactory effect on an L-band GEO SAR system with an azimuth resolution around 5 m [32,33]. Although this model achieved good fitting results for 660 s exposure time, the invalidation of "stop-and-go" assumption was not considered. In 2014, Tian et al. proposed an "Equivalent Position" model, which considered the range error introduced by the "stopand-go" assumption [34]. However, this equivalent model was derived on the basis of a circular orbit, which limited the usage of this model. Liu et al. presented an iterative approximation range model, which also took the assumption failure into account and had high computational efficiency and accuracy [35]. However, this model was obtained by iterative calculation, so it could not be expressed by analytical expression and further spectrum derivation.
From the above, the existing GEO SAR slant range models still have the following problems that are not yet resolved: (1) As the azimuth resolution increases, the low-order Taylor expansion range model shows insufficient accuracy in ultralong integration time. (2) The iterative approximation range model has a high accuracy but cannot be expressed by an analytical expression.
(3) The range models lack the flexibility of precision adjustment for different exposure times.
The above problems are what this paper is dedicated to solving. This paper is organized as follows. In Section 2, there is a brief explanation of the invalidation of "stop-andgo" assumption, then the integration time of GEO SAR at different azimuth resolutions and fitting effects of several "non-stop-and-go" range models are shown. In Section 3, a novel range model based on the mth-order Taylor expansion is put forward, and the derivation of the pulse transmitting and receiving distance are given, respectively. Section 4 shows the simulation results of fitting errors and imaging results, based on several slant range models, which validate the effectiveness of the proposed range model in this paper. Section 5 gives conclusions.

Invalidation of "Stop-and-Go" Assumption
The "stop-and-go" assumption is widely used in classic SAR imaging algorithms, which assumes that the radar platform and targets are relatively stationary, during pulse propagation. According to this hypothesis, the total propagation distance can be calculated as twice the instantaneous slant range between the antenna phase center and the target, at the time of pulse transmission. Nevertheless, this assumption ignores the range errors introduced by two forms of motion, one is the movement of platform relative to the target during pulse transmission and reception, and the other is the movement of platform relative to the target during pulse propagation. Figure 1 illustrates the GEO SAR pulse propagation model in an Earth-centered Earth-fixed (ECEF) coordinate system.
The first kind of error is determined by the radar system, which can be ignored when the pulse width is short enough. In GEO SAR signal processing, the second movement has a greater impact on the calculation of the slant range between the radar and target, especially when the integration time is very long. The propagation delays of the geosynchronous orbit satellite are always longer than 200 ms, during this period, the satellite runs along the orbit and the target is also displaced due to the rotation of the earth and its own motion. Thus, the range error caused by the "stop-and-go" assumption cannot be ignored. GEO SAR pulse propagation model in the ECEF coordinate system. P is a point target on the earth, S0 is the satellite position when the pulse is transmitted, S1 is the satellite position when the pulse reaches the target, S2 is the satellite position when the pulse is received. r1 and r2 represent the propagation distance of pulse transmission and reception, respectively.
The first kind of error is determined by the radar system, which can be ignored when the pulse width is short enough. In GEO SAR signal processing, the second movement has a greater impact on the calculation of the slant range between the radar and target, especially when the integration time is very long. The propagation delays of the geosynchronous orbit satellite are always longer than 200 ms, during this period, the satellite runs along the orbit and the target is also displaced due to the rotation of the earth and its own motion. Thus, the range error caused by the "stop-and-go" assumption cannot be ignored.

Ultralong Integration Time
The azimuth resolution of SAR is determined by the angle of synthetic aperture, so a SAR system always needs a longer integration time to obtain the same resolution on a higher orbit. A classic formula for calculating the integration time of different azimuth resolution can be expressed as follows: where    c R denotes the slant range at the center of the aperture,  syn is the angle of synthetic aperture, and ecef v is the velocity of the platform in the ECEF coordinate system.
Unlike LEO SAR and airborne SAR, the velocity of GEO SAR satellite is quite different in the Earth-centered inertial (ECI) coordinate and the ECEF coordinate. This is mainly because the velocity of the satellite decreases as the orbital height increases, which causes the effect of the earth's rotation to become more and more obvious. However, Equation (1) is inaccurate for GEO SAR because of the complicated space trajectory and velocity changes. To illustrate this problem, we simulated the GEO SAR space trajectory and velocity under two orbital configurations. The simulation parameters are shown in Table 1. GEO SAR pulse propagation model in the ECEF coordinate system. P is a point target on the earth, S 0 is the satellite position when the pulse is transmitted, S 1 is the satellite position when the pulse reaches the target, S 2 is the satellite position when the pulse is received. r 1 and r 2 represent the propagation distance of pulse transmission and reception, respectively.

Ultralong Integration Time
The azimuth resolution of SAR is determined by the angle of synthetic aperture, so a SAR system always needs a longer integration time to obtain the same resolution on a higher orbit. A classic formula for calculating the integration time of different azimuth resolution can be expressed as follows: where R(η c ) denotes the slant range at the center of the aperture, θ syn is the angle of synthetic aperture, and v ece f is the velocity of the platform in the ECEF coordinate system. Unlike LEO SAR and airborne SAR, the velocity of GEO SAR satellite is quite different in the Earth-centered inertial (ECI) coordinate and the ECEF coordinate. This is mainly because the velocity of the satellite decreases as the orbital height increases, which causes the effect of the earth's rotation to become more and more obvious. However, Equation (1) is inaccurate for GEO SAR because of the complicated space trajectory and velocity changes. To illustrate this problem, we simulated the GEO SAR space trajectory and velocity under two orbital configurations. The simulation parameters are shown in Table 1. In Table 1, the orbital inclination of 53° and the eccentricity of 0.07 corresponds to the classic "8" shaped orbit, and the inclination of 7.4° and the eccentricity of 0.1 corresponds to a near-circular orbit, which is introduced in [12]. The simulation results of trajectory and velocity under these two orbit configurations are shown in Figures 2 and 3. It could be seen that the trajectory of GEO SAR in the ECEF coordinate always presents to be a complex spatial curve, and the shape of the trajectory is determined by the orbital parameters. Common shapes include "8" shaped, drop-shaped, and ellipse. From Figure 3, the velocity of GEO SAR in the ECEF coordinate is much lower than that of LEO SAR, which means that the integration time is seriously increased. In Table 1, the orbital inclination of 53° and the eccentricity of 0.07 corresponds to the classic "8" shaped orbit, and the inclination of 7.4° and the eccentricity of 0.1 corresponds to a near-circular orbit, which is introduced in [12]. The simulation results of trajectory and velocity under these two orbit configurations are shown in Figures 2 and 3. It could be seen that the trajectory of GEO SAR in the ECEF coordinate always presents to be a complex spatial curve, and the shape of the trajectory is determined by the orbital parameters. Common shapes include "8" shaped, drop-shaped, and ellipse. From Figure 3, the velocity of GEO SAR in the ECEF coordinate is much lower than that of LEO SAR, which means that the integration time is seriously increased. Considering that (1) will introduce an intolerable error when used for calculating the integration time of GEO SAR, we propose an accurate search method based on dichotomy. The basic steps can be briefly described as follows: 1 Set an initial search region [Tmin, Tmax] according to the estimated integration time by (1), and this region is required to include the real integration time. Through the above steps, an accurate integration time can be obtained, its accuracy is controlled by the interval threshold, which can be selected according to the error requirements. The interval threshold is selected as 1 s in this experiment, and the simulation results are shown in Figure 4. It could be seen that the trajectory of GEO SAR in the ECEF coordinate always presents to be a complex spatial curve, and the shape of the trajectory is determined by the orbital parameters. Common shapes include "8" shaped, drop-shaped, and ellipse. From Figure 3, the velocity of GEO SAR in the ECEF coordinate is much lower than that of LEO SAR, which means that the integration time is seriously increased.
Considering that (1) will introduce an intolerable error when used for calculating the integration time of GEO SAR, we propose an accurate search method based on dichotomy. The basic steps can be briefly described as follows:

1.
Set an initial search region [T min , T max ] according to the estimated integration time by (1), and this region is required to include the real integration time.

2.
Define T mid = (T min + T max )/2. If T max − T min is smaller than the threshold, then return the T real = T mid and stop the iteration.
Through the above steps, an accurate integration time can be obtained, its accuracy is controlled by the interval threshold, which can be selected according to the error requirements. The interval threshold is selected as 1 s in this experiment, and the simulation results are shown in Figure 4.
Considering that (1) will introduce an intolerable error when used for calculating the integration time of GEO SAR, we propose an accurate search method based on dichotomy. The basic steps can be briefly described as follows: 1 Set an initial search region [Tmin, Tmax] according to the estimated integration time by (1), and this region is required to include the real integration time. 2 Define Tmid = (Tmin + Tmax)/2. If Tmax − Tmin is smaller than the threshold, then return the Treal = Tmid and stop the iteration. 3 Calculate the synthetic aperture angle of different exposure time Tmin, Tmax, and Tmid, denoted as min, max, and mid.

4
Define the new search region: if mid > syn, then Tmax = Tmid, else let Tmin = Tmid. Return to step 2.
Through the above steps, an accurate integration time can be obtained, its accuracy is controlled by the interval threshold, which can be selected according to the error requirements. The interval threshold is selected as 1 s in this experiment, and the simulation results are shown in Figure 4.  According to Figure 4, the integration time of the classic "8" shaped orbit is shorter than that of the near-circular orbit, but the maximum integration time also reaches 1086 s at 5 m azimuth resolution. As for the near-circular orbit, the maximum integration time even reaches 2400 s at the same resolution. Regarding 2 m azimuth resolution, the integration time becomes much longer. Therefore, whether the existing range models can be applied to such a long imaging time remains to be verified.

"Non-Stop-and-Go" Range Models
With the continuous deepening of GEO SAR research, many scholars proposed their own range models, which took the "stop-and-go" assumption failure into account [31][32][33][34][35]. According to Figure 4, the integration time of the classic "8" shaped orbit is shorter than that of the near-circular orbit, but the maximum integration time also reaches 1086 s at 5 m azimuth resolution. As for the near-circular orbit, the maximum integration time even reaches 2400 s at the same resolution. Regarding 2 m azimuth resolution, the integration time becomes much longer. Therefore, whether the existing range models can be applied to such a long imaging time remains to be verified.

"Non-Stop-and-Go" Range Models
With the continuous deepening of GEO SAR research, many scholars proposed their own range models, which took the "stop-and-go" assumption failure into account [31][32][33][34][35]. In order to analyze the applicability of these models for long integration times, a series of fitting experiments were conducted. The parameters used in simulation are shown in Table 1, the observation time was chosen as 2000 s, and the "8" shaped orbit was selected for this experiment. The phase error simulation results of the whole orbit are shown in Figure 5 and Table 2. Among these figures, x-axis is the observation time, y-axis is the true anomaly, and z-axis is the phase error between these range models and the real propagation distance.
In order to get the fitting error of these range models, we also need to know the real distance of pulse propagation. Considering that the equation about the real distance is a transcendental equation and has no analytical solution, a dichotomy search method is used to approximate it. The specific calculation process is described in [36].
In order to analyze the applicability of these models for long integration times, a series of fitting experiments were conducted. The parameters used in simulation are shown in Table 1, the observation time was chosen as 2000 s, and the "8" shaped orbit was selected for this experiment. The phase error simulation results of the whole orbit are shown in Figure  5 and Table 2. Among these figures, x-axis is the observation time, y-axis is the true anomaly, and z-axis is the phase error between these range models and the real propagation distance.
(a) (b) (c)  In order to get the fitting error of these range models, we also need to know the real distance of pulse propagation. Considering that the equation about the real distance is a transcendental equation and has no analytical solution, a dichotomy search method is used to approximate it. The specific calculation process is described in [36].
According to Figure 5, a simple summary of these two "non-stop-and-go" range models can be made: (1) Fourth-order Taylor expansion model considering the invalidation of "stop-and-go" assumption. This model can be expressed in the form of a power series, which can obtain the 2-D spectrum to form an efficient frequency domain imaging algorithm. However, limited by the order of expansion, this model only has a satisfactory fitting effect in a short observation time. As the exposure time increases, the fitting error accumulates rapidly at both edges of the aperture. (2) Iterative approximation range model. The main idea of this range model is to use twice the transmitting delay to approximately calculate the satellite position in the orbit, then an iteration can be formed to calculate a relatively accurate result. The calculation of this model is simple, and only one iteration already has a high accuracy, which can be used in echo generation and time-domain imaging algorithms. However, this range model cannot be expressed by an analytical expression, so its usage for frequency domain imaging algorithms is limited.

Proposed Range Model
Based on the two "non-stop-and-go" GEO SAR range models introduced in Section 2, an improved range model is proposed, which mainly has the following two advantages:  According to Figure 5, a simple summary of these two "non-stop-and-go" range models can be made: (1) Fourth-order Taylor expansion model considering the invalidation of "stop-and-go" assumption. This model can be expressed in the form of a power series, which can obtain the 2-D spectrum to form an efficient frequency domain imaging algorithm. However, limited by the order of expansion, this model only has a satisfactory fitting effect in a short observation time. As the exposure time increases, the fitting error accumulates rapidly at both edges of the aperture. (2) Iterative approximation range model. The main idea of this range model is to use twice the transmitting delay to approximately calculate the satellite position in the orbit, then an iteration can be formed to calculate a relatively accurate result. The calculation of this model is simple, and only one iteration already has a high accuracy, which can be used in echo generation and time-domain imaging algorithms. However, this range model cannot be expressed by an analytical expression, so its usage for frequency domain imaging algorithms is limited.

Proposed Range Model
Based on the two "non-stop-and-go" GEO SAR range models introduced in Section 2, an improved range model is proposed, which mainly has the following two advantages: (1) A general calculation method for mth-order expansion of pulse transmitting distance is given, and the expansion order can be adjusted according to the integration time and error accuracy requirement. (2) An accurate pulse receiving distance is obtained by using the thought of iterative approximation, and the analytical expression is obtained by Taylor expansion in the ECEF coordinate system.

Pulse Transmitting Distance
This part introduces the mth-order Taylor expansion method of the pulse transmitting distance, and all the following derivation is completed in the ground-fixed coordinate system. Without loss of generality, assume that the ground target has a moving speed. Figure 6 shows the propagation history of the transmitted pulse at the nth pulse repetition time (PRT).

Pulse Transmitting Distance
This part introduces the mth-order Taylor expansion method of the pulse transmitting distance, and all the following derivation is completed in the ground-fixed coordinate system. Without loss of generality, assume that the ground target has a moving speed. Figure 6 shows the propagation history of the transmitted pulse at the nth pulse repetition time (PRT). According to Figure 6, the pulse transmitting distance vector can be expressed as: To derive the modulus of r1n, we can first square the two sides of Equation (2), Taking the derivative of Equation (3), the first-derivative and second-derivative of the transmitting distance can be obtained: where Asn and Atn are the acceleration of satellite and target, respectively. According to Figure 6, the pulse transmitting distance vector can be expressed as: To derive the modulus of r 1n , we can first square the two sides of Equation (2), Taking the derivative of Equation (3), the first-derivative and second-derivative of the transmitting distance can be obtained: where A sn and A tn are the acceleration of satellite and target, respectively. Then, the high order derivatives of transmitting distance can be achieved by taking the derivative of Equation (5) step by step. In order to make the formula concise and easy to find the rules, the results of previous derivative are substituted into the next expression. The derivatives from the third-order to the fifth-order can be expressed as follows: 1n · r 1n + 3(r 1n ) 2 r 1n (7) r 1n 1n ·r 1n +10r (3) 1n ·r 1n r 1n (8) Remote Sens. 2021, 13, 255 9 of 26 According to the above equations and binomial expansion, the general formula of the mth-order derivative of the transmitting distance can be written in the following form: where C(m − 1, i) represents the combination formula. After obtaining the high order derivatives, the mth-order Taylor expansion of the transmitting distance can be performed at the beam center moment: where nT is the azimuth time, and k m represents the mth-order coefficient of azimuth time.
In order to calculate the mth-order derivative of transmitting distance, first, we need to find the first-order to mth-order derivatives of R sn and R tn . Considering that the motion state of the target is known already, the key point to this problem lies in the derivative of the satellite position vector. During the calculation process, we found that the mthorder derivative of R sn has a form that is easy to represent and iteratively calculate in the orbital coordinate system. Therefore, the derivatives of R sn is first calculated in the orbital coordinate system, and are then rotated into the ECEF coordinate system. For convenience, the three coordinate systems used in the subsequent derivation are defined here-E v represents the orbital coordinate system, E o represents the ECI coordinate system, and E g represents the ECEF coordinate system.
According to the orbit equation, the satellite velocity vector in the orbit coordinate system E v can be expressed as: where f is the true anomaly, µ is the geocentric gravitational constant, µ = 3.986005 × 10 14 m 3 /s 2 , P and e are the semi-latus rectum and the eccentricity of the orbit. Taking sin f and cos f as a whole, the mth-order derivative of the satellite position vector can be expressed as: According to Equations (12) and (13), if the 1st to (m − 1)-th order derivatives of true anomaly are obtained, the mth-order derivative of the satellite position vector in E v can be obtained. The first-order to fifth-order derivatives of f are given in Appendix A.
Next is the coordinate rotation part. First, rotate the mth-order derivative of R s_v from E v to E o . Considering that the rotation matrix A ov is constant, it can be directly written as: Second, rotate the mth-order derivative of R s_o from E o to E g . However, the rotation matrix A go has the variable orbit time t s , so we also need to take the derivative of A go . The rotation matrix A go can be written as: where H g is the Greenwich hour angle, ω e is the angular velocity of the earth's rotation, t s is the Greenwich mean time. The mth-order derivative of A go can be expressed as: Therefore, the mth-order derivative of satellite position vector in E g can be written as follows: From the above, the complete calculation process of the mth-order Taylor expansion of transmitting distance is given, and Figure 7 summarizes the process of calculation.

Pulse Receiving Distance
According to the analysis and simulation results in Section 2, the "stop-and-go" hypothesis is invalid in GEO SAR. Therefore, it is inaccurate to approximate the pulse receiving distance with the transmitting distance. In view of the fact that the iterative approximation range model has a high accuracy but cannot be expressed by an analytical expression, a new calculation method of pulse receiving distance is proposed in this paper, which combines the Taylor expansion with one iteration approximation. Figure 8 shows the propagation history of transmitting and receiving pulse at the nth PRT in the ECEF coordinate system.

Pulse Receiving Distance
According to the analysis and simulation results in Section 2, the "stop-and-go" hypothesis is invalid in GEO SAR. Therefore, it is inaccurate to approximate the pulse receiving distance with the transmitting distance. In view of the fact that the iterative approximation range model has a high accuracy but cannot be expressed by an analytical expression, a new calculation method of pulse receiving distance is proposed in this paper, which combines the Taylor expansion with one iteration approximation. Figure 8 shows the propagation history of transmitting and receiving pulse at the nth PRT in the ECEF coordinate system. pothesis is invalid in GEO SAR. Therefore, it is inaccurate to approximate the pulse receiving distance with the transmitting distance. In view of the fact that the iterative approximation range model has a high accuracy but cannot be expressed by an analytical expression, a new calculation method of pulse receiving distance is proposed in this paper, which combines the Taylor expansion with one iteration approximation. Figure 8 shows the propagation history of transmitting and receiving pulse at the nth PRT in the ECEF coordinate system. The relationship between pulse transmitting distance and receiving distance can be expressed as: Figure 8. The propagation history of transmitting and receiving pulse at the nth PRT. P is a point target on the earth. r 1n , r 2n represent the transmitting distance vector and receiving distance vector, respectively. τ n is the total propagation time of the pulse.
The relationship between pulse transmitting distance and receiving distance can be expressed as: Similar to the transmitting distance, the first-order and second-order derivatives of the receiving distance can be expressed as: The general formula of the mth-order derivative of the pulse receiving distance can be written in the following form: Considering that the pulse propagation time τ n is subsecond, its higher order terms have little effect on the fitting error. Therefore, this paper only performs a second-order Taylor expansion of the receiving distance on fast time, and the subsequent simulations in Section 4 also prove that the second-order expansion can achieve sufficiently high accuracy.
Performing the fast-time Taylor expansion of the receiving distance at a pulse transmitting time, that is, τ equal to 0, the expansion can be written as: In view of the fact that the actual value of the total pulse propagation time τ n cannot be solved, the idea of one iteration approximation is employed here. The total propagation time is replaced by double the transmitting delay: It should be pointed out that this approximation step implies the upper limit of the fitting performance of this method, its fitting error would not be better than that of the iterative approximation range model.
Substituting the approximate propagation time into Equation (24), the first-order term can be expressed as: Then, the slow-time Taylor expansion of ∆r 1 is performed at the beam center moment: Considering that the first-order term is the main source of the fitting error, it should be expanded to a specified order, according to the accuracy requirements in actual use. In the following simulation of this paper, the fifth-order expansion of ∆r 1 proved to have a sufficiently high accuracy.
Next, repeating the same steps, the second-order term can be expressed as: Similar to the expansion of ∆r 1 , the slow-time Taylor expansion of ∆r 2 is performed at the beam center moment: Given that ∆r 2 is the second-order term of fast-time, it has less impact on the fitting error. Therefore, only the first-order expansion is sufficient to ensure the fitting accuracy. The calculation results of the ∆r 20 and ∆r 21 are as follows: From the above, after considering the failure of "stop-and-go" assumption, the pulse receiving distance can be expanded as: The first part of this section introduces a general method to calculate the mth-order Taylor expansion of pulse transmitting distance in GEO SAR, whose expansion order can be adjusted by the accuracy requirement and integration time. The second part illustrates an approach to approximately calculate the pulse receiving distance based on twice Taylor expansion, the first expansion is for fast-time and the other is for slow-time. In addition, all derivation of this range model is based on a ground moving target, so this model is also applicable to the moving target.
Based on the derivation and analysis above, an accurate range model in GEO SAR, which considers the invalidation of "stop-and-go" assumption can be expressed as: r 1n ≈ r 10 + k 1 (nT) + k 2 (nT) 2 + k 3 (nT) 3 + · · · + k m (nT) m ∆r ≈ 2 ∆r 10 + ∆r 20 + (∆k 11 + ∆k 21 )(nT) + (∆k 12 where r 1n and ∆r represents the Taylor expansion of the pulse transmitting distance and the compensation term for the "stop-and-go" error, respectively. The fitting precision of this range model for the ultralong integration time is mainly determined by r 1n , whose accuracy can be adjusted by the expansion order m. The calculation of ∆r mainly uses the ideas of iterative approximation and twice Taylor expansion, since the total propagation time is subsecond, there is often no need to expand ∆r to a very high order in actual use.

Imaging Method
Considering the curved trajectory and invalidation of the "stop-and-go" assumption of GEO SAR, the conventional imaging algorithms cannot be directly used for focusing. Based on the range model proposed in this paper, an improved back projection (BP) focusing algorithm is used for imaging [37,38]. Figure 9 shows the flowchart of the improved BP algorithm.
where r1n and r represents the Taylor expansion of the pulse transmitting distance and the compensation term for the "stop-and-go" error, respectively. The fitting precision of this range model for the ultralong integration time is mainly determined by r1n, whose accuracy can be adjusted by the expansion order m. The calculation of r mainly uses the ideas of iterative approximation and twice Taylor expansion, since the total propagation time is subsecond, there is often no need to expand r to a very high order in actual use.

Imaging Method
Considering the curved trajectory and invalidation of the "stop-and-go" assumption of GEO SAR, the conventional imaging algorithms cannot be directly used for focusing. Based on the range model proposed in this paper, an improved back projection (BP) focusing algorithm is used for imaging [37,38]. Figure 9 shows the flowchart of the improved BP algorithm. The main process of the algorithm is the same as that of traditional BP, and the differences are marked in the flowchart with red boxes. The main different steps are as follows: (1) Calculating the true anomaly and its derivatives from 1st to (m − 1)-th order based on the satellite ephemeris data and orbit configuration; (2) Determining the Taylor expansion order m by the SAR system parameters and error accuracy, and then performing the mth-order Taylor expansion to obtain the transmitting distance for each grid point. (3) Calculating the compensation term of the "stop-and-go" assumption failure, and then bringing the complete pulse propagation distance into the azimuth compression. The main process of the algorithm is the same as that of traditional BP, and the differences are marked in the flowchart with red boxes. The main different steps are as follows:

Simulation Results and Discussion
(1) Calculating the true anomaly and its derivatives from 1st to (m − 1)-th order based on the satellite ephemeris data and orbit configuration; (2) Determining the Taylor expansion order m by the SAR system parameters and error accuracy, and then performing the mth-order Taylor expansion to obtain the transmitting distance for each grid point. (3) Calculating the compensation term of the "stop-and-go" assumption failure, and then bringing the complete pulse propagation distance into the azimuth compression.

Simulation Results and Discussion
The range model proposed in this paper can be divided into two parts-pulse transmitting distance and compensation term for "stop-and-go" error. In order to verify the effectiveness of the proposed range model and analyze the impact of different terms on the fitting accuracy, this section gives some simulation results on pulse propagation distance fitting and point target focusing. The simulation parameters are shown in Table 1, and if there are none specifically mentioned below, the satellite orbit is an "8" shaped orbit. All targets are in the stationary state.

Fitting Error of the Pulse Transmitting Distance
This subsection focuses on simulating the fitting error of each order of Taylor expansion for different integration time, the simulation results can be helpful for selecting the expansion order of the range model. According to Equation (34), the final slant distance is composed of twice the transmitting distance and a compensation term, so the phase error of the transmitting distance should be at least less than π/8. Figure 10 and Table 3 show the fitting error of the pulse transmitting distance of the whole orbit from fourth-order to sixth-order Taylor expansion. According to the simulation results, it was obvious that the fitting error would gradually decrease as the expansion order increased. With an observation time of 2000 s, the maximum phase error of the fourth-order and fifth-order expansion were both greater than π/8, only the sixth-order Taylor expansion could meet the accuracy requirement.  Table 1, and if there are none specifically mentioned below, the satellite orbit is an "8" shaped orbit. All targets are in the stationary state.

Fitting Error of the Pulse Transmitting Distance
This subsection focuses on simulating the fitting error of each order of Taylor expansion for different integration time, the simulation results can be helpful for selecting the expansion order of the range model. According to Equation (34), the final slant distance is composed of twice the transmitting distance and a compensation term, so the phase error of the transmitting distance should be at least less than /8  . Figure 10 and Table 3 show the fitting error of the pulse transmitting distance of the whole orbit from fourth-order to sixth-order Taylor expansion. According to the simulation results, it was obvious that the fitting error would gradually decrease as the expansion order increased. With an observation time of 2000 s, the maximum phase error of the fourth-order and fifth-order expansion were both greater than /8  , only the sixth-order Taylor expansion could meet the accuracy requirement.   Furthermore, it could be seen from Figure 10 that the positions corresponding to the minimum and maximum fitting errors changed with the expansion orders, and showed obvious odd and even changes. This was mainly because the reminder terms corresponding to the different order expansion were different forms, so their extreme points would appear in different positions.

Expansion Order Mean Error (rad) Max Error (rad) Standard Deviation (rad)
In order to further analyze the selection of the expansion order with different integration times, the maximum fitting phase error of the whole orbit was simulated. The observation time was selected as 6000 s, and the results are shown in Figure 11.  Furthermore, it could be seen from Figure 10 that the positions corresponding to the minimum and maximum fitting errors changed with the expansion orders, and showed obvious odd and even changes. This was mainly because the reminder terms corresponding to the different order expansion were different forms, so their extreme points would appear in different positions.
In order to further analyze the selection of the expansion order with different integration times, the maximum fitting phase error of the whole orbit was simulated. The observation time was selected as 6000 s, and the results are shown in Figure 11. It was obvious that the observation time corresponding to the error bound would change as the satellite orbit configuration changed. Therefore, in order to analyze the influence of orbit configuration on the fitting error, we simulated both phase error in the two types of orbit. It should be pointed out that since the pitch angle of some positions in the near-circular orbit was too large, the Zero-Doppler attitude guidance in [39,40] was no longer applicable, and the attitude steering was only performed on the former.
According to the symmetry of Figure 11a, the observation time corresponding to the error bound /8  of the 3rd-order to 7th-order Taylor expansion were 328 s, 870 s, 1866 s, 3050 s, and 4744 s. As for Figure 11b, the observation time corresponding to error bound /8  were 516 s, 1146 s, 2180 s, 3646 s, and 5534 s, which were slightly larger than that of the "8" shaped orbit. This was mainly because the velocity of the satellite in the nearcircular orbit was smaller, so the accumulated error was less than that of the "8" shaped orbit, for the same duration.
From the above simulation results, in order to select an appropriate expansion order to reduce the computation cost, it was necessary to analyze the fitting errors of different Taylor expansion orders in advance.

Fitting Error of Compensation Term
The error of "stop-and-go" assumption mainly comes from replacing the pulse receiving distance with the transmitting distance. However, the equation on the real receiving distance was a transcendental equation and had no analytical solution, so it could only be represented using approximate methods. Considering that the compensation term was derived from iterative approximation and Taylor expansion, its fitting effect would not be better than that of the iterative approximation range model. Figures 12 and 13 show the fitting errors of these two methods for calculating the pulse receiving distance for different integration times. The expansion order used in the proposed method was the same as Equation (33), which is a fifth-order expansion for the first-order term and a first-order expansion for the second-order term. It was obvious that the observation time corresponding to the error bound would change as the satellite orbit configuration changed. Therefore, in order to analyze the influence of orbit configuration on the fitting error, we simulated both phase error in the two types of orbit. It should be pointed out that since the pitch angle of some positions in the near-circular orbit was too large, the Zero-Doppler attitude guidance in [39,40] was no longer applicable, and the attitude steering was only performed on the former.
According to the symmetry of Figure 11a, the observation time corresponding to the error bound π/8 of the 3rd-order to 7th-order Taylor expansion were 328 s, 870 s, 1866 s, 3050 s, and 4744 s. As for Figure 11b, the observation time corresponding to error bound π/8 were 516 s, 1146 s, 2180 s, 3646 s, and 5534 s, which were slightly larger than that of the "8" shaped orbit. This was mainly because the velocity of the satellite in the near-circular orbit was smaller, so the accumulated error was less than that of the "8" shaped orbit, for the same duration.
From the above simulation results, in order to select an appropriate expansion order to reduce the computation cost, it was necessary to analyze the fitting errors of different Taylor expansion orders in advance.

Fitting Error of Compensation Term
The error of "stop-and-go" assumption mainly comes from replacing the pulse receiving distance with the transmitting distance. However, the equation on the real receiving distance was a transcendental equation and had no analytical solution, so it could only be represented using approximate methods. Considering that the compensation term was derived from iterative approximation and Taylor expansion, its fitting effect would not be better than that of the iterative approximation range model. Figures 12 and 13 show the fitting errors of these two methods for calculating the pulse receiving distance for different integration times. The expansion order used in the proposed method was the same as Equation (33), which is a fifth-order expansion for the first-order term and a first-order expansion for the second-order term.
be represented using approximate methods. Considering that the compensation term was derived from iterative approximation and Taylor expansion, its fitting effect would not be better than that of the iterative approximation range model. Figures 12 and 13 show the fitting errors of these two methods for calculating the pulse receiving distance for different integration times. The expansion order used in the proposed method was the same as Equation (33), which is a fifth-order expansion for the first-order term and a first-order expansion for the second-order term. Comparing Figures 12a and 13a, the fitting error of the proposed compensation term was close to the iterative approximation range model for short integration time. When the observation time increased to 1000 s and 2000 s, the fitting error of the proposed method would be slightly larger than that of the iterative approximation range model. However, the errors were still in the magnitude of 10 −5 and 10 −4 rad, which verified the effectiveness of the proposed compensation term for a long integration time.
According to the above simulation results, it could be seen that the compensation term proposed in this paper had a satisfactory effect for compensating the errors introduced by the "stop-and-go" assumption. Furthermore, considering that the motion of the target was taken into account when deriving the compensation term, it was also applicable to moving targets.

Focusing Results at Different Orbit Position
Since the motion state of the GEO satellite varied during the flight, the Taylor expansion range model had different fitting effects at different orbit positions. Figure 14 shows the relationship between the fitting error of the fourth-order Taylor expansion and the orbit positions; the observation time was 1000 s. Comparing Figures 12a and 13a, the fitting error of the proposed compensation term was close to the iterative approximation range model for short integration time. When the observation time increased to 1000 s and 2000 s, the fitting error of the proposed method would be slightly larger than that of the iterative approximation range model. However, the errors were still in the magnitude of 10 −5 and 10 −4 rad, which verified the effectiveness of the proposed compensation term for a long integration time.
According to the above simulation results, it could be seen that the compensation term proposed in this paper had a satisfactory effect for compensating the errors introduced by the "stop-and-go" assumption. Furthermore, considering that the motion of the target was taken into account when deriving the compensation term, it was also applicable to moving targets.

Focusing Results at Different Orbit Position
Since the motion state of the GEO satellite varied during the flight, the Taylor expansion range model had different fitting effects at different orbit positions. Figure 14 shows the relationship between the fitting error of the fourth-order Taylor expansion and the orbit positions; the observation time was 1000 s.
According to Figure 14, the maximum fitting error of the fourth-order expansion appeared at 45 • and 315 • of true anomaly, the minimum fitting error appeared in multiple locations, including the perigee and apogee. Therefore, to illustrate the focusing effect at different orbit positions, we conducted the imaging simulations based on fourth-order expansion at perigee and 45 • of true anomaly. The simulation parameters are shown in Table 1, and the orbit was chosen as the "8" shaped orbit, the integration time was chosen as 1000 s.

Focusing Results at Different Orbit Position
Since the motion state of the GEO satellite varied during the flight, the Taylor expansion range model had different fitting effects at different orbit positions. Figure 14 shows the relationship between the fitting error of the fourth-order Taylor expansion and the orbit positions; the observation time was 1000 s. According to Figure 14, the maximum fitting error of the fourth-order expansion appeared at 45° and 315° of true anomaly, the minimum fitting error appeared in multiple locations, including the perigee and apogee. Therefore, to illustrate the focusing effect at  Figure 15 shows the fitting error and the focusing results of the fourth-order Taylor expansion at perigee. From Figure 15a, it can be seen that the fitting error of fourthorder Taylor expansion was much smaller than π/4, which could lead to a satisfactory focusing results in Figure 15b,c. There were three quality parameters in Figure 15cimpulse response width (IRW), peak sidelobe ratio (PSLR), integrated sidelobe ratio (ISLR). According to the satellite synthetic aperture angle, the theoretical azimuth resolution was 2.2641 m, which was close to the experimental result.  Table 1, and the orbit was chosen as the "8" shaped orbit, the integration time was chosen as 1000 s. Figure 15 shows the fitting error and the focusing results of the fourth-order Taylor expansion at perigee. From Figure 15a, it can be seen that the fitting error of fourth-order Taylor expansion was much smaller than / 4 π , which could lead to a satisfactory focusing results in Figure 15b,c. There were three quality parameters in Figure 15c-impulse response width (IRW), peak sidelobe ratio (PSLR), integrated sidelobe ratio (ISLR). According to the satellite synthetic aperture angle, the theoretical azimuth resolution was 2.2641 m, which was close to the experimental result.  Figure 16 shows the fitting error and focusing result of the fourth-order Taylor expansion at 45° of true anomaly. From Figure 16b,c, it can be seen that the azimuth focusing results of the fourth-order Taylor expansion show obvious asymmetry. The main reason for azimuth-defocusing was that the radial velocity of satellite changed rapidly around this orbit position, and the fourth-order Taylor expansion was difficult to meet the accuracy requirement. According to Figure 16a, the fitting error showed obvious odd function characteristics, which was consistent with the fact that the main term of its reminder was a fifth-order polynomial. The maximum error was about 1.57 rad, which was larger than the error bound / 4 π , thus resulting in the defocusing of azimuth.  Figure 16 shows the fitting error and focusing result of the fourth-order Taylor expansion at 45 • of true anomaly. From Figure 16b,c, it can be seen that the azimuth focusing results of the fourth-order Taylor expansion show obvious asymmetry. The main reason for azimuth-defocusing was that the radial velocity of satellite changed rapidly around this orbit position, and the fourth-order Taylor expansion was difficult to meet the accuracy requirement. According to Figure 16a, the fitting error showed obvious odd function characteristics, which was consistent with the fact that the main term of its reminder was a fifth-order polynomial. The maximum error was about 1.57 rad, which was larger than the error bound π/4, thus resulting in the defocusing of azimuth.
According to Figure 17a, the fitting error of sixth-order Taylor expansion remained at a very low level at the same position, which could achieve relatively satisfactory focusing results. The imaging results are shown in Figure 17b,c. It should be noted that since the satellite runs faster at this orbit position, this position would obtain a larger synthetic aperture angle and a higher azimuth resolution for the same integration time. The theoretical azimuth resolution was 1.6069 m, which was also close to the experimental result.
for azimuth-defocusing was that the radial velocity of satellite changed rapidly around this orbit position, and the fourth-order Taylor expansion was difficult to meet the accuracy requirement. According to Figure 16a, the fitting error showed obvious odd function characteristics, which was consistent with the fact that the main term of its reminder was a fifth-order polynomial. The maximum error was about 1.57 rad, which was larger than the error bound /4  , thus resulting in the defocusing of azimuth. According to Figure 17a, the fitting error of sixth-order Taylor expansion remained at a very low level at the same position, which could achieve relatively satisfactory focusing results. The imaging results are shown in Figure 17b,c. It should be noted that since the satellite runs faster at this orbit position, this position would obtain a larger synthetic aperture angle and a higher azimuth resolution for the same integration time. The theoretical azimuth resolution was 1.6069 m, which was also close to the experimental result. Through the above simulation results, it could be seen that since the motion state of satellite changes all the time, the range model of a certain expansion order in different orbit positions would have different fitting effects. The perigee and apogee are always positions with the minimum fitting error, so it is insufficient to verify models or algorithms at only those positions. In conclusion, the expansion order should be determined by the maximum fitting error of the whole orbit, or using different expansion order at different orbit positions for the purpose of reducing the computational burden.

Focusing Results with Ultralong Integration Time
In order to verify the effectiveness of the mth-order Taylor expansion for target focusing with ultralong integration time, a series of simulation experiments for point targets imaging were conducted. The simulation parameters are shown in Table 1, and the orbit was chosen as the "8" shaped orbit. The integration time was 2000 s in these experiments, and the expansion order of pulse transmitting distance was sixth. Figure 18 shows the relationship between the fitting error of the sixth-order Taylor expansion and orbit position for 2000 s observation time. The imaging positions chosen are 55° of true anomaly and perigee, which were the locations of the maximum and minimum fitting error. The simulation scene was 40 km × 40 km, which had twenty-five point targets established in it, and the intervals between the two adjacent points were both 10 km in the range and azimuth direction. The distribution of point targets in the latitude-longitude grid at perigee is shown in Figure 19. Through the above simulation results, it could be seen that since the motion state of satellite changes all the time, the range model of a certain expansion order in different orbit positions would have different fitting effects. The perigee and apogee are always positions with the minimum fitting error, so it is insufficient to verify models or algorithms at only those positions. In conclusion, the expansion order should be determined by the maximum fitting error of the whole orbit, or using different expansion order at different orbit positions for the purpose of reducing the computational burden.

Focusing Results with Ultralong Integration Time
In order to verify the effectiveness of the mth-order Taylor expansion for target focusing with ultralong integration time, a series of simulation experiments for point targets imaging were conducted. The simulation parameters are shown in Table 1, and the orbit was chosen as the "8" shaped orbit. The integration time was 2000 s in these experiments, and the expansion order of pulse transmitting distance was sixth. Figure 18 shows the relationship between the fitting error of the sixth-order Taylor expansion and orbit position for 2000 s observation time. The imaging positions chosen are 55 • of true anomaly and perigee, which were the locations of the maximum and minimum fitting error. The simulation scene was 40 km × 40 km, which had twenty-five point targets established in it, and the intervals between the two adjacent points were both 10 km in the range and azimuth direction. The distribution of point targets in the latitude-longitude grid at perigee is shown in Figure 19.
relationship between the fitting error of the sixth-order Taylor expansion and orbit position for 2000 s observation time. The imaging positions chosen are 55° of true anomaly and perigee, which were the locations of the maximum and minimum fitting error. The simulation scene was 40 km × 40 km, which had twenty-five point targets established in it, and the intervals between the two adjacent points were both 10 km in the range and azimuth direction. The distribution of point targets in the latitude-longitude grid at perigee is shown in Figure 19.  It needs to be pointed out that these point targets were all placed on the earth surface according to satellite motions and the WGS84 coordinates, rather than being placed on the slant plane, which avoided the problem of targets separated from the earth surface. Targets from P1 to P5 were arranged along the range direction, and the beam scanning velocity was along the longitude at the perigee, according to Figure 19.
The focusing results of twenty-five targets are shown in Figure 20, the horizontal and vertical direction were the azimuth and range directions, respectively. Through the five refined imaging results, the targets located at different positions in the scene could all be well focused. Table 4 lists the evaluation measures of the five point targets, including PSLR, ISLR, and IRW. The positions in Table 4 indicate the coordinates of the point targets in the imaging scene, and assuming that the center-point was the origin of the coordinates. As shown in Table 4, the range PSLR and azimuth PSLR of different position targets were all close to the theoretical value −13.26 dB, and the ISLR of range direction were also close to the theoretical value −10 dB. It could be seen that the azimuth ISLR were slightly lower than the ideal value, which was mainly because of the curved trajectory of satellite. As for the IRW, the theoretical azimuth resolution of the center-point was 1.1346 m, which was in good agreement with experimental results. It needs to be pointed out that these point targets were all placed on the earth surface according to satellite motions and the WGS84 coordinates, rather than being placed on the slant plane, which avoided the problem of targets separated from the earth surface. Targets from P1 to P5 were arranged along the range direction, and the beam scanning velocity was along the longitude at the perigee, according to Figure 19.
The focusing results of twenty-five targets are shown in Figure 20, the horizontal and vertical direction were the azimuth and range directions, respectively. Through the five refined imaging results, the targets located at different positions in the scene could all be well focused. Table 4 lists the evaluation measures of the five point targets, including PSLR, ISLR, and IRW. The positions in Table 4 indicate the coordinates of the point targets in the imaging scene, and assuming that the center-point was the origin of the coordinates. As shown in Table 4, the range PSLR and azimuth PSLR of different position targets were all close to the theoretical value −13.26 dB, and the ISLR of range direction were also close to the theoretical value −10 dB. It could be seen that the azimuth ISLR were slightly lower than the ideal value, which was mainly because of the curved trajectory of satellite. As for the IRW, the theoretical azimuth resolution of the center-point was 1.1346 m, which was in good agreement with experimental results.  Next, performing the focusing experiment at the true anomaly of 55°. Figure 21 shows the distribution of point targets in the latitude-longitude coordinates. Similar to the previous experiment, the targets from P1 to P5 were arranged along the range direction. However, the difference was that the satellite velocity was no longer parallel to the longitude at this position, so the targets showed an oblique distribution in the latitudelongitude coordinates.  Next, performing the focusing experiment at the true anomaly of 55 • . Figure 21 shows the distribution of point targets in the latitude-longitude coordinates. Similar to the previous experiment, the targets from P1 to P5 were arranged along the range direction. However, the difference was that the satellite velocity was no longer parallel to the longitude at this position, so the targets showed an oblique distribution in the latitude-longitude coordinates.
The focusing results and evaluation parameters are shown in Figure 22 and Table 5. According to Table 5, most evaluation measures are close to the ideal values. Only the azimuth ISLR was slightly lower than −10 dB, which resulted from the cured trajectory of GEO satellite. The azimuth IRW was 0.76 m at this orbit position, which was close to the theoretical value of 0.7567 m. Remote Sens. 2021, 13, x FOR PEER REVIEW 22 of 27 The focusing results and evaluation parameters are shown in Figure 22 and Table 5. According to Table 5, most evaluation measures are close to the ideal values. Only the azimuth ISLR was slightly lower than −10 dB, which resulted from the cured trajectory of GEO satellite. The azimuth IRW was 0.76 m at this orbit position, which was close to the theoretical value of 0.7567 m.    The focusing results and evaluation parameters are shown in Figure 22 and Table 5. According to Table 5, most evaluation measures are close to the ideal values. Only the azimuth ISLR was slightly lower than −10 dB, which resulted from the cured trajectory of GEO satellite. The azimuth IRW was 0.76 m at this orbit position, which was close to the theoretical value of 0.7567 m.

The Influence of Different Sources of Distortion on Imaging
With an ultralong synthesis time and high resolution in GEO SAR, the sources of distortions in the phase distribution over the antenna synthesis aperture would have many forms. In this section, we briefly discuss the impact of two sources of distortion on imaging: (1) Trajectory distortions caused by errors in measuring the parameters of the spacecraft motion.
(2) Hardware distortions and instabilities during the formation and reception of probing signals (in particular, phase instability, instability of the microwave path parameters).
The simulation parameters are shown in Table 1, and the orbit was chosen as the "8" shaped orbit. The integration time was 2000 s in these experiments, and the orbital position was located at the perigee.
The first was the impact of motion parameter measurement errors. In order to simplify the model, we divided the form of measurement error into common forms such as constant value, quadratic polynomial, and sinusoidal function. The simulation results are shown in Figure 23.

The Influence of Different Sources of Distortion on Imaging
With an ultralong synthesis time and high resolution in GEO SAR, the sources of distortions in the phase distribution over the antenna synthesis aperture would have many forms. In this section, we briefly discuss the impact of two sources of distortion on imaging: (1) Trajectory distortions caused by errors in measuring the parameters of the spacecraft motion.
(2) Hardware distortions and instabilities during the formation and reception of probing signals (in particular, phase instability, instability of the microwave path parameters).
The simulation parameters are shown in Table 1, and the orbit was chosen as the "8" shaped orbit. The integration time was 2000 s in these experiments, and the orbital position was located at the perigee.
The first was the impact of motion parameter measurement errors. In order to simplify the model, we divided the form of measurement error into common forms such as constant value, quadratic polynomial, and sinusoidal function. The simulation results are shown in Figure 23. Through the above simulation experiments, it could be seen that the measurement error of satellite motion parameters would cause the target focusing distortion. The form of measurement errors determined the impact on imaging. Constant errors caused the target's offset, while quadratic polynomial and sinusoidal errors would cause the target to defocus in the azimuth direction.
Next, we discuss the impact of phase errors caused by hardware distortions and instabilities. Similar to the previous experiment, we divided the forms of phase distortion Through the above simulation experiments, it could be seen that the measurement error of satellite motion parameters would cause the target focusing distortion. The form of measurement errors determined the impact on imaging. Constant errors caused the target's offset, while quadratic polynomial and sinusoidal errors would cause the target to defocus in the azimuth direction.
Next, we discuss the impact of phase errors caused by hardware distortions and instabilities. Similar to the previous experiment, we divided the forms of phase distortion into common forms, such as constant value, quadratic polynomial, and sinusoidal function. The simulation results are shown in Figure 24. into common forms, such as constant value, quadratic polynomial, and sinusoidal function. The simulation results are shown in Figure 24. According to Figure 24, the phase error caused by hardware distortions and instabilities also had a certain impact on target focusing. Constant phase errors would not cause obvious changes in the image, while the quadratic polynomial and sinusoidal errors would cause varying degrees of defocus.
As a summary, a slight measurement error and hardware distortion would introduce an intolerant error with the ultralong integration time, which would result in deterioration of the images. Therefore, methods such as calibration and phase autofocus are needed to reduce their influence in the actual system.

Conclusions
In view of the wide application prospects, GEO SAR attracts more and more attention from scholars and related organizations. However, the extremely high orbit also brings great challenges to the establishment of accurate slant range models, especially for the range calculation for ultralong time. There are two main difficulties in solving this problem, one is the "stop-and-go" assumption failure caused by the long pulse propagation time, and the other is that the fitting error of traditional low-order Taylor expansion would quickly accumulate at both ends of the aperture.
To solve the problems induced by extremely high orbit and ultralong integration time, an accurate range model based on mth-order Taylor expansion was put forward in this paper. This model was mainly composed of an mth-order Taylor expansion of transmitting distance and an accurate compensation term of the error introduced by the invalidation of "stop-and-go" assumption. The superiority and effectiveness of this new model for ultralong integration time was verified through the simulation experiments. The imaging results demonstrated that this range model was effective for GEO SAR with 2000 s integration time and a 40 km wide imaging scene. The evaluation parameters of point targets were close to the ideal values, and the azimuth resolution could reach 1 meter. According to Figure 24, the phase error caused by hardware distortions and instabilities also had a certain impact on target focusing. Constant phase errors would not cause obvious changes in the image, while the quadratic polynomial and sinusoidal errors would cause varying degrees of defocus.
As a summary, a slight measurement error and hardware distortion would introduce an intolerant error with the ultralong integration time, which would result in deterioration of the images. Therefore, methods such as calibration and phase autofocus are needed to reduce their influence in the actual system.

Conclusions
In view of the wide application prospects, GEO SAR attracts more and more attention from scholars and related organizations. However, the extremely high orbit also brings great challenges to the establishment of accurate slant range models, especially for the range calculation for ultralong time. There are two main difficulties in solving this problem, one is the "stop-and-go" assumption failure caused by the long pulse propagation time, and the other is that the fitting error of traditional low-order Taylor expansion would quickly accumulate at both ends of the aperture.
To solve the problems induced by extremely high orbit and ultralong integration time, an accurate range model based on mth-order Taylor expansion was put forward in this paper. This model was mainly composed of an mth-order Taylor expansion of transmitting distance and an accurate compensation term of the error introduced by the invalidation of "stop-and-go" assumption. Acknowledgments: The authors would like to thank all reviewers and editors for their comments on this paper.

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

Appendix A
The orbit equation can be written as: where R s is the separation distance between the satellite and earth center, f is the true anomaly. P, e, and a are the semi-latus rectum, eccentricity, and semi-major axis of the orbit, respectively. The derivative of R s and f can be expressed as: Since Equation (A1) is an implicit equation, the mth-order derivative of f cannot be obtained directly, here are two solutions for this problem.
One way is using the chain rule for derivation with composite functions to derive. We can differentiate t on both sides of the equation, and after a series of derivation, the second-order to fifth-order derivatives of f can be expressed as: 12e 3 + 14e − 144e 3 + 32e sin 2 f + 25e 2 + 1 cos f − 163e 2 sin 2 f cos f + 192e 3 cos 4 f R 6 s (A7) Another way is to achieve the different order derivatives of f through numerical calculation, that is, to calculate the finite difference step-by-step. Experiments proved that this method has sufficient accuracy when the time intervals are short and the orders of derivative are not too high.