Integrity monitoring for precise orbit determination of LEO satellites

Due to an increasing requirement for high accuracy orbital information for low Earth orbit (LEO) satellites, precise orbit determination (POD) of LEO satellites is a topic of growing interest. To assure the safety and reliability of the applications requiring high accuracy LEO orbits in near-real-time, integrity monitoring (IM) is an essential operation of the POD process. In this contribution, the IM strategy for LEO POD in both the kinematic and reduced-dynamic modes is investigated. The overbounding parameters of the signal-in-space range error are investigated for the GPS products provided by the International GNSS Service’s Real-Time Service and the Multi-GNSS Advanced Demonstration of Orbit and Clock Analysis service. Benefiting from the dynamic models used and the improved model strength, the test results based on the data of the LEO satellite GRACE FO-1 show that the average-case mean protection levels (PLs) can be reduced from about 3–4 m in the kinematic mode to about 1 m in the reduced-dynamic mode in the radial, along-track and cross-track directions. The overbounding mean values of the SISRE play the dominant role in the final PLs. In the reduced-dynamic mode and average-case projection, the IM availabilities reach above 99% in the radial, along-track and cross-track directions with the alert limit (AL) set to 2 m. The values are still above 98% with the AL set to 4 m, when the duty cycle of tracking is reduced to 40%, e.g., in the case of power limits for miniature satellites such as CubeSats.


Introduction
Precise orbit determination (POD) of low Earth orbit (LEO) satellites has been intensively investigated in the past decades (Montenbruck et al. 2005;Gu et al. 2017). By applying the reduced-dynamic model integrated with GNSS observations tracked onboard and the use of comprehensive dynamic models, the accuracy of LEO POD can reach a few centimeters in "near-real-time," as a result of the availability of precise GNSS orbital and clock products provided by realtime services (Montenbruck et al. 2013;Allahvirdi-Zadeh et al. 2021). Such (near)-real-time LEO POD is required for applications such as Interferometric Synthetic Aperture Radar (InSAR) (Montenbruck 2017) and radio occultation analyses (Montenbruck et al. 2005). In recent years, LEO mega-constellations have been proposed as an augmentation to GNSS-based positioning, navigation and timing (PNT) services (Reid et al. 2018;Han et al. 2020). Compared to traditional GNSS-based positioning, LEO-aided GNSS positioning can benefit from significantly increased satellite numbers, improved satellite geometry, stronger signal strengths (Reid et al. 2018) and greater satellite velocities, which reduce the precise point positioning (PPP) convergence time (Ge et al. 2018). As one of the conditions for LEO-aided positioning, the precise orbits of the LEO satellites are essential information.
Integrity monitoring (IM) in the near-real-time of the LEO satellite POD is important to ensure reliable positioning, thus ensuring safety and reliability for diverse user applications. With a pre-defined probability of misleading information (PMI), the protection level (PL) is a safety parameter that is computed to bound the orbital positioning errors (PEs). When the PL exceeds a pre-defined alert limit (AL), a warning message can be sent to the LEO system so that the affected satellite can be, e.g., marked as "unhealthy" and excluded from further processing. While LEO POD has been intensively studied, the IM of LEO POD is a new topic worthy of study. IM has been investigated for different positioning methods for airborne and ground-based users in the last few decades. In aviation, the Advanced Receiver Autonomous Integrity Monitoring (ARAIM) algorithm (Blanch et al. 2012) was developed for dual-frequency multi-constellation (DFMC) signals to facilitate the vertical and horizontal guidance of aircraft. The algorithm was extended to ground-based applications using other positioning methods (El-Mowafy and Kubo 2018; Hassan et al. 2020). At the same time, the traditional single-frequency single-constellation satellite-based augmentation system (SBAS) is transitioning to a new generation that supports DFMC users (Wu et al. 2020). Its IM method is defined for aeronautical users (EUROCAE 2019), and extended for ground-based positioning of users for the more complex measurement environments encountered on the ground .
As a user of the GNSS signals in space, the IM of the dual-frequency LEO POD shares similar principles with the aeronautical users using dual-frequency GNSS signals. For instance, LEO satellites are flying under "open-sky" conditions and experience a relatively simple multipath environment. Furthermore, the first-order ionospheric delays are eliminated through the ionosphere-free (IF) combination. However, situations are also different compared to, e.g., those in aviation using the ARAIM or SBAS methods. The tropospheric residuals (after modeling) do not need to be considered in the IM anymore due to the altitude of the LEO satellites. Instead of the carrier-smoothed code observations used in the ARAIM and SBAS, the phase observations are directly used in LEO POD and hence influence the final POD precision and accuracy.
Near-real-time LEO POD relies on high accuracy realtime GNSS orbital and clock products generated by different providers. Service such as the International GNSS Service (IGS)'s Real-Time Service (RTS) (https:// www. igs. org/ rts/ produ cts/; El-Mowafy et al. 2017) is one example of ground-based processing that provides real-time POD. Services with corrections broadcast through GNSS or GEO satellites, e.g., the Multi-GNSS Advanced Demonstration of Orbit and Clock Analysis (MADOCA) service (Zhang et al. 2019) based on the Japanese Quasi-Zenith Satellite System (QZSS), and commercial services provided by, e.g., Fugro (Hauschild et al. 2016), also demonstrate possibilities for future onboard high accuracy POD. As the signalin-space range errors (SISRE) of these real-time products do not strictly follow a zero-mean normal distribution (Rife et al. 2006) as assumed in the POD least-squares processing, the overbounding parameters of the SISRE need to be determined using appropriate methods. The impact of different overbounding parameters on the LEO orbital PLs also requires investigation. Finally, LEO POD in the reduceddynamic mode involves the estimation of orbital dynamic parameters and the numerical integration of the orbits. As a result, numerically integrated partial derivatives of the Cartesian orbits with respect to the estimable orbital parameters are needed in the IM procedure for reduced-dynamic POD.
In this study, an IM procedure is proposed for LEO POD. The overbounding parameters are then investigated for the SISRE of real-time GNSS products under the average-case and worst-case for the orbital and clock error projection in the signal direction. Subsequently, the POD results and the PLs are assessed using real data from the GRACE Follow-On mission. The impacts of the overbounding parameters of the SISRE on the PLs are studied, and the PLs are investigated for miniature satellites for scenarios of non-continuous satellite tracking. The conclusions are then drawn and presented at the end.

Processing strategy
LEO POD and IM are performed for the reduced-dynamic mode, combining dynamic models with the GPS L1 and L2 phase observations. The code observations are used to pre-process the a priori clock synchronization (Dach et al. 2015). The IF combination is formed to remove the firstorder ionospheric delays so that only the orbital parameters, the LEO satellite clocks Δt L , and the IF float-valued phase ambiguities n IF remain as estimable parameters in the linearized observation equations: With where E(•) is the expectation operator, and Δ IF denotes the observed-minus-computed (O-C) terms of the IF phase observations. c denotes the speed of light, and f j , n j , j ( j = 1, 2 ) denotes the frequency, ambiguity and wavelength on frequency j , respectively. As the absolute values of LEO satellite clocks and the IF float ambiguities are not of concern in this study, a reference ambiguity n 0 IF (with the value of zero) was set to avoid the rank deficiency between the LEO satellite clocks and the ambiguities. The terms G IF and d G IF denote the GPS IF satellite phase and code bias vector, respectively, and L,IF refers to the LEO IF satellite phase bias. ΔK and ΔD are the increment vectors of the six Keplerian elements at the initial state and the time-constant orbital dynamic parameters with respect to the a priori values, respectively. denotes the vector of the stochastic piece-wise constant accelerations. ΔD and are given in the radial (R), along-track (S) and cross-track (W) directions with the stochastic accelerations estimated, in this study, within 6 min time intervals and constrained to zeros with a standard deviation ( ) of 5 × 10 −6 [s −2 ] × 0 . The a priori standard deviation of unit weight ( 0 ) refers to the zenith direction on L1. A K , A D and A contain the partial derivatives of the observations with respect to the Keplerian elements, the dynamic parameters and the stochastic accelerations, respectively, which can be further expressed as: where IF ∕ r I denotes the partial derivatives of the observation with respect to the LEO position vector in the inertial system, i.e., the GPS-to-LEO unit direction vectors. r I ∕ K and r I / D are the partial derivatives of the LEO position vector with respect to the Keplerian elements and the dynamic parameters, respectively. These partial derivatives are numerically integrated based on the variational equations (Montenbruck and Gill 2000). r I / , the partial derivatives of the position vector with respect to the stochastic accelerations, can be computed as linear combinations of r I ∕ K (Dach et al. 2015).
The unknowns, contained in the increment vector Δx containing all the orbital dynamic parameters, the LEO satellite clock parameters of all epochs and the ambiguity parameters of all epochs and all satellites, are estimated in a least-squares adjustment. The variance-covariance Q l of all the observations can be formulated as a block diagonal matrix constructed from the variance-covariance matrices of the phase O-C terms at each time point t i , denoted as Q Δ t i , and the variance matrix of the constraints for stochastic accelerations: where the operator blkdiag (•) forms the block diagonal matrix. For epoch t i , the variance-covariance matrix of the IF O-C terms is expressed as a diagonal matrix of the corresponding variances for each satellite s (1, ⋯ ,m ), with m denoting the number of satellites: In this study, as was done for Sentinel satellite POD products and in the Bernese LEO campaign, equal weighting is applied for all the phase IF O-C terms with 2 s = (f 4 1 + f 4 2 )∕ f 2 1 − f 2 2 2 × 2 0 , as the GNSS receiver onboard LEO satellite experiences less elevation-dependency for the noise and multipath effects than on the ground, so that traditional elevation-dependent weighting models might not lead to better results.
For LEO POD, the least-squares adjustment is performed using 24 h of GPS observations with a sampling interval of 30 s, which is sufficient for LEO POD ). The dynamic models used for numerical integration are summarized in Table 1. The LEO orbits are numerically integrated in the inertial system with all the orbital parameters solved in a least-squares adjustment.
In addition to the reduced-dynamic mode, for comparison purposes, the POD is also performed in the kinematic mode (i.e., without using any dynamic model). In such a case, the orbits are directly estimated in the form of the Cartesian coordinates of the satellite in the Earth-fixed system.

Integrity monitoring
For the IM, the variance-covariance matrix of the observations is constructed assuming Gaussian noise that overbounds all existing errors remaining in the O-C terms. The variance-covariance Q l of the observations for the purpose of integrity can be formulated as:  (Lyard et al. 2006) Taking into account the factors considered for IM in aviation using the DFMC SBAS (EUROCAE 2019) that also operates in an open-sky measurement environment and uses undifferenced IF observations, the variance element ̃2 s ( s = 1, ⋯ , m ) can be defined considering the variance of the noise and multipath ( ̃2 air ), the variance of the GNSS orbital and clock errors projected into the signal direction ( ̃2 GNSS ) and the variance of the ionospheric residuals after forming the IF cothe overbounding standardmbination ( ̃2 iono ), which is assumed to account for miss-modeled ionospheric residuals that will be discussed later: Note that LEO satellites are above the troposphere, and thus, there are no tropospheric delays to consider. The variance where ̃N ,DFMC and ̃M ,DFMC are the overbounding standard deviations of the IF carrier-smoothed observation noise and multipath after 100 epochs using L1 and L5, applied for aeronautical users with the DFMC SBAS as defined in EUROCAE (2019). The factor (13) is set to 10, assuming the standard deviation of the raw code noise and multipath is 100 times that of the phase noise and multipath, and the code-smoothing over 100 epochs reduces this amplification factor to 10. The calculation of ̃i ono follows that for the case of DFMC SBAS, i.e., an elevation-angle dependent function accounting for the high-order ionospheric delays, ray bending and excess Total Electron Content (TEC). The ionosphere differently influences the LEO satellites at different altitudes. Some of them could even be above the ionosphere; therefore, we made a conservative assumption using the same formula as for the aircraft. The overbounding standard deviation for the orbital and clock errors projected into the signal direction ( ̃G NSS ) is dependent on the service provider of the real-time GNSS products.
Recall that stochastic accelerations are constrained to zero with the standard deviation = 5 × 10 −6 [s −2 ] × 0 . To keep the same weight of constraint for the stochastic accelerations as that in the POD, the standard deviation for the constraint ( ̃2 in m∕s 2 ) for integrity purposes is expressed as: where ̃0 denotes the overbounding standard deviation of unit weight on L1 (in the unit of meters), calculated with the in the zenith direction divided by the IF amplification factor (10), the variance-covariance matrix Q Δx of the unknown vector is: where A is the total design matrix (including the part of the absolute constraint) for the reduced-dynamic POD. Next, the variance-covariance matrix of the orbital parameters denoted as Q O , is derived from Q Δx as the first 729 × 729 elements in the reduced-dynamic mode, accounting for six Keplerian elements at the initial condition, one set of 3D time-constant dynamic parameters, and 240 sets of 3D piece-wise constant stochastic accelerations (within 6 min intervals). The variance-covariance matrix of the orbital parameters Q O obtained through the least-squares adjustment is transformed to that of the Cartesian coordinates in the radial, along-track and cross-track directions at each epoch using the error propagation law, such that: with the transformation matrix F t i containing the partial derivatives of the RSW coordinates with respect to the orbital parameters expressed as: where r I t i ∕ O represents the partial derivatives of the position vector in the inertial system at t i with respect to all the orbital parameters O , which need to be computed by numerical integration based on the variational equations and their linear combinations (as mentioned earlier). R I2T t i denotes the rotation matrix from the inertial system to the RSW-system at t i . In this study, R I2T at t i is approximated by: where r I and v I represents the LEO position and velocity vectors, respectively, and ‖ • ‖ is the norm of the vector. For IM, the PL is computed covering the positional errors with a pre-defined PMI (Shively 2014) as mentioned in the Introduction. Assuming that the PMI is equally distributed in the three directions, the PLs of the orbital element q can then be formulated as: is the inverse CDF of the standard normal distribution, ̃R SW,q t i is the overbounding standard deviation of the q-th orbital element at epoch t i , with q=1, 2, 3 representing the orbital component in the radial, along-track and cross-track directions, respectively, which can be derived from Q RSW in Eq. 16. m GNSS denotes here the scalar overbounding mean value of the GNSS orbital and clock error projections that are computed along with the ̃G NSS (Eq. 12) (which will be explained in the next section). J l Δ is a vector of ones with length l Δ , which denotes the number of used phase IF observations. S RSW,q projects the bias vector b from the observation to the position domain at the corresponding epoch. The matrix F O selects the orbital parameters from the entire unknown vector containing the orbital parameters, the LEO satellite clock parameters and the float ambiguity parameters (see Eq. 1), and e q is a 3 × 1 vector with the q-th element set to the value of one and the other elements set to zeros.
In the proposed algorithm, large cycle slips or biases will be detected by the fault detection and exclusion (FDE) step. Concerning the remaining smaller biases that could vary over time, compared to the ambiguity-fixed solutions, the float ambiguities could absorb satellite-specific errors that vary little with time. The values of the float ambiguities themselves are not of concern. Other small miss-detected biases errors could remain in the observation residuals considered in b (Eq. 19), which represents possible nominal biases under the system normal operation, and their random variabilities could be considered in the stochastic term ̃2 s (Eq. 12). These parameters can be properly set based on a future comprehensive study.
Compared to the reduced-dynamic mode, the procedure for the IM in the kinematic mode is simpler, as the variance-covariance matrix of the Cartesian orbits can be directly obtained within the least-squares adjustment as: where the design matrix A k can be expressed with the help of the combined GNSS-to-LEO unit direction vectors for all observations ( A I,n ), the design matrix for the LEO clock parameters A T,n , as well as that for the ambiguity vector ( Λ IF ) as: The variance-covariance matrix of the orbits Q r I t i ( i = 1, ⋯ , n ) are contained in the first 3n × 3n elements of Q Δx k . To obtain the variance-covariance matrix of the orbits in the radial, along-track and cross-track directions, i.e., Q RSW,k t i , similar to the reduced-dynamic mode, transformation is required from the inertial system to the RSW-system: with the overbounding standard deviation of the orbits obtained in the radial, along-track and cross-track directions ( ̃R SW,k,q ), the corresponding protection levels in the kinematic mode ( PL k,q ) can be computed: with where F k selects the kinematic orbital positions of the corresponding epoch from the entire vector of unknowns. In summary, compared to the existing IM algorithm for the DFMC SBAS, the LEO IM are subject to the following changes: -The reduced-dynamic LEO POD requires a more complicated scheme for the projection from the observation domain to the position domain as, instead of directly computing the Cartesian coordinates, the dynamic orbital parameters are estimated. In this regard, the objective of the least-squares adjustment in the reduced-dynamic POD is different from that of DFMC SBAS positioning.
In the reduced-dynamic POD, the GNSS observations are used to improve the dynamic parameters in a leastsquares adjustment so that the final orbits are obtained through numerical integration using an improved dynamic model. -Concerning the real-time GNSS satellite orbits and clocks, the LEO POD has more options, which may exhibit higher accuracy than using SBAS corrections. -The model for LEO POD is stronger than that used by ground-based or airborne users, as LEO satellites are above the troposphere (i.e., it uses less term when com-

SISRE analysis of real-time GNSS services
As mentioned in the previous section, the overbounding standard deviation ( ̃G NSS ) and mean values ( m GNSS ) of the GNSS orbital and clock errors, projected into the signal direction, are important components in the computation of PLs. In this study, two different cases are considered: the average-case and the worst-case.

Average-case projection
In the average-case, the SISRE assumes a uniform distribution of the points over the sphere of the approximate LEO orbital radius, with a projection factor R defined for the radial direction, and a factor SW for the along-track and cross-track directions, respectively (Heng et al. 2011). The overbounding standard deviation can thus be expressed as: where ̃R Δt G is the overbounding standard deviation of the term R × Δr R − c × Δt G considering the correlation between the radial orbital error Δr R and the satellite clock error Δt G of the GNSS satellites. ̃S and ̃W represent the overbounding standard deviations of the along-track and cross-track orbital errors of the GNSS satellites, respectively. The overbounding mean value of the projected GNSS orbital and clock errors can be expressed as: where m RΔt G is the overbounding mean value of the term R × Δr R − c × Δt G , and m S and m W refer to those of the along-track and cross-track orbital errors of the GNSS satellites, respectively. The factors R and SW are related to both the orbital heights of the LEO satellite and the GNSS satellites and do not vary much for LEO orbital altitudes from 500 to 1000 km . The radial orbital error and the clock error play the dominant role in the overbounding standard deviation and mean values. In this study, the overbounding standard deviation ̃G NSS and mean values m GNSS are computed for two types of realtime GPS products, i.e., those from the IGS RTS and the MADOCA Service. The IGS final orbits and clocks are used as the reference for the assessment of the errors. Note that the real-time GPS satellite clocks are re-referenced to those of the IGS final clocks. The satellite-specific timeconstant mean values are removed from Δt G , as they can be absorbed by the float ambiguities (Zhang et al. 2019). The errors (all the single components and the radial-clock combined component) exceeding 0.5 m are excluded from the analysis, and the remaining errors exceeding m err ± 10 × err are further screened as outliers, where m err and err denote the mean value and standard deviation of the corresponding error type, respectively. These outliers are also screened in the phase-preprocessing during the POD process.
A two-step overbounding strategy is used for the computation of the overbounding standard deviations and mean values (Blanch et al. 2018). Compared to the paired overbounding method (Rife et al. 2006), the two-step strategy has the advantage of generating a smaller overbounding mean value, which is essential for reducing the LEO PLs, and thus may improve solution availability, as will be explained in Section "Test results." In this study, no excess mass is set for the overbounding to avoid the high inflation of the actual probability of the overbounding distribution due to a large number of observations in the batch leastsquares POD process.  Table 2 lists the contributions of different errors to the overbounding standard deviation ̃G NSS and mean value m GNSS for the orbital and clock projection errors assuming the average-case over the period September 2018. The LEO GRACE FO-1 satellite, with approximately 503 km orbital height (Wen et al. 2019), for the period from August 14 to 20, 2018 was used as an example. To show the correlation between the terms ( R × Δr R ) and (c × Δt G ) , the contributions are also given for the corresponding single terms with ̃R and m R denoting the overbounding standard deviation and mean value of the radial orbital error, and ̃Δ t G and m Δt G denoting those of the clock error. For the purpose of comparison, all epochs having valid components are used for the analysis.
From Table 2 it can be seen that ̃R Δt G is smaller than the case when treating the radial orbital errors and the clocks as independent terms, which confirms the correlation between the two parameters. The situation is similar for the overbounding mean values m RΔt G . The m GNSS of both products are within 5 cm, and their ̃G NSS are at the sub-dm to dm-level.

Worst-case projection
The worst-case projection assumes the largest projection of the orbital and clock errors. Based on Montenbruck et al. (2015), the worst-case for the orbital contribution in SISRE is expanded from the case for ground-based users to the case for LEO satellites as shown in Fig. 1.
As the LEO orbits are near-circular, they are assumed to be circles with a radius R L . When the orbital error vector Δr intersects with the sphere of radius R L (see Δr in in the right panel of Fig. 1), the worst-case orbital error Δr is assumed to have the absolute value of ‖Δr‖: With where < • > is the dot product of two vectors, and R L and R G are approximated by the mean LEO orbital height and the mean GPS orbital height of 20,200 km (with the addition of the Earth radius of 6371 km as used in the previous subsection). e G denotes the unit direction vector from the Earth to the GPS satellite, and Δr R refers to the radial component of Δr . In the case when the orbital error vector does not intersect the sphere of radius R L (see Δr out in the right panel of Fig. 1), the orbital error vector Δr is projected into the nearest tangent line from the GNSS satellite to the sphere for LEO orbits using the relations: As the worst-case orbital contributions to SISRE might not be normally distributed, similar to the case in the previous section, an overbounding standard deviation ( ̃Δ̃r ) and mean value ( m Δr ) need to be determined for Δr . As a pessimistic (31) Δr = ‖Δr‖, for the yellow area with ≤ (32) Δr = ‖Δr‖, for the blue area with ≤ − Δr = ‖Δr‖ × cos ( − ), for the grey area with, < ≤ 2 (36) Δr = ‖Δr‖ × cos ( + ), for the orange area with, 2 < ≤ − Fig. 1 Average-case (left) and the worst-case (right) of the orbital contribution in SISRE for LEO satellites assumption, without considering the correlation between the radial orbits and the satellite clock errors (Montenbruck et al. 2015), the overbounding standard deviation and mean value of the projected GNSS orbital and clock errors can be expressed as: Using the same data set as before, the overbounding standard deviations and mean values in the single and combined terms are given in Table 3. For the purpose of comparison, ̃Δ̃r Δt G and m ΔrΔt G are given for the case when the correlations between Δr and Δt G are considered, i.e., for the term ( Δr − c × Δt G ). Outlier rejection was performed as in the previous section for the single and combined items of Δr and Δt G . Epochs having valid values for all terms, i.e., the clocks, the radial orbits, and the combined clock and radial orbits, are used for the analysis. From Table 3, it can be seen that ̃Δ̃r and m Δr in the worstcase are higher than the corresponding orbital contributions under the average-case (as listed in Table 2), i.e., the terms √ 2 R ×̃2 R + 2 SW × ̃2 S +̃2 W a n d R ×m R + SW × m S +m W . Likewise, the overbounding standard deviations and mean values of the combined orbital and clock errors ̃Δ̃r Δt G and m ΔrΔt G are lower than those computed based on the single terms, i.e., the final ̃G NSS and m GNSS . Note that the slight differences in ̃Δ t G and m Δt G compared to those in the average-case are caused by the different epochs used in the analysis due to the outlier rejection of different parameters. Overall, the ̃G NSS and m GNSS are higher than those in the average-case, i.e., around 0.1 m for ̃G NSS and sub-dm to dmlevel for m GNSS . The results are largely consistent with the User Range Accuracy (URA) analysis performed by Cheng et al. (2018) for the CNES (National Center for Space Studies) realtime products, stating that a URA of 0.15 m is a reasonable value for the GPS satellites. Note that the use of an overbounding mean value helps in reducing the overbounding standard deviation.

Test results
In this section, the orbital PEs and the PLs are processed for the LEO satellite GRACE FO-1 during the period August both illustrated for the PLs in Fig. 2 as red and blue dots, respectively. Note that the bias contribution needs to be considered in IM, as without the mean value m GNSS , a much larger overbounding standard deviation ̃G NSS could be needed to satisfy the overbounding characteristic of a Gaussian distribution. The blue dots in Fig. 2 are shown to illustrate the contribution of the bias term.
From Fig. 2, it can be seen that the PLs overbound the PEs in all three directions, and in both modes, for the test day. Due to the increased model strength, although having the same overbounding parameters, the PLs are much lower in the reduced-dynamic mode than in the kinematic mode. In addition, the highest PLs in the radial direction in the kinematic mode, caused by the high correlation between the radial component and the satellite clocks, do not appear to be higher than those in the other directions in the reduceddynamic mode due to de-correlation between the radial orbital component and the satellite clock parameter after applying the orbital dynamic models. Comparing Eq. (19) and Eq. (26), both the formal standard deviation ̃R SW,q and projection matrix S RSW,q are increased in the kinematic mode due to the weaker model strength without using the dynamic model. More concretely, compared to the kinematic POD that solely relies on the GNSS observations, the reduced-dynamic model gains more information from both the GNSS observations and the dynamic models.
Comparing the red and the blue dots in Fig. 2 i.e., with and without considering biases in the PLs, it can be seen that the bias contribution plays a significant role in the final PLs. This is partially caused by the fact that the absolute values of the projection matrices S RSW,q t i and S RSW,k,q t i (see Section "Integrity Monitoring") were used in the IM, which pessimistically (as a conservative measure) assumes that the bias term always has the same sign as the corresponding values in the projection matrix. Figure 3 shows the absolute PEs and the PLs for the average-condition applying the MADOCA overbounding parameters for all three directions. The LEO data of GRACE FO-1 from August 14 to 19, 2018, were processed with the ̃G NSS and the m GNSS taken from Tables 2 and 3 for the averagecase and the MADOCA products. It can be observed that the reduced-dynamic mode has significantly reduced both of the PEs and the PLs compared to the kinematic mode. The figure shows that the PLs bound all PEs in the three directions and during the entire test period (see the 0.2 m cap in the right panel). Table 4 lists the mean value of the absolute PEs and the PLs, and the corresponding IM availabilities with an assumed AL of 4 m. This is an example presented for demonstration purposes, showing the availability differences for different estimation modes and under the average-case and the worst-case projections. The availability is defined as the ratio when PL is smaller than AL in the corresponding direction. The values in Table 4 distinguish between the average-case and the worst-case.
While the PLs in the average-case projection can indeed bound all the PEs in the three directions and during the entire test period, the worst-case projection provides safer but overall pessimistic PLs. Benefiting from the increased model strength, as shown in Table 4, applying the dynamic model can effectively reduce the mean PLs from about and PLs assuming the average-condition applying the overbounding parameters of the MADOCA products (right). Note that the y-scales of the two sub-figures are different 3-4 m in the kinematic mode to around 1 m in the reduceddynamic mode for all three directions for the average-case assumption. Under such conditions, the availabilities reach 100% for an AL of 4 m as shown in Table 4, and are above 99% even with the AL decreased to 2 m.

PLs with different overbounding parameters
To test the impact of different overbounding parameters, the PLs were computed using the satellite geometry provided by the GRACE FO-1 and the MADOCA products on August 14, 2018, with the value of ̃G NSS varying from 0.05 to 0.15 m and the value of m GNSS varying from 0 to 0.15 m. Figure 4 shows the change of the mean PLs with respect to the m GNSS in the kinematic (top) and reduced-dynamic modes (bottom).
As shown in Fig. 4, compared to the ̃G NSS , the mean PLs are more influenced by the m GNSS of the projected GNSS orbital and clock errors. The results for LEO POD motivate investigations into further reducing the overbounding mean values. With m GNSS within 5 cm, the mean values of PLs in the reduced-dynamic mode are within 2 m for the averagecase even with a large ̃G NSS up to 0.15 m. With m GNSS within 2 cm, the mean values of PLs are limited to 1 m under the same conditions.

PLs with different duty cycles
For mini-satellites such as CubeSats, the limited onboard power might not allow for GNSS observations to be continuously collected over an entire period, as power might need to be used for tasks other than positioning. As such, duty cycles smaller than 100% are often applied for the GNSS tracking (Lantto and Gross 2018). The reduced data decreases the model strength and the precision of the estimated orbital parameters. However, in the reduced-dynamic mode the orbital positions and PLs can still be computed over the entire period, benefiting from the dynamic model and numerical integration .
The mean PLs are computed and shown in Fig. 5 for duty cycles varying from 40 to 100% in the reduced-dynamic mode. The geometry of GRACE FO-1 on August 14, 2018, was used for the tests, and the overbounding parameters of the MADOCA products in the average-case (Table 2) were applied for the analysis.
From Fig. 5 it can be seen that the mean value of PLs increase, i.e., tending to reduced availability, in all three directions with decreasing duty cycles. This is due to the reduced model strength when decreasing the amount of data used for POD. However, the mean PLs after the increase are limited to about 2 m even when reducing the duty cycle to 40%. With an AL set to 4 m, the availabilities are still above 98% when the duty cycle is reduced to 40%.

Concluding remarks
Integrity monitoring (IM) is an essential task in guaranteeing the safety and reliability of a positioning service. This study proposed new IM strategies for GNSS-based near-realtime LEO POD, investigating how different error sources and model strength could influence their final PLs.
As a condition for high accuracy near-real-time LEO POD, the overbounding standard deviation and mean value were studied for the SISRE of real-time GPS orbits and clocks using the IGS RTS and MADOCA products. The investigation distinguished between the average-case and worst-case projections. Using the geometry of the LEO satellite GRACE FO-1, the results showed that for the averagecase assumption, the overbounding standard deviations are at the sub-dm to dm-level, and the overbounding mean values are within 5 cm, while in the worst-case the overbounding standard deviations are about 0.1 m and the corresponding mean values are at the sub-dm to dm-level.
The PLs were computed using both the real-time GNSS products and real data for the LEO satellite GRACE FO-1. Assuming a PMI of 10 −5 , the strong dynamic models and the improved model strength bring the mean PLs from about 3-4 m in the kinematic mode down to about 1 m in the reduced-dynamic mode for the average-case assumption in the radial, along-track and cross-track directions, where the PL safely bound all the orbital PEs. The worstcase projection leads to safer, but overall more conservative (pessimistic) PLs, having mean values of about 1-3 m in the reduced-dynamic mode and over 5 m in the kinematic mode. The IM availabilities in the reduced-dynamic mode and for the average-case projection are above 99% in all three directions when the AL is set to 2 m. The analysis further showed that the overbounding mean value of the SISRE plays a determining role in the final PLs. For small satellites, with duty cycles down to 40%, the mean values of PLs (in the reduced-dynamic mode and the average-case projection) increase compared to those with a 100% duty cycle, but are limited to about 2 m using the overbounding parameters of the MADOCA products. The corresponding IM availabilities are shown to be above 98% when the AL is set to 4 m, even with the duty cycle reduced to 40%.