**1. Introduction**

Real-time synthetic aperture radar (SAR) imaging has always been a research focus in certain applications. Real-time imaging enhances the effectiveness and widens the application range of the SAR technique. In airborne SAR, real-time sensing plays an important role, and there are many well-established systems [1–5] and techniques [6–11]. Real-time processing techniques are also available for unmanned aerial vehicle (UAV) SAR applications [12–14]. For spaceborne SAR, because of its unique configuration and onboard computation capabilities, real-time SAR imaging systems and techniques are still being developed for both digital [15–17] and optical [18] methods.

For the purpose of implementing onboard real-time spaceborne SAR imaging, a set of Doppler parameters with high precision and accuracy must be calculated, and accurate position and velocity vectors are needed for the Doppler parameter calculation. In typical methods used to process SAR echo data, a high-accuracy orbit ephemeris is applied in the ground-based processing system, and state

vectors, i.e., position and velocity vectors of the SAR satellite, are determined with high accuracy. However, this kind of high-accuracy orbit ephemeris is only generated hours, days, or even weeks after a single SAR observation [19]. Thus, this kind of high-accuracy orbit ephemeris cannot be obtained during real-time SAR imaging. Nevertheless, we are still able to achieve onboard orbit determination using GNSS-based receivers and appropriate algorithms [20,21]. It is acknowledged that this approach to orbit determination has relatively low accuracy in the state vectors of SAR satellites compared with the high-accuracy orbit ephemeris. However, it remains the best, and possibly the only, onboard orbit determination method that is available for real-time SAR imaging.

SAR mission designers must account for the fact that errors in such onboard orbit determination data will adversely affect the quality of the SAR imaging product through inaccurate Doppler parameter calculation. Although estimation algorithms have been developed to generate the Doppler centroid and Doppler rate with high accuracy for SAR imaging, they are time-consuming in practice. Real-time spaceborne SAR imaging certainly does not belong on a mission that is not time-sensitive. The product of real-time SAR imaging may be the source data for other mission types, such as onboard deformation monitoring or target recognition. Given these requirements, the estimation algorithm should not be built into the processing workflow because of its long processing time. Thus, the concept of quadrature phase error (QPE) is introduced to evaluate the influence of onboard orbit determination errors on spaceborne SAR imaging products. The azimuth impulse response width of a SAR image has a certain relationship with QPE, and it provides a more intuitive evaluation standard of the errors.

Errors of an onboard orbit determination system are often given in the form of probabilities, along with the type of probability distributions and the corresponding parameters. The most common case for these errors is a normal distribution with an expected value of zero and a standard deviation *σ*. Therefore, the probability distribution of QPE needs to be examined. In practice, there exist two major methods that have been applied to QPE analysis, i.e., the extreme value method and the numerical simulation method. The extreme value method is the simplest way to investigate QPE; it includes the maximum possible error in the calculation and ignores the probability distribution of the errors. This method faces challenges in revealing the QPE probability distribution, which means it only reflects the extreme situation in practice. For the numerical simulation methods, Monte Carlo simulations are often applied [22]. However, this approach requires a huge number of generated error samples and repeated experiments.

In this paper, an analytical approximation model of the QPE introduced by orbit determination errors is proposed. With the a priori probability of the onboard orbit determination system and SAR satellite orbit elements, together with the observing geometry of the spaceborne SAR and appropriate approximations, a series of equations describing the parameters of the QPE probability distribution are presented. This analytical approximation model can provide approximation results at two granularities, i.e., the approximations with the satellite's true anomaly as the independent variable and the approximations for all positions of the satellite during its entire orbit.

Following this section, Section 2 describes the QPE analytical approximation model in general. Section 3 derives the key variables in the analytical approximation model and reveals the appropriate approximations in the model. This is followed by the evaluation of the analytical approximation model by comparing it with numerical simulation results in Section 4, and a related discussion is presented in Section 5. Finally, Section 6 provides a brief conclusion regarding the analytical approximation model.

#### **2. Proposed Model: QPE Analytical Approximation Model**

The QPE analytical approximation model requires four groups of parameters as input, i.e., the Earth's physical parameters, the satellite platform orbit parameters, the SAR payload key parameters, and the a priori probability of the onboard orbit determination system. In this section, the coordinate systems, vector in the line-of-sight method, and QPE approximation method in the proposed model are described.

#### *2.1. Coordinate Systems for Modeling*

Two main coordinate systems are used in the derivation of the analytical approximation model. These two coordinate systems can be helpful in generating a clear description of the relative movement of the spaceborne SAR and the point of interest on the surface of the Earth.

#### 2.1.1. Earth-Centered Inertial Coordinate System

Generally, the QPE calculation requires the position, velocity, and acceleration coordinates of the spaceborne SAR platform and the SAR antenna aiming point on the Earth's surface, together with other variables. A widely used coordinate system in modeling the relative motion between the spaceborne SAR and the point of interest is the Earth-Centered Inertial (ECI) coordinate system. This system has an *x*-axis and *z*-axis aligned with the mean equinox and the celestial North Pole, respectively. The *y*-axis forms a right-handed coordinate system that joins the *x*-axis and *z*-axis.

In common practice, the onboard GNSS-based orbit determination system of a spaceborne satellite calculates its positions in a given coordinate system using the GNSS broadcast message [21]. Each GNSS in service now has its own unique coordinate system, e.g., WGS 84 in GPS, CGCS2000 in BDS, and GTRF in Galileo. For modeling the orbit determination errors in the proposed analytical approximation model, all real-time-measured position and velocity coordinates of the spaceborne SAR platform are transformed into the ECI system in this paper.

#### 2.1.2. Orbital Plane Coordinate System

The movement parameters of a spaceborne SAR satellite can be described in the orbital plane coordinate system. In this system, the Earth is located at one of the foci of the elliptical orbit, and the satellite can be regarded as a point. Polar coordinates are used to describe the movement parameters.

There exist six orbital elements that uniquely identify a specific orbit of the satellite: the semi-major axis *a*, eccentricity *e*, and true anomaly *ν* are used in the orbital plane coordinate system to determine the satellite's movement at a certain time. The inclination *i*, the longitude of the ascending node Ω, and the argument of periapsis *ω* are necessary to construct the transformation matrix between the orbital plane coordinate system and the ECI system.

The orbital plane coordinate system and the ECI system are summarized in Figure 1.

**Figure 1.** The orbital plane (OP) coordinate system and the ECI system. Elements related to the OP coordinate system are in red. *ω* is the argument of periapsis, *i* is the orbit eccentricity, *ν* is true anomaly, and *θ<sup>L</sup>* is the center off-nadir angle in SAR.

#### *2.2. Vector in the Line-of-Sight Method*

In the analytical approximation model presented in this paper, the attitude errors of the spaceborne SAR satellite are not taken into consideration. The true value of the attitude of the platform is used in both the analytical approach and the following numerical verification.

With this assumption, the vector in the line-of-sight (VLS) method can be used to calculate the state vectors of the point of interest. The relative range from the SAR antenna phase center to the point of interest can be derived from the time delay of the echo signal. For spaceborne SAR geometry, the echo signal enters the receiver after several pulse repetition intervals (PRIs) until the moment of transmission, and then it is demodulated and sampled by the receiver. In each range profile of the sampled data, each sampled time delay can be determined by the index number of its range gate. The signal transmission delay from the antenna phase center to the sampling ADCs in the receiver can be accurately measured in the state-of-the-art spaceborne SAR system. *tdelay* represents the delayed time of the echo of the point-of-interest signal from its transmission time, and *c* represents the velocity of light, whose relative range is described by

$$
\rho = \frac{c t\_{delay}}{2} \tag{1}
$$

A unit line-of-sight vector **u**ˆ *<sup>a</sup>*,*ECI* in the ECI system can be created; it is perpendicular to the aperture of the SAR antenna and points to the surface of the Earth. The position vector **P***tar*,*ECI* and velocity vector **V***tar*,*ECI* of the point on the surface of the Earth at which the antenna points can be shown in stripmap SAR geometry as

$$\mathbf{P}\_{\text{tar},\text{EC},I} = \mathbf{P}\_{\text{sat},\text{EC},I} + \rho \cdot \hat{\mathbf{u}}\_{a,\text{EC},I} = \begin{bmatrix} P\_{\text{tar},\text{EC},\text{ $x$ }} \\ P\_{\text{tar},\text{EC},\text{ $y$ }} \\ P\_{\text{tar},\text{EC},\text{ $z$ }} \end{bmatrix} \tag{2}$$

$$\mathbf{V}\_{\rm tar,ECI} = \begin{bmatrix} -\omega\_{\rm c} \cdot P\_{\rm tar,ECI,y} \\ \omega\_{\rm c} \cdot P\_{\rm tar,ECI,x} \\ 0 \end{bmatrix} \tag{3}$$

in which *ω<sup>e</sup>* is the angular speed of the Earth's rotation.

Because there are errors in the satellite's position and velocity measurements from the onboard GNSS-based orbit determination system, the calculated **P***tar*,*ECI* and **V***tar*,*ECI* will also contain errors at the same time. The relative range *ρ* remains unaffected because errors exist in measurements while the real geometry remains unchanged, which means that the time delay remains equal to the accuracy value.

#### *2.3. QPE Approximation*

QPE can be expressed in radians as [23]

$$\text{QPE} = \pi \Delta f\_r \left(\frac{T\_{int}}{2}\right)^2 \tag{4}$$

in which *fr* is the Doppler rate of the echo signal, and *Tint* represents the integration time, also known as the synthetic-aperture time. In stripmap SAR geometry, the integration time is the time it takes the point of interest to enter and exit the entire 3 dB edge of the SAR antenna's beam illumination.

The variables **R***st*, **V***st*, and **A***st* represent the relative position, velocity, and acceleration of the spaceborne SAR satellite, respectively. In the ECI system, for the antenna aiming point on the ground, the commonly used *fr* expression in stripmap mode when the squint angle equals zero is

$$f\_r = \frac{2}{\lambda} \left[ \frac{\mathbf{V}\_{st} \cdot \mathbf{V}\_{st}}{\left| \mathbf{R}\_{st} \right|} + \frac{\mathbf{A}\_{st} \cdot \mathbf{R}\_{st}}{\left| \mathbf{R}\_{st} \right|} - \frac{\left( \mathbf{V}\_{st} \cdot \mathbf{R}\_{st} \right)^2}{\left| \mathbf{R}\_{st} \right|^3} \right] \tag{5}$$

The *fr* expression in Equation (5) is an approximation expression, and it can achieve high precision if more terms remain. In the analytical approximation model, only the first two terms in Equation (5) are analyzed. This simplifies the model while retaining a relatively high precision.

The orbit determination system measurements contain errors, which will lead to errors in **R***st*, **V***st*, and **A***st*. The terms that define the differences between the state vectors derived from measurements and the true state vectors are

$$
\Delta \mathbf{R}\_{\text{sf}} = \mathbf{R}\_{\text{st},\text{ct}} - \mathbf{R}\_{\text{sf}} \tag{6}
$$

$$
\Delta \mathbf{V}\_{st} = \mathbf{V}\_{st, \varepsilon} - \mathbf{V}\_{st} \tag{7}
$$

$$
\Delta \mathbf{A}\_{st} = \mathbf{A}\_{st, \varepsilon} - \mathbf{A}\_{st} \tag{8}
$$

where the subscript *e* indicates that the vectors contain measurement errors. Thus, the difference between the measured *fr* and true *fr* value is

$$
\Delta f\_r = f\_{r, \varepsilon} - f\_r \tag{9}
$$

The parameter *fr* in Equation (5) is often used to determine the precise value for SAR raw data simulation or imaging algorithms. For the approximate expression of QPE, considering the magnitude of each term in *fr*, some terms can be omitted while still achieving an acceptable approximation of QPE. The VLS method in Section 2.2 maintains the real value of the relative range *ρ* = |**R***st*|. Thus, Equation (4) combined with Equation (9) gives the QPE approximation expression

$$\begin{split} \text{QPE} &= \pi \left( f\_{r,t} - f\_r \right) \left( \frac{T\_{\text{int}}}{2} \right)^2 \\ &\approx \frac{2\pi}{\lambda \rho} \left( 2\Delta \mathbf{V}\_{st} \cdot \mathbf{V}\_{st} + \Delta \mathbf{A}\_{st} \cdot \mathbf{R}\_{st} \right) \left( \frac{T\_{\text{int}}}{2} \right)^2 \end{split} \tag{10}$$

For convenience, the terms Δ**V***st* · **V***st* and Δ**A***st* · **R***st* are named *velocity* and *acceleration vector term* in the following sections.

The approximation of QPE in Equation (10) is the basis of the QPE analytical approximation function and maximum QPE derivation.

#### **3. Key Variables in the Analytical Approximation Model**

For a given true anomaly *ν*, the distance between the satellite and the Earth's center is

$$r = \frac{a\left(1 - \varepsilon^2\right)}{1 + \varepsilon\cos\nu} = \frac{p}{1 + \varepsilon\cos\nu} \tag{11}$$

where *p* = *a* - <sup>1</sup> − *<sup>e</sup>*<sup>2</sup> in Equation (11) is the semi-latus rectum of the elliptic orbit. The position, velocity, and acceleration vectors of the satellite in the orbital plane (OP) are

$$\mathbf{P}\_{\text{sat}, OP} = r \begin{bmatrix} \cos \nu \\ \sin \nu \\ 0 \end{bmatrix} \tag{12}$$

$$\mathbf{V}\_{\text{sat}, OP} = \sqrt{\frac{\mu}{p}} \begin{bmatrix} -\sin\nu\\ \varepsilon + \cos\nu\\ 0 \end{bmatrix} \tag{13}$$

$$\mathbf{A}\_{\rm sat,OP} = -\frac{\mu \left(1 + e \cos \nu\right)^2}{p^2} \begin{bmatrix} \cos \nu\\ \sin \nu\\ 0 \end{bmatrix} \tag{14}$$

where *μ* is the standard gravitational parameter of the Earth.

The transformation matrix *AO*<sup>2</sup>*<sup>E</sup>* is used to convert the vector in the orbital plane to the ECI system:

$$\mathbf{A}\_{O2E} = \begin{bmatrix} \cos \omega & -\sin \omega & 0\\ \sin \omega & \cos \omega & 0\\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0\\ 0 & \cos i & -\sin i\\ 0 & \sin i & \cos i \end{bmatrix} \begin{bmatrix} \cos \Omega & -\sin \Omega & 0\\ \sin \Omega & \cos \Omega & 0\\ 0 & 0 & 1 \end{bmatrix} \tag{15}$$

The parameter Ω in **A***O*2*<sup>E</sup>* results in rotation with the *z*-axis in the ECI system. Considering the geometry and the Earth's ellipsoid reference, in the following discussion, we set Ω = 0 for convenience, and this will not result in an error value in the QPE approximation.

The onboard orbit determination system measurements are given in ECI terms, which can be regarded as the true state vector plus an error vector. With the assumption that the measurement errors of the satellite's position and velocity follow a normal distribution with an expected value of 0, the error vectors can be described as

$$
\Delta \mathbf{P} = \begin{bmatrix} \Delta p\_x \\ \Delta p\_y \\ \Delta p\_z \end{bmatrix} \sim \begin{bmatrix} N\left(0, \sigma\_{p,x}^2\right) \\ N\left(0, \sigma\_{p,y}^2\right) \\ N\left(0, \sigma\_{p,z}^2\right) \end{bmatrix} \tag{16}
$$

$$
\Delta \mathbf{V} = \begin{bmatrix} \Delta v\_x \\ \Delta v\_y \\ \Delta v\_z \end{bmatrix} \sim \begin{bmatrix} N\left(0, \sigma\_{v,x}^2\right) \\ N\left(0, \sigma\_{v,y}^2\right) \\ N\left(0, \sigma\_{v,z}^2\right) \end{bmatrix} \tag{17}
$$

Then, the state vector measurements of the satellite are

$$\mathbf{P}\_{\rm sat,ECI,\varphi} = \mathbf{A}\_{O2E} \mathbf{P}\_{\rm sat,OP} + \Delta \mathbf{P} \tag{18}$$

$$\mathbf{V}\_{\text{sat},ECI} = \mathbf{A}\_{\text{O2E}} \mathbf{V}\_{\text{sat},OP} + \Delta \mathbf{V} \tag{19}$$

From Equations (2), (3), (18) and (19), the position vectors of the antenna aiming point are

$$\mathbf{P}\_{\text{tar},ECI,\varepsilon} = \mathbf{P}\_{\text{sat},ECI,\varepsilon} + \rho \cdot \boldsymbol{\mathcal{d}}\_{a,ECI} \tag{20}$$

With the help of Equations (3), (19) and (20), the difference in the relative velocity from the satellite to the antenna aiming point can be given as

$$
\Delta \mathbf{V}\_{\rm st} = (\mathbf{V}\_{\rm sat, ECI, \varepsilon} - \mathbf{V}\_{\rm tar, ECI, \varepsilon}) - (\mathbf{V}\_{\rm sat, ECI} - \mathbf{V}\_{\rm tar, ECI}) = \begin{bmatrix} \Delta v\_{\rm x} + \omega\_{\rm t} \Delta p\_{\rm y} \\ \Delta v\_{\rm y} - \omega\_{\rm t} \Delta p\_{\rm x} \\ \Delta v\_{\rm z} \end{bmatrix} \tag{21}
$$

#### *3.1. The Velocity Vector Term in the Approximate Doppler Rate*

With Equation (21), we have Δ**V***st*. However, we still need the expression of **V***st* to calculate the velocity vector term Δ**V***st* · **V***st*. It is worth recalling that we need the approximate Doppler rate expression. To calculate the velocity vector term in Equation (10), we propose using the velocity vector of a satellite traveling in a circular orbit **V***<sup>s</sup>* in our analytical approximation model rather than using **V***st*. Thus, we have

$$\begin{split} \mathbf{V}\_{st} &= \mathbf{A}\_{O2E} \mathbf{V}\_{sat, OP}|\_{\varepsilon=0} \\ &= \sqrt{\frac{\mu}{p}} \begin{bmatrix} -\sin\left(\nu + \omega\right) \\ \cos\left(\nu + \omega\right)\cos i \\ \cos\left(\nu + \omega\right)\sin i \end{bmatrix} \end{split} \tag{22}$$

The velocity vector term results in a Doppler rate difference by

$$\begin{split} \Delta f\_{\tau,\nu} &= \frac{2}{\lambda \rho} \cdot 2\Delta \mathbf{V}\_{st} \cdot \mathbf{V}\_{st} \\ &= \frac{4}{\lambda \rho} \sqrt{\frac{\mu}{p}} \left[ - \left( \Delta v\_x + \omega\_t \Delta p\_y \right) \sin \left( \nu + \omega \right) \\ &\quad + \left( \Delta v\_y - \omega\_t \Delta p\_x \right) \cos \left( \nu + \omega \right) \cos i + \Delta v\_z \cos \left( \nu + \omega \right) \sin i \right] \end{split} \tag{23}$$

Assuming that the variances of the position and velocity errors are the same standard deviation in the *x*-, *y*-, and *z*-direction in the ECI system, i.e., *σp*,*<sup>x</sup>* = *σp*,*<sup>y</sup>* = *σp*,*<sup>z</sup>* = *σ<sup>p</sup>* and *σv*,*<sup>x</sup>* = *σv*,*<sup>y</sup>* = *σv*,*<sup>z</sup>* = *σv*, we can find the expected value and standard deviation of Δ*fr*,*<sup>v</sup>* by

$$\mathbb{E}\left[\Delta f\_{r,v}\right] = 0\tag{24}$$

$$\sigma \left[ \Delta f\_{r, \nu} \right] = \frac{4}{\lambda \rho} \sqrt{\frac{\mu}{p}} \sqrt{\sigma\_v^2 + \omega\_\varepsilon^2 \left[ \sin^2 \left( \nu + \omega \right) + \cos^2 \left( \nu + \omega \right) \cos^2 i \right]} \sigma\_p^2 \tag{25}$$

We also need an approximate maximum value of the standard deviation of Δ*fr*,*<sup>v</sup>* to calculate the approximate maximum value of QPE. We suppose that cos *i* has a maximum value of cos *i* = 1, and

$$
\sigma \left[ \Delta f\_{r,v} \right]\_{\max} = \frac{4}{\lambda \rho} \sqrt{\frac{\mu}{p}} \cdot \sqrt{\omega\_\varepsilon^2 \sigma\_p^2 + \sigma\_v^2} \tag{26}
$$

#### *3.2. True Anomaly*

Equation (21) gives the expression of difference vectors in the velocity vector term in Equation (10). However, the difference vector in the acceleration vector term, i.e., Δ**A***st*, is given its expression with the help of true anomaly calculation.

*Remote Sens.* **2019**, *11*, 1663

The accurate acceleration vector term in the orbital plane is given in Equation (14). True anomaly *ν* should be calculated using the onboard-measured satellite's position and velocity vectors. A common way to solve *ν* is given in [24]

$$\nu = \text{atan2}\left(\sqrt{\frac{p}{\mu}} \left(\mathbf{V}\_{\text{sat}, ECI} \cdot \mathbf{R}\_{\text{sat}, ECI}\right), p - |\mathbf{R}\_{\text{sat}, ECI}|\right) \tag{27}$$

where the atan2 (*Y*, *X*) function returns an unambiguous result of arctan (*Y*/*X*) with a range of (−*π*, *π*]. Obviously, the measured state vectors of the satellite **P***sat*,*ECI*,*<sup>e</sup>* and **V***sat*,*ECI*,*<sup>e</sup>* will result in a value of *ν<sup>e</sup>* that is not accurate. The probably distribution parameters of *ν<sup>e</sup>* should be determined in order to analyze the term Δ**A***st* and its influence on QPE.

The detailed derivation is in Appendix A. Here, the standard deviation and its maximum value of the difference between the calculated *ν<sup>e</sup>* and its accurate value Δ*ν* = *ν<sup>e</sup>* − *ν* are given.

$$\left| \sigma \left[ \Delta \nu \right] \right| = \left| \cos \nu \right| \cdot \frac{1 + \varepsilon \cos \nu}{\sqrt{\mu a \epsilon^2 \left( 1 - \epsilon^2 \right)}} \cdot \sqrt{\frac{\mu}{p} \left( 1 + \epsilon^2 + 2\varepsilon \cos \nu \right) \sigma\_p^2 + \frac{p^2}{\left( 1 + \varepsilon \cos \nu \right)^2} \sigma\_v^2} \tag{28}$$

$$
\sigma \left[ \Delta \nu \right]\_{\text{max}} = \frac{1}{c} \sqrt{\frac{1}{a^2} \sigma\_p^2 + \frac{a}{\mu} \sigma\_v^2} \tag{29}
$$

We assume that the variances of the position and velocity errors are the same standard deviation in the *x*-, *y*-, and *z*-direction in the ECI system, i.e., *σp*,*<sup>x</sup>* = *σp*,*<sup>y</sup>* = *σp*,*<sup>z</sup>* = *σ<sup>p</sup>* and *σv*,*<sup>x</sup>* = *σv*,*<sup>y</sup>* = *σv*,*<sup>z</sup>* = *σv*. Equation (29) can provide a good approximation when orbit eccentricity *e* is rather small.

It is important to point out that *σ* [*ν*]max is not equal to the exact maximum value of *σ* [*ν*] for *ν* ∈ [0, 2*π*]. This approximate maximum value is subsequently used to calculate the maximum value of QPE. With this assumption, it would be clearer to identify the core error terms of each variable. A more general expression without the assumption can be found in Appendix A.

#### *3.3. The Acceleration Vector Term in the Approximate Doppler Rate*

With the expression *σ* [Δ*ν*], the standard deviation of the difference between *ν* calculated from state vector measurements and its true value is addressed, and the acceleration vector term can now can be dealt with. As **R***st* needs to be analyzed, the traditional yaw steering method [25] is applied here.

$$\theta\_{yuv} = \arctan\left[\frac{\sin i}{N - \cos i} \cos \left(\nu + \omega\right)\right] \tag{30}$$

where *N* is the number of revolutions per day.

In this approximation model, we replace Δ**A***st* with Δ**A***<sup>s</sup>* only and ignore the errors in the target acceleration vector while keeping the satellite's part. In the meantime, *ν* is replaced by *ν* + Δ*ν* in the row vector of Equation (14). We can combine Equations (14), (18) and (20) with

$$\hat{\mathbf{u}}\_{\pi,\mathrm{E}\boldsymbol{\mathcal{I}}} = \begin{bmatrix} -\sin\nu & -\cos\nu & 0\\ \cos\nu & -\sin\nu & 0\\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \cos\theta\_{y\mathrm{uv}} & 0 & -\sin\theta\_{y\mathrm{uv}}\\ 0 & 1 & 0\\ \sin\theta\_{y\mathrm{uv}} & 0 & \cos\theta\_{y\mathrm{uv}} \end{bmatrix} \begin{bmatrix} 1 & 0 & 0\\ 0 & \cos\theta\_{\mathrm{L}} & \sin\theta\_{\mathrm{L}}\\ 0 & -\sin\theta\_{\mathrm{L}} & \cos\theta\_{\mathrm{L}} \end{bmatrix} \begin{bmatrix} 0\\ 1\\ 0 \end{bmatrix} \tag{31}$$

where *θ<sup>L</sup>* is the center off-nadir angle in stripmap SAR. The acceleration vector term results in a Doppler rate difference of

$$\begin{split} \Delta f\_{r,\mu} &= \frac{2}{\lambda \rho} \Delta \mathbf{A}\_{sl} \cdot \mathbf{R}\_{sl} \\ &= k\_{a} \left[ \cos \theta\_{L} + \sqrt{\sin^{2} \theta\_{L} \sin^{2} \theta\_{yuw} + \cos^{2} \theta\_{L}} \sin \left( \Delta \nu - \arctan \frac{\cos \theta\_{L}}{\sin \theta\_{L} \sin \theta\_{yuw}} \right) \right] \end{split} \tag{32}$$

in which

$$k\_d = \frac{2\mu \left(1 + \varepsilon \cos \nu\right)^2}{\lambda p^2} \tag{33}$$

In common SAR operations, *θyaw* is usually small. Then, we have the following expansion without cubic and higher-order terms.

$$
\arctan \frac{\cos \theta\_L}{\sin \theta\_L \sin \theta\_{yaw}} = \frac{\pi}{2} - \arctan \frac{\sin \theta\_L \sin \theta\_{yaw}}{\cos \theta\_L}
$$

$$
\approx \frac{\pi}{2} - \tan \theta\_L \sin \theta\_{yaw} \tag{34}
$$

With the result from Appendix B and Equations (32)–(34), the expected value and standard deviation of Δ*fr*,*<sup>a</sup>* are

$$\mathbb{E}\left[\Delta f\_{\mathbb{f},\mathfrak{a}}\right] = k\_{4}\cos\theta\_{L} + k\_{4}\sqrt{\sin^{2}\theta\_{L}\sin^{2}\theta\_{\text{yuw}} + \cos^{2}\theta\_{L}}\left(1 - \frac{\text{Var}\left[\Delta\nu\right]}{2}\right)\left(\frac{1}{2}\tan^{2}\theta\_{L}\sin^{2}\theta\_{\text{yuw}} - 1\right) \tag{35}$$

$$\sigma\left[\Delta f\_{7,\mu}\right] = k\_d \sqrt{\sin^2\theta\_L \sin^2\theta\_{yuw} + \cos^2\theta\_L} \sqrt{\text{Var}\left[\Delta\nu\right] \left(\tan^2\theta\_L \sin^2\theta\_{yuw} + \frac{\text{Var}\left[\Delta\nu\right]}{2}\right)}\tag{36}$$

Equation (35) shows an interesting result. Instead of the expected zero value, we have E [Δ*fr*,*a*], which means there will be bias in the expected value of the Doppler rate *fr* due to the calculated true anomaly. This bias can be and should be removed when calculating the Doppler parameters for the imaging algorithm.

However, we need the approximate maximum value of *σ* [Δ*fr*,*a*]. If we put aside the quadratic term of Var[Δ*ν*] in Equation (36), the leftover part is modulated by Var[Δ*ν*] and sin2 *θyaw* at each *ν* point. When *x* is relatively small, arctan *x* ∝ *x*. Combined with Equation (30), we can ascertain that sin2 *θyaw* ∝ sin2 [cos (*ν* + *ω*)] ∝ cos2 (*ν* + *ω*). It can also be found that, from Equation (28), if we take the appropriate approximation, then *σ* [Δ*ν*] ∝ cos *ν*. From these approximations, *ν* in *f* (*ν*) = cos (*ν* + *ω*) cos *ν* when it reaches its maximum value *f* (*ν*)*max*, and the very same *ν* will also make *σ* [Δ*fr*,*a*] reach its maximum. Obviously, at this moment,

$$
\Delta \nu = k \pi - \frac{\omega}{2}, \text{ for } k \in \mathbb{Z}\_{>0} \tag{37}
$$

If *k* = 1, then

$$\begin{split} \text{Var}\left[\left.\Delta\nu\right]\_{\text{max}}\right|\_{\nu=\pi-\omega/2} &\approx \left( \left| \cos\frac{\omega}{2} \right| \sigma \left[\Delta\nu\right]\_{\text{max}} \right)^{2} \\ &= \frac{\cos^{2}\frac{\omega}{2}}{e^{2}} \left( \frac{1}{a^{2}}\sigma\_{p}^{2} + \frac{a}{\mu}\sigma\_{v}^{2} \right) \end{split} \tag{38}$$

together with

$$k\_{a,max} = \frac{2\mu \left(1 + \varepsilon\right)^2}{\lambda p^2} \tag{39}$$

and *θyaw* = *θyaw*,*max* in Equation (36). Thus, we have

$$\begin{aligned} \sigma \left[ \Delta f\_{r,a} \right]\_{\max} &= k\_{a,\max} \sqrt{\sin^2 \theta\_L \sin^2 \theta\_{yaw,\max} + \cos^2 \theta\_L} \\ &\cdot \sqrt{\text{Var} \left[ \Delta \nu \right]\_{\max} \left( \tan^2 \theta\_L \sin^2 \theta\_{yaw,\max} + \frac{\text{Var} \left[ \Delta \nu \right]\_{\max}}{2} \right)} \end{aligned} \tag{40}$$

### *3.4. Integration Time, Slant Range, and Approximate QPE*

For a given SAR configuration and true anomaly, the integration time is used to calculate QPE, while the slant range contributes both to the integration time and QPE. In our proposed analytical approximation model for QPE, both the accurate and approximate integration time and slant range are considered. The former ones are applied to generate the approximate QPE with true anomaly as the independent variable, and the latter ones are merged into the maximum of the QPE calculation.

In the accurate calculation, the slant ranges are obtained by solving a quadratic equation:

$$\frac{\left(\rho \hat{\mathfrak{u}}\_{a,\text{ECI},x} + \mathcal{R}\_{\text{sat},\text{ECI},x}\right)^2 + \left(\rho \hat{\mathfrak{u}}\_{a,\text{ECI},y} + \mathcal{R}\_{\text{sat},\text{ECI},y}\right)^2}{E\_a^2} + \frac{\left(\rho \hat{\mathfrak{u}}\_{a,\text{ECI},z} + \mathcal{R}\_{\text{sat},\text{ECI},z}\right)^2}{E\_b^2} = 1\tag{41}$$

where *Ea* = 6,378,137.000 m and *Eb* = 6,356,752.314 m are the semi-major axis and semi-minor axis of the ellipsoid reference in WGS 84, respectively. Then, the integration time is

$$T\_{int} = \frac{\lambda}{L\_{\text{ar}}} \frac{|\mathbf{R}\_{st}|}{|\mathbf{V}\_{st}|} \frac{|\mathbf{R}\_{\text{s}}|}{|\mathbf{R}\_{t}|} \tag{42}$$

In the approximate maximum QPE calculation, we use a slant range that is based on the circular Earth with a radius of *Emean* = (2*Ea* + *Eb*) /3 and

$$\rho\_{\text{mcan}} = E\_{\text{mcan}} \frac{\sin \left[ \theta\_L + \pi - \arcsin \left( a \sin \theta\_L / E\_{\text{mcan}} \right) \right]}{\sin \theta\_L} \tag{43}$$

With the satellite's velocity at apogee only instead of |**V***st*|, together with *a* instead of |**R***s*|, we get the mean integration time.

$$T\_{int,mean} = \frac{a\lambda\rho\_{mean}}{E\_{mean}L\_a} \sqrt{\frac{a\left(1+\varepsilon\right)}{\mu\left(1-\varepsilon\right)}}\tag{44}$$

Thus, we have the approximate expected value and standard deviation of QPE and its maximum expression:

$$\mathbb{E}\left[\text{QPE}\right] = \pi \mathbb{E}\left[\Delta f\_{r,\mu}\right] \left(\frac{T\_{\text{int}}}{2}\right)^2 \tag{45}$$

$$\text{'comp.'} \qquad \sqrt{\frac{\gamma\_{\text{int},\gamma\_{\text{int},\gamma\_{\text{int},\gamma\_{\text{int},\gamma\_{\text{int},\gamma\_{\text{int},\gamma\_{\text{int},\gamma\_{\text{int},\gamma}}}}}}}{\gamma\_{\text{int}}} \left(T\_{\text{int}}\right)^2$$

$$
\sigma\left[\text{QPE}\right] = \pi\sqrt{\sigma^2\left[\Delta f\_{r,v}\right] + \sigma^2\left[\Delta f\_{r,a}\right]}\left(\frac{T\_{int}}{2}\right)^2
$$

$$
=\sqrt{\sigma^2\left[\text{QPE}\_v\right] + \sigma^2\left[\text{QPE}\_a\right]}\tag{46}
$$

$$
\sigma \left[ \text{QPE} \right]\_{\text{max}} = \pi \sqrt{\sigma^2 \left[ \Delta f\_{r,\nu} \right]\_{\text{max}} + \sigma^2 \left[ \Delta f\_{r,\mu} \right]\_{\text{max}}} \left( \frac{T\_{int,\text{mean}}}{2} \right)^2 \tag{47}
$$

where *σ* [QPE*v*] and *σ* [QPE*a*] are the standard deviations of QPE introduced by the velocity vector term and acceleration vector term, respectively.

#### **4. Evaluation and Results**

In this section, the analytical approximation model is evaluated by comparing it with the Monte Carlo simulation approach.

#### *4.1. Evaluation Model*

For the purpose of validating the proposed analytical approximation model for QPE, a Monte Carlo-based numerical simulation is introduced as a reference for validating the proposed method. In the simulation approach, large samples of error data are generated first according to the specified probability distribution parameters. Then, these error samples are input to a simulation program. For each sample of error data in the orbit determination data, the simulation program calculates the Doppler parameters of the SAR imaging procedure with and without the given error sample set and records the difference between these two situations. The probability distribution of the difference from the Doppler parameter simulation is analyzed to provide a reference for the validation of the analytical approximation model.

In this simulation program, there are no approximations such as those in the analytical approximation model. For beam steering, the simulation uses the total zero Doppler steering (TZDS) method [26] of TerraSAR-X, which has a joint yaw and pitch steering to minimize the Doppler centroid residence. The ellipsoid reference in WGS 84 is chosen to be the Earth's surface mode in the simulation.

The QPE approximation has the same slant range and integration time data as the Monte Carlo simulation in order to get a more precise approximation result for *ν* ∈ [0, 2*π*]. The QPE maximum value calculation works with the mean slant range and integration time for a much simpler computation.

The Monte Carlo method generates 30,000 error samples of the position and velocity vectors per true anomaly, and there are 1000 true anomaly calculation points in *ν* ∈ [0, 2*π*].

The inputs of the simulation are provided in Table 1. These parameters can be divided into four groups. The first group includes the Earth's physical parameters *Ea* and *Eb*, which represent the semi-major and semi-minor axes of the Earth ellipsoid, respectively. The second group includes the satellite platform orbit parameters, namely, the semi-major axes of the orbit *a*, the argument of periapsis *ω*, the orbital inclination *i*, and the eccentricity *e*. The third group contains the SAR payload parameters, i.e., the center frequency of the transmitting signal, the center off-nadir angle, and the azimuth length of the antenna. The last group represents the probability distribution of the errors in the real-time orbit determination data. The errors have a normal distribution with an expected value of zero, thus only the standard deviations are listed in the table. *σ<sup>p</sup>* and *σ<sup>v</sup>* represent the standard deviations of the errors in the position and velocity data from the real-time orbit determination system. A comparison of the 3D positioning residuals (RMS) using broadcast ephemerides results in around 1.5 m in the reference [21]. We selected a 3 m standard deviation in the real-time positioning errors. The standard deviation in the real-time velocity error is approximated by the results in [20].


#### *4.2. Approximation of the Velocity Vector Term*

Figure 2 shows the result of comparing the statistical parameters of Δ*fr*,*<sup>v</sup>* and QPE*<sup>v</sup>* between the Monte Carlo simulation and the analytical approximation model. E [QPE*v*] shows oscillations around zero, and these are due to the Monte Carlo simulation method. With more error samples processed in the numerical simulation, the oscillation would be much smaller, which confirms the effectiveness of the analytical approximation model. The mean differences in the results between the numerical simulation and analytical approximation model in Figure 2c,d are 1.64 × <sup>10</sup>−<sup>3</sup> Hz/s and 0.093◦. The difference of *σ* [Δ*fr*,*v*] may be high for the algorithm method, but it is a satisfactory result for the approximate QPE calculation. The integration time is not a consistent value, which leads to the differently shaped curves in Figure 2c,d.

**Figure 2.** Velocity vector term: numerical simulation results (blue) and analytical approximation results (red). The figures show that the results from the analytical approximation model are consistent with the numerical simulation results. (**a**,**b**) The expected values have no significant relationship with the variation in the true anomaly and are positioned at a relatively low level, i.e., around zero. The differences between the two sets of results in (**c**,**d**) are due to the ignored terms in the approximation derivations.

#### *4.3. True Anomaly Calculation*

*σ* [Δ*ν*] plays an important role in the analytical approximation model of QPE, and Figure 3 presents the results from both the numerical simulation and the analytical approximation model. Figure 3a indicates that the results of the analytical approximation model are almost identical to those of the numerical simulation. However, in the analytical approximation model, at 90◦ and 270◦, *σ* [Δ*ν*] reaches zero according to Equation (28), while the numerical simulation result is not zero at that position. In the proposed model, we assume that there are no errors in the term *p* − |**R***sat*,*ECI*| in Equation (27), which leads to a value of zero for *σ* [Δ*ν*] when the satellite is positioned at true anomalies 90◦ and 270◦. In fact, at those positions, the errors in the measurements start playing a major role in the *σ* [Δ*ν*] calculation. However, at those two specific positions, *σ* [Δ*ν*] reaches its minimum, thus it does not introduce much error to the analytical approximation model.

**Figure 3.** The true anomaly result from measurements between numerical simulation and the analytical approximation model. (**a**) The analytical approximation model (red) fits the numerical simulation (blue) in most areas of true anomaly. However, at certain true anomalies, which are the unique positions of the satellite in its orbit, the difference between the two sets of results increases. (**b**) At around 90◦ and 270◦, significant peaks arise for approximately 0.023◦ in the standard deviation of Δ*ν*.

#### *4.4. Approximation of the Acceleration Vector Term*

The expected value of the Doppler rate by the acceleration vector term E [Δ*fr*,*a*] shows different characteristics from those of the velocity vector terms. Figure 4a no longer has an expected value of zero compared with Figure 2a. This non-zero E [Δ*fr*,*a*] value is due to the calculation of true anomaly from onboard state vector measurements, which also result in a non-zero E [QPE*a*]. In addition, the analytical approximation model fits the numerical simulation result for all *ν* ∈ [0, 2*π*].

**Figure 4.** Acceleration vector term: numerical simulation results (blue) and analytical approximation results (red). The analytical approximation results are consistent with the numerical simulation results. The expected values in (**a**,**b**) no longer have an expected value of zero, indicating different characteristics from those of the velocity vector terms. The differences around 90◦ and 270◦ of true anomaly between the two sets of results in (**c**,**d**) are due to the ignored terms in the calculated true anomaly approximation derivations.

For the comparison between *σ* [Δ*fr*,*a*] and *σ* [QPE*a*], at most positions, the analytical approximation model fits the numerical simulation point by point. Apart from some simulation points, it maintains accuracy. From the comparison of Figure 4c,d, the influence of the integration time may not seem significant in contrast to the comparison of Figure 2c,d. In fact, both *σ* [QPE*a*] and *σ* [QPE*v*] are calculated with the same integration time series. This major difference in the curve appearance is caused by the relatively narrow fluctuation interval of *σ* [QPE*v*].

### *4.5. QPE Approximation*

The Monte Carlo simulation and analytical approximation model are illustrated in Figure 5. Figure 5a shows exactly the same result as that in Figure 4b, which confirms Equation (45) in the analytical approximation model. The results of *σ* [QPE] in Figure 5b, however, are presented in a slightly complex form. This is because of the combination of the properties of *σ* [QPE*v*] and *σ* [QPE*a*].

**Figure 5.** QPE result from numerical simulation (blue) compared with that from the analytical approximation model (red). The expected value and standard deviation of QPE are shown in (**a**) and (**b**), respectively. The dashed magenta line in (**b**) is the value of *σ* [QPE]*max* from Equation (47).

The dashed magenta line in Figure 5b is *σ* [QPE]*max*. This maximum value does not equal the exact maximum of *σ* [QPE]; instead, it has a 1.80% relative error with the given simulation parameters. The most important advantage of *σ* [QPE]*max* is that researchers can calculate this value easily with the aid of only a calculator and the appropriate parameters, combined with Equation (47) and related equations, and an acceptable approximation maximum value can still be determined.

#### **5. Discussion**

We designed an analytical approximation model of the QPE introduced by orbit determination errors. Our aim is to reduce the calculations and the processing time while revealing the core of the process by which error is transferred to QPE calculations. The technical parameter iteration period during early-stage development for an onboard real-time SAR imaging mission will benefit from the proposed methods.

SAR imaging requires the relative position, velocity, and acceleration vectors between the radar antenna and the targets to estimate the Doppler centroid and Doppler rate. Currently, in many ground-based processing systems for existing SAR missions, high-accuracy orbit determination data are applied to estimate the Doppler parameters. Such data have a time delay from the SAR raw data acquisition period, which is usually in the range of several hours or more. For real-time onboard SAR imaging missions, it is not realistic to wait for such high-accuracy orbit determination data to meet the "real-time" requirement. Real-time onboard orbit determination data have to be applied in the real-time onboard SAR imaging procedure. Current spaceborne SAR missions are usually based on a professional remote sensing satellite platform, and they are not as sensitive to perturbation. For a spaceborne SAR payload mounted on a small satellite platform to perform a constellation mission, such as MirrorSAR proposed by the DLR [27], there will be some differences (i.e., it is not as steady as a large remote sensing satellite platform) in the movement of the actual target satellite. In these situations, the errors in the real-time onboard orbit determination data should be evaluated in order to examine whether they will jeopardize the new real-time SAR imaging mission concept.

There exist other terms of phase error in the context of SAR data processing [28]. Besides QPE, there are linear, higher-order (more than cubic), sinusoidal, and random phase errors, among others. We selected QPE as the error of interest because it has defocus and loss-of-resolution effects on the

image. An acceptable QPE should be seen as a tool that is used out of necessity but is not adequate for new SAR missions.

Research on the prediction of SAR imaging quality from the precision of the orbit ephemeris has been carried out [29]: the upper limit expressions of multiple phase errors and their sensitivity to the orbit ephemeris covariance have been previously established. However, this method has two major limitations when applied to spaceborne real-time SAR imaging missions. Firstly, those upper limits only reflect the extreme situation, which will not always occur during a mission. We need to examine the whole situation of the phase error, as well as its probability form. Currently, the numerical simulation method can produce the probability of the phase error, but it is time-consuming. Secondly, in spaceborne real-time SAR imaging missions, it is impossible to get the orbit ephemeris directly, and it is calculated from onboard orbit determination data. The direct link between the precision of the onboard orbit determination data and the phase error in SAR imaging is still unclear. The proposed analytical approximation model is a potential solution to the two limitations above. With the precision of the onboard orbit determination system, which can be measured on the ground, the proposed analytical approximation model gives the direct QPE probability distribution result, which can simplify the analysis.

Compared with existing methods, the proposed analytical approximation model reduces the total calculations and processing time while revealing the core of the process by which errors are transferred to QPE calculations. This proposed analytical approximation model can reveal the QPE distribution, compared with the extreme value method mentioned in Section 1. For the numerical simulation method, with *M* error samples for each true anomaly calculation point and a total of *N* true anomaly calculation points, the calculation complexity reaches *O* (*MN*). It is important to point out that *M* is usually a large number for better statistical results. No error sample set is generated in the analytical approximation model; thus, the calculation complexity is reduced. The QPE distribution result for *N* true anomaly calculation points is *O* (*N*), and for the distribution parameters, the validated true anomaly *ν* ∈ [0, 2*π*] is *O* (1), which can significantly reduce the total calculation and processing time. Spaceborne real-time SAR imaging mission designers can quickly determine the maximum QPE with the help of only a calculator and the appropriate parameters. Therefore, this model will accelerate the technical parameter iteration procedure during the early-stage development of an onboard real-time SAR imaging mission and provide strong proof for the feasibility of such spaceborne real-time SAR missions.

There will always be some processing blocks in the design of a spaceborne SAR system. For each processing block, the errors introduced by processing must be limited for optimal performance of the whole system. From the design practice of the existing system, an azimuth impulse response width of 1.02 can be an acceptable limitation for a single processing block, which equals a QPE of no more than 0.27*π* [23]. The results in Section 4 are used here as an example. With the analytical approximation model, the result, and the three-sigma rule, it is clear that, with the simulation parameters given, if E [Δ*fr*,*a*] is removed during the Doppler rate calculation, then the QPE introduced by the onboard orbit determination system will be less than 40.02◦, with a probability of 99.73%. This result equals an azimuth impulse response width of less than 2%, with a probability of 99.73%.

In its current form, the proposed analytical approximation model is given when the attitude error is ignored; i.e., the orbit determination error and the attitude error of the platform are decoupled. The next generation of the analytical approximation model will take the attitude error into consideration in order to examine the coupled relationship. In addition, for spaceborne bistatic real-time SAR imaging missions, a bistatic form of the model will be needed. These two directions will be our future research topics.

### **6. Conclusions**

As the errors in the state vectors from an onboard orbit determination system for a spaceborne SAR affect real-time SAR imaging quality, this paper proposes an analytical approximation model to identify the introduction of QPE by state vector errors. This analytical approximation model can provide approximation results at two granularities, i.e., the approximation with the satellite's true anomaly as the independent variable and the approximation for *ν* ∈ [0, 2*π*]. Compared with numerical simulation methods such as Monte Carlo, the proposed analytical approximation model reduces the total calculations and processing time. Another benefit of the proposed analytical approximation model is that it reveals the core of the error transfer, so each parameter's role in the introduced QPE is clear. Spaceborne real-time SAR imaging mission designers can quickly determine QPE with the help of only a calculator and the appropriate parameters. In some ways, this will accelerate the technical parameter iteration period during early-stage development of an onboard real-time SAR imaging mission and provide convincing evidence of the feasibility of spaceborne real-time SAR missions.

**Author Contributions:** Conceptualization, X.Y. and J.C.; Funding acquisition, J.C. and O.L.; Investigation, X.Y. and H.N.; Methodology, X.Y.; Project administration, J.C. and O.L.; Resources, X.Y., J.C., H.N., and O.L.; Supervision, J.C. and O.L.; Validation, X.Y. and H.N.; Visualization, X.Y. and H.N.; Writing—original draft, X.Y.; and Writing—review & editing, X.Y. and H.N.

**Funding:** This research was funded by the National Natural Science Foundation of China under Grant 61861136008 and Grant 61671043. The author Xiaoyu Yan was funded by China Scholarship Council (CSC) under Grant 201706020033.

**Acknowledgments:** The authors thank Yue Fang from Beihang University for useful discussions. The authors are grateful to Xi (Cathy) Chen from Northwestern University for her understanding and useful suggestions.

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

#### **Appendix A. Errors in Calculating True Anomaly**

Let us have

$$\mathbf{Y} = \sqrt{\frac{p}{\mu}} \mathbf{V}\_{\text{sat}, ECI} \cdot \mathbf{R}\_{\text{sat}, ECI} \tag{A1}$$

$$X = p - |\mathbf{R}\_{sat,ECI}|\tag{A2}$$

with the assumption that the position error has no influence on term *X* in Equation (A2). Thus, we have

$$
\Delta \frac{\chi}{X} = \sqrt{\frac{p}{\mu}} \cdot \frac{1}{p-r} \Delta \mathbf{V} \cdot \mathbf{R} \tag{A3}
$$

The term Δ**V** · **R** is the difference in the term **V***sat*,*ECI* · **R***sat*,*ECI* in Equation (A1) between the measured data and the accurate value.

$$\begin{split} \Delta \mathbf{V} \cdot \mathbf{R} &= \mathbf{V}\_{\mathrm{sat}, \mathrm{ECI}, \ell} \cdot \mathbf{R}\_{\mathrm{sat}, \mathrm{ECI}, \ell} - \mathbf{V}\_{\mathrm{sat}, \mathrm{ECI} \cdot \mathbf{R}} \cdot \mathbf{R}\_{\mathrm{sat}, \mathrm{ECI}} \\ &= (\mathbf{A}\_{\mathrm{O2E}} \mathbf{V}\_{\mathrm{sat}, \mathrm{OP}} + \Delta \mathbf{V}) \cdot (\mathbf{A}\_{\mathrm{O2E}} \mathbf{P}\_{\mathrm{sat}, \mathrm{OP}} + \Delta \mathbf{P}) - \mathbf{A}\_{\mathrm{O2E}} \mathbf{V}\_{\mathrm{sat}, \mathrm{OP}} \cdot \mathbf{A}\_{\mathrm{O2E}} \mathbf{P}\_{\mathrm{sat}, \mathrm{OP}} \\ &\approx c\_{\mathrm{P}} \sqrt{\frac{\mu}{p}} + c\_{\mathrm{U}} \frac{a \left(1 - \varepsilon^{2}\right)}{1 + \varepsilon \cos \nu} \end{split} \tag{A4}$$

where the position error times the velocity error are ignored and

$$\mathcal{L}\_p = -\Delta p\_x \left[ \sin \left( \nu + \omega \right) + \varepsilon \sin \omega \right] + \Delta p\_y \left[ \cos \left( \nu + \omega \right) + \varepsilon \cos \omega \right] \cos i + \Delta p\_z \left[ \cos \left( \nu + \omega \right) + \varepsilon \cos \omega \right] \sin i \tag{A5}$$

$$c\_{\upsilon} = \Delta v\_x \cos \left(\nu + \omega \right) + \Delta v\_y \sin \left(\nu + \omega \right) \cos i + \Delta v\_z \sin \left(\nu + \omega \right) \sin i \tag{A6}$$

The *cp* and *cv* parameters in Equations (A5) and (A6) can be regarded as the core function of the error introduced by position and velocity errors when calculating true anomaly. The calculated true anomaly's variance depends directly on the variance of the terms *cp* and *cv*, which link directly to the variance of the measurement errors of position and velocity in each axis of the ECI system.

Considering the derivative of Equation (27), with *Y*/*X* = tan *ν*, we have the difference in calculated *ν* and its real value.

$$\begin{split} \Lambda \nu &= \left( \operatorname{atan2} \frac{Y}{X} \right)' \cdot \Lambda \frac{Y}{X} \\ &= \frac{1}{1 + \left( Y/X \right)^{2}} \cdot \sqrt{\frac{p}{\mu}} \cdot \frac{1}{p - r} \cdot \Delta \mathbf{V} \cdot \mathbf{R} \\ &= \frac{\cos \nu \left( 1 + \epsilon \cos \nu \right)}{\varepsilon \sqrt{\mu a \left( 1 - \epsilon^{2} \right)}} \left[ c\_{p} \sqrt{\frac{\mu}{p}} + c\_{\upsilon} \frac{a \left( 1 - \epsilon^{2} \right)}{1 + \epsilon \cos \nu} \right] \end{split} \tag{A7}$$

The variance of Δ*ν* is

$$\text{Var}\left[\Delta\nu\right] = \cos^2\nu \cdot \frac{\left(1 + \varepsilon\cos\nu\right)^2}{\mu a e^2 \left(1 - e^2\right)} \cdot \left[\frac{\mu}{p}\text{Var}\left[c\_P\right] + \left(\frac{a\left(1 - e^2\right)}{1 + \varepsilon\cos\nu}\right)^2 \text{Var}\left[c\_U\right]\right] \tag{A8}$$

A general expression of the variance of *cp* and *cv* is

$$\begin{split} \text{Var}\left[\mathfrak{c}\_{p}\right] = \text{Var}\left[\Delta p\_{x}\right]\left[\sin\left(\nu+\omega\right)+\varepsilon\sin\omega\right]^{2}+\text{Var}\left[\Delta p\_{y}\right]\left[\cos\left(\nu+\omega\right)+\varepsilon\cos\omega\right]^{2}\cos^{2}i \\ &+ \text{Var}\left[\Delta p\_{z}\right]\left[\cos\left(\nu+\omega\right)+\varepsilon\cos\omega\right]^{2}\sin^{2}i \end{split} \tag{A9}$$

$$\text{Var}\left[\varepsilon\_{\nu}\right] = \text{Var}\left[\Delta v\_{x}\right]\cos^{2}\left(\nu + \omega\right) + \text{Var}\left[\Delta v\_{y}\right]\sin^{2}\left(\nu + \omega\right)\cos^{2}i + \text{Var}\left[\Delta v\_{z}\right]\sin^{2}\left(\nu + \omega\right)\sin^{2}i \tag{A10}$$

If all the variances of position and velocity errors are the same standard deviation in the *x*-, *y*-, and *z*-directions in the ECI system, i.e., *σp*,*<sup>x</sup>* = *σp*,*<sup>y</sup>* = *σp*,*<sup>z</sup>* = *σ<sup>p</sup>* and *σv*,*<sup>x</sup>* = *σv*,*<sup>y</sup>* = *σv*,*<sup>z</sup>* = *σv*, simpler expressions for *cp* and *cv* can be given by

$$\text{Var}\left[\mathfrak{c}\_{p}\right] = \sigma\_{p}^{2}\left(1 + e^{2} + 2\varepsilon\cos\nu\right) \tag{A11}$$

$$\text{Var}\left[c\_v\right] = \sigma\_v^2 \tag{A12}$$

With Equations (A8), (A11) and (A12), we can derive Equation (28) in Section 3.2. This expression can be used to determine the dependent variable *σ* [Δ*ν*] with the independent variable *ν*. In addition, we need a simpler expression of the maximum *σ* [Δ*ν*] in order to support the quick calculation method for QPE. Beginning with Equation (28), we can ignore some terms in the equation related to orbit eccentricity *e* when its value is rather small (less than 0.01). Thus, we have

$$
\sigma \left[ \Delta \nu \right]\_{\text{max}} = \frac{1}{\varepsilon} \sqrt{\frac{\sigma\_p^2}{a^2} + \frac{a \sigma\_v^2}{\mu}} \tag{A13}
$$

which is Equation (29) in Section 3.2.
