A Sensitivity Study of POD Using Dual-Frequency GPS for CubeSats Data Limitation and Resources

Making use of dual-frequency (DF) global navigation satellite system (GNSS) observations and good dynamic models, the precise orbit determination (POD) for the satellites on low earth orbits has been intensively investigated in the last decades and has achieved an accuracy of centimeters. With the rapidly increasing number of the CubeSat missions in recent years, the POD of CubeSats were also attempted with combined dynamic models and GNSS DF observations. While comprehensive dynamic models are allowed to be used in the postprocessing mode, strong constraints on the data completeness, continuity, and restricted resources due to the power and size limits of CubeSats still hamper the high-accuracy POD. An analysis of these constraints and their impact on the achievable orbital accuracy thus needs to be considered in the planning phase. In this study, with the focus put on the use of DF GNSS data in postprocessing CubeSat POD, a detailed sensitivity analysis of the orbital accuracy was performed w.r.t. the data continuity, completeness, observation sampling interval, latency requirements, availability of the attitude information, and arc length. It is found that the overlapping of several constraints often causes a relatively large degradation in the orbital accuracy, especially when one of the constraints is related to a low duty-cycle of, e.g., below 40% of time. Assuming that the GNSS data is properly tracked except for the assumed constraints, and using the International GNSS Service (IGS) final products or products from the IGS real-time service, the 3D orbital accuracy for arcs of 6 h to 24 h should generally be within or around 1 dm, provided that the limitation on data is not too severe, i.e., with a duty-cycle not lower than 40% and an observation sampling interval not larger than 60 s.


Introduction
Combining strong dynamic models and observations of global navigation satellite system (GNSS) collected onboard, the processing strategies used for the low Earth orbit (LEO) reduced-dynamic precise orbit determination (POD) have been intensively investigated in the last decades [1][2][3][4][5]. Compared to the dynamic orbits based on the dynamic models only and the kinematic orbits based on the GNSS observations only, the reduced-dynamic orbit determination has the advantage that it is more robust against model deficiencies for the former case, and possible data constraints, e.g., data in-continuity for the latter case. Following the need for the POD by a range of applications, such as the studies in the Earth gravity field [6,7], the atmospheric sounding [8], the altimetry [9], and the Interferometric Synthetic Aperture Radar (InSAR) [10], the LEO orbital accuracy based on undifferenced GNSS observations can reach a few centimeters [11,12], and even better at mm-level in case of baseline processing for formation-flying [12][13][14]. For undifferenced POD, the highest accuracy is often realized in postprocessing mode applying good dynamic models, precise GNSS satellite orbits, and clocks, and under the condition that the GPS dual-frequency phase and code data are tracked and collected with a good continuity and completeness. The antenna attitude information is often well monitored by extra sensors like star cameras and inertial measurement unit (IMU), and eventually will be used to correctly apply the antenna sensor offsets, the antenna phase centre offsets (PCOs), and variations (PCVs) [15].
The number of CubeSat missions has dramatically increased in this decade. By April 2020, there are over 1000 CubeSats successfully launched into their low earth orbits [16]. In addition to educational purposes, CubeSats are nowadays also used for the demonstration of new technologies and scientific researches [17,18]. Compared to large or medium LEO satellites, the CubeSats are mostly cost-effective and occupy a much smaller size and weight, i.e., typically from below 1 kg (0.25U) to more than 20 kg (12U). This often leads to smaller battery volume and smaller size of the solar panels. In addition to the two-line elements (TLEs), other sensors like sun sensors and magnetometers have been used for orbit determination of CubeSats with relatively low accuracy, i.e., at km-level [19]. In the case of higher positioning accuracy requirements, the GNSS-based processing has been demonstrated as a useful method for the CubeSat POD, for which successful postprocessed position fix was achieved with a meter-level accuracy [20]. While the onboard positioning normally utilizes simplified dynamic models due to the limitation of the computational load [21,22], the data postprocessing allows the usage of comprehensive dynamic models and the estimation of enough additional dynamic parameters to compensate remaining model deficiencies. However, even for the postprocessing mode, strongly constrained satellite size and power need to be shared among all sensors onboard. The capacity for data tracking, storage, and transfer is thus restricted, which limits the accuracy of orbit determination. The GNSS satellites may not be able to be fully tracked continuously due to these limitations, which could yield a duty-cycle smaller than 100% [20,21,23], and lead to fewer data available for the processing and more ambiguities to be set up. The size of the data collected between each data dump could also be limited [20]. The collected GNSS observations are most likely to be temporarily stored onboard and transferred back to the ground processing centre (GPC) every orbital cycle or a few cycles during a limited ground contact time, while information of tasks other than positioning is required to be stored and transferred at the same time. This may implicitly limit the GNSS data amount during the downlink. Due to these constraints and the potential limit on computational load in case of the real-time onboard processing, a relatively low sampling rate might be preferable for CubeSat POD, or low in general but only high in specific areas [24]. Furthermore, one may need to prepare for a low number of tracked satellites in case of assigning channels for the purpose of initial signal acquisitions [20] and in case of improper data tracking due to possible antenna control problems, which will be discussed later in this section.
In addition to the data continuity and availability, other information required for the POD may also be limited for CubeSat missions. Due to the limitation of the CubeSats in terms of their size, power, and capacity to install extra sensors, the antenna attitude information may not be well measured, estimated, documented, and transferred back to the GPC with sufficient details. For example, due to power management needs, the IMU under its own high-accuracy experiment [25] onboard the CubeSat Simulation To Flight-1 (STF-1), investigated by NASA's Katherine Johnson Independent Verification and Validation (IV&V) Facility and West Virginia University [21,26], was unable to be activated during the GNSS data collection. Depending on the latency requirements of the CubeSat orbits, the high-accuracy GNSS satellite orbits and clocks might not be available for the processing, e.g., the international GNSS service (IGS) final products [27,28], which has a latency of 12-18 days. As explained in [11], this could lead to dramatic accuracy differences already in case of good data continuity and availability, mainly due to the decreased sampling rates of the satellite clocks. Moreover, instead of processing orbits with long arcs of 24 h or 12 h (which accounts for about 16 and 8 orbital cycles for an orbital period of 1.5 h), rapid processing may also be used to compute the orbit with shorter arcs of several hours. While some conditions mentioned above may not significantly influence the orbital accuracy when individually being considered, large degradation of different levels could happen when several conditions happen together. Therefore, at the planning phase of a CubeSat mission, it is important to assess the orbital accuracy under different scenarios, so that one can properly compromise between the limits on diverse resources and the expected orbital accuracy.
The GNSS-based orbit determination of small satellites is often performed with single-frequency signals [23,24,29,30], which leads to an accuracy of decimeters to meters depending on the data availability and the dynamic models used. In recent years, to enable a higher POD accuracy, the usage of dual-frequency code and phase GPS data has been attempted and planned in CubeSat missions, even for a multi-constellation scenario [31]. Having dual-frequency phase and code observations at hand, the first-order ionospheric delays, which amounts to about 99% of the total ionospheric delays [32], can be eliminated by forming the ionosphere-free (IF) combination instead of the group and phase ionospheric correction (GRAPHIC) combination, which is often formed in the single-frequency case involving both the code and phase signals. While having the second frequency should theoretically enables better POD accuracy, a higher challenge for data tracking, storage, and transfer is also associated with the use of the dual-frequency GNSS in CubeSat missions. In addition to the planned constraints on diverse resources, the actual data collection of the dual-frequency CubeSat missions may not fulfil the expectation due to the lower robustness compared to the large LEO satellites. As an example, Figure 1 shows the real data condition of the 3U CubeSat STF-1 experimented by NASA's Katherine Johnson IV&V Facility and West Virginia University, which was expected to collect dual-frequency GPS data with a duty-cycle of 100% during a specific 12 h period. The actual data condition, however, did not meet the expectation, possibly due to the lack of the antenna attitude control. Due to the limitations in the antenna system, reports of short data pieces and frequent cycle slips, which hampers the high-accuracy POD, can also be found for the triple-CubeSat CanX-2 [20]. As such, not only for the planning purpose but also to get prepared for the possible sub-optimal real data collection, investigations of the CubeSat POD accuracy with diverse constraints on data and resources are important. A proper compromise is very possibly needed between the limiting resources, the suboptimal data conditions, and the requirements on orbital accuracy.
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 21 so that one can properly compromise between the limits on diverse resources and the expected orbital accuracy. The GNSS-based orbit determination of small satellites is often performed with single-frequency signals [23,24,29,30], which leads to an accuracy of decimeters to meters depending on the data availability and the dynamic models used. In recent years, to enable a higher POD accuracy, the usage of dual-frequency code and phase GPS data has been attempted and planned in CubeSat missions, even for a multi-constellation scenario [31]. Having dual-frequency phase and code observations at hand, the first-order ionospheric delays, which amounts to about 99% of the total ionospheric delays [32], can be eliminated by forming the ionosphere-free (IF) combination instead of the group and phase ionospheric correction (GRAPHIC) combination, which is often formed in the single-frequency case involving both the code and phase signals. While having the second frequency should theoretically enables better POD accuracy, a higher challenge for data tracking, storage, and transfer is also associated with the use of the dual-frequency GNSS in CubeSat missions. In addition to the planned constraints on diverse resources, the actual data collection of the dual-frequency CubeSat missions may not fulfil the expectation due to the lower robustness compared to the large LEO satellites. As an example, Figure 1 shows the real data condition of the 3U CubeSat STF-1 experimented by NASA's Katherine Johnson IV&V Facility and West Virginia University, which was expected to collect dual-frequency GPS data with a duty-cycle of 100% during a specific 12 h period. The actual data condition, however, did not meet the expectation, possibly due to the lack of the antenna attitude control. Due to the limitations in the antenna system, reports of short data pieces and frequent cycle slips, which hampers the high-accuracy POD, can also be found for the triple-CubeSat CanX-2 [20]. As such, not only for the planning purpose but also to get prepared for the possible sub-optimal real data collection, investigations of the CubeSat POD accuracy with diverse constraints on data and resources are important. A proper compromise is very possibly needed between the limiting resources, the suboptimal data conditions, and the requirements on orbital accuracy. For the dual-frequency CubeSat POD, which is the focus of this study, it was pointed out by previous studies that the discontinuity of the observation data with a low duty-cycle may significantly degrade the orbital accuracy in both the real-time and the postprocessing modes [21,33]. In this study, we restrict focus on postprocessing reduced-dynamic dual-frequency CubeSat POD, which relieves the strong constraints of the computational load onboard the CubeSat and allows for the usage of good dynamic models. In the experiments of CanX-2 [20] and STF-1 [21], geodetic grade For the dual-frequency CubeSat POD, which is the focus of this study, it was pointed out by previous studies that the discontinuity of the observation data with a low duty-cycle may significantly degrade the orbital accuracy in both the real-time and the postprocessing modes [21,33]. In this study, we restrict focus on postprocessing reduced-dynamic dual-frequency CubeSat POD, which relieves the strong constraints of the computational load onboard the CubeSat and allows for the usage of good dynamic models. In the experiments of CanX-2 [20] and STF-1 [21], geodetic grade receivers NovAtel OEM4-G2 and OEM615 tracking dual-frequency GPS signals were mounted onboard the CubeSats. With the expectation that in future CubeSat missions, dual-frequency measurements can be properly tracked in each duty-cycle by geodetic receivers with a mean average satellite number of at least 5 to 6, this study performs a sensitivity analysis of different limitations on the data availability, continuity, resources, latencies, and arc lengths with respect to the orbital accuracy. The overlapping effect of different constraints is investigated, and a compromise between these factors and the orbital accuracy is discussed. In addition to CubeSats, the analysis should also benefit the POD of other small satellites that need to bare similar constraints on data, power, and resources.

Processing Strategy
In this study, the final LEO orbits (where CubeSats are classified as LEO) are postprocessed following the reduced-dynamic LEO processing scheme of the Bernese GNSS Software V5.2 [34]. After correcting the antenna PCOs, PCVs, and the offsets between the antenna reference point (ARP) and the satellite center-of-mass (CoM), the determined final orbit in this study is referred to the satellite CoM. For the sensitivity analysis purpose, the final reduced-dynamic orbits are produced assuming different observation sampling rates, duty cycling, varying mean number of satellites by changing the elevation mask as will be explained in the next section, arc lengths, using different types of the IGS products, with and without the antenna attitude information. The processing procedure can generally be summarized in the following steps, which is also given in Figure 2 as a flow diagram:

1.
Compute kinematic orbits using single point positioning (SPP) employing the IF combination of the code observations. These kinematic orbits, denoted as vectorr K , are discrete and have an accuracy of meters.

2.
Computation of the code-based reduced-dynamic orbits. The reduced-dynamic orbits are computed with accelerations based on a series of gravitational and nongravitational terms, such as the Earth gravitational terms, the Earth tidal terms, the gravitational attraction from the sun, moon, and other planets, as well as the general relativistic term. Note that mis-modeled effects like the solar radiation pressure and the air drag will be largely absorbed by the estimated dynamic parameters and the stochastic velocity changes or accelerations set up later in the processing [35]. Details of the processing and the dynamic models are given in Table 1. Making use of the kinematic code orbits from the first step, the six Keplerian elements at the initial condition (the semi-major axis of the orbit, the orbital eccentricity, the inclination of the orbital plane, the right ascension of the ascending node, the argument of perigee, and the argument of latitude at the initial condition), and a remaining part of the dynamic models are estimated with a batch least-squares adjustment, which includes at this step nine flight-oriented dynamic parameters. These estimable dynamic parameters contain three constant terms (a R0 , a S0 , and a W0 ) and six periodic terms (a RC , a SC , a WC , a RS , a SS , and a WS ) in the radial (R), along-track (S) and cross-track (W) directions. The total acceleration a can then be distributed into the term a 0 , which is assumed known by applying the models given in Table 1, and an additional dynamic term a dyn that is to be adjusted: With a dyn = a R e R + a S e S + a W e W , a R = a R0 + a RC cos(U) + a RS sin(U), a S = a S0 + a SC cos(U) + a SS sin(U), where e R , e S , and e W represent the unit vectors in the radial, along-track, and cross-track directions, respectively. U denotes the satellite argument of latitude. The code-based kinematic orbitsr K obtained from the first step are used as observations to adjust the 15 parameters mentioned above in a least-squares sense. The linearized observation equation at the epoch t i can be formulated as: wherer 0 is the a priori orbit vector obtained based on numerical integration on hand the a 0 , and theâ dyn and the Keplerian elements estimated from the last iteration. The vectors x k and x d contain the increments of the six Keplerian elements and the nine dynamic parameters, respectively, and the design matrices A rk and A rd contain the partial derivatives of the position vectors with respect to x k and x d . The partial derivatives are computed with numerical integration of the variational equations [36,37]. E is the expectation operator. The reduced-dynamic orbit can be interpolated for time epochs with higher sampling rates and produced for periods with data gaps.
In this study, all orbits are resampled into time epochs with 10 s sampling interval for assessment, regardless of which duty-cycles and observation sampling rates are applied in the processing.

3.
Phase preprocessing and orbit improvements. This step preprocesses the raw phase observations to detect cycle slips and mark bad observations. The preprocessing goes through several iterations to improve the LEO orbit quality. The orbit improvement is realized through estimation of stochastic velocity changes [38] in addition to the 15 parameters mentioned in the second step, and is performed in a least-squares adjustment making use of the IF combination of the phase observations. One set of the stochastic velocity changes is considered in each predefined time interval, e.g., every 15 min. The linearized phase observation equation at the epoch t i can be expressed as: With where ∆ϕ IF represents the observed-minus-computed (O-C) term of the IF phase observations. c and ∆t r denote the speed of light and the receiver clock error, respectively. λ j , f j , and N j represent the wavelength, the frequency, and the ambiguity on frequency j (j = 1, 2), respectively. Note that N I F is not an integer. The receiver clock error is estimated epoch-wise independently, and the ambiguity is assumed constant before the detection of a cycle slip. Note that new ambiguities are setup for estimation at the beginning of each round of duty cycling. x v stands for the vector containing all stochastic velocity changes in the RSW directions from the first to the current epoch and note that x v is constrained to zero with a predefined a priori standard deviation. The design matrices A lk , A ld , and A lv contain the partial derivatives of the O-C terms with respect to the x k , x d , and x v , respectively. To be estimated are the vector [x k , x d , x v ]T, the receiver clock offset ∆t r , and the term N IF . Note that very little code observations are used to avoid the problem of matrix singularity between the receiver clock offset and the ambiguity terms. Note that the ambiguities on L1 and L2 are not attempted to be fixed in this study.

4.
Generation of final orbits. With the preprocessed phase observations, the six Keplerian elements, the three constant dynamic parameters (a R0 , a S0 , and a W0 ) are estimated together with stochastic accelerations in the RSW directions. The accelerations are set up in shorter time intervals compared to those in Step 3. The linearized phase observation equation is thus formulated as: where x d0 and x a denote the increment vector of the three constant dynamic parameters and all the stochastic accelerations from the first to the current epoch, respectively. x a is constrained to zero with a predefined a priori standard deviation (selected as 5 × 10 −9 m/s 2 in all the three directions based on the default setting for GRACE satellites in the Bernese software, as real data from the GRACE Follow-on mission is used for test purposes in this study, which will be explained later. This value may vary for satellites of other missions). A ld0 and A la correspond to the partial derivatives of the phase O-C terms with respect to the x d0 and x a , respectively. Note that very little code observations are used to avoid singularity mentioned above.
where 0 and denote the increment vector of the three constant dynamic parameters and all the stochastic accelerations from the first to the current epoch, respectively.
is constrained to zero with a predefined a priori standard deviation (selected as 5× 10 −9 m/s 2 in all the three directions based on the default setting for GRACE satellites in the Bernese software, as real data from the GRACE Follow-on mission is used for test purposes in this study, which will be explained later. This value may vary for satellites of other missions). 0 and correspond to the partial derivatives of the phase O-C terms with respect to the 0 and , respectively. Note that very little code observations are used to avoid singularity mentioned above.

Coordinate Transformation
Nutation and precession: IAU2000R06 [34] Sub-daily pole variations: IERS Conventions 2010 [41] Earth rotation parameters: IGS final, rapid, ultra-rapid products In this study, the data from the LEO satellite of the GRACE Follow-on mission [43], the GRACE FO-1, which continuously tracked dual-frequency code and phase GPS observations on L1 and L2, are used for different scenarios simulating conditions with limited data and resources on 22 August 2018. For the purpose of consistency, the processed orbits are differenced with the best-case scenario in this study (as the reference for comparison), i.e., the phase-based reduced dynamic orbits with a duty-cycle of 100%, an observation sampling interval of 10 s, an elevation mask of 5 degrees (default setting in Bernese for LEO processing), using the IGS final products and processed in a 24 h arc. The use of our best-case scenario as a reference for comparing different tested cases is justified by comparing it with the reference orbits provided by the Jet Propulsion Laboratory (JPL) [44] for the same day, where good consistency with a 3-dimentional (3D) root-mean-squared error (RMSE) of below 2 cm was achieved. Nevertheless, it should be noted that in this study, the 3D RMSE computed under different scenarios refer to the accuracy difference compared with the best-case scenario, but not the absolute orbital accuracy. Note that for LEO missions, independent technics like the satellite laser ranging (SLR) and the K-band ranging (KBR) are often utilized to validate the accuracy of reference orbits. As a large number of variants with respect to different observation sampling rates, duty-cycles, satellite numbers, latencies, availability of the attitude information, and arc lengths are tested, this study does not attempt to study long-term variation of results. As a consistency check, the POD results of the best-case scenario from 14 to 21 August 2018 were compared with the JPL reference orbits, and no significant difference from those of the test day in this study can be observed, i.e., at a few millimeters. As such, we consider the simulations based on the test day are representative.

Orbit Determination under Different Scenarios
As a start point of our analysis, let us look at the results of code-based processing and see the impact of using a reduced-dynamic model. Figure 3 shows the 3D positional errors of the kinematic orbits (red) based on the code observations on P1 and P2 forming the IF combination (Step 1 in Section 2), and when computing the same data with the reduced-dynamic orbits (blue, Step 2 in Section 2) for GRACE FO-1 for the tested 1 day of data. As complete dual-frequency CubeSats data for a long period needed in our analysis are currently not available, the data from a typical LEO satellite that may share similar orbit conditions with CubeSats, in this study GRACE FO-1, were used for simulation of the best-case scenario of the CubeSat processing [33]. For the tested LEO satellite, the dual-frequency GPS data are tracked in a complete and continuous manner with an observation sampling interval of 10 s. A duty-cycle of 100% and a mean satellite number of about 9 above the elevation mask of 5 degrees were available for the processing of a daily arc using the IGS final products. The 3D RMSE of the code-based kinematic (Step 1 in Section 2) and reduced-dynamic orbits (Step 1 in Section 2) amount to about 1.6 m and 0.2 m, respectively, which clearly shows the value of employing strong dynamic models. When reducing the data availability and continuity, the accuracy of the reduced-dynamic orbit could degrade to the sub-meter level.
For high-accuracy POD, as discussed in Section 2, the code-based orbits are only used as the a priori orbits for the phase processing. Therefore, in the rest of the sections, the sensitivity analysis will focus on the final phase-based reduced-dynamic orbits.

Duty Cycling and Satellite Numbers
As mentioned before, due to the constraints of the CubeSats on the power, the data storage, and the time of active attitude control possibly during the sunlight, the GNSS observations might have breaks, and are collected in pieces of different length [20]. In this study, observation data are assumed to be collected with different duty-cycles, i.e., with the power turned on only for a part of the time. The duty-cycle is defined here as a percentage of time data was available out of the total time considered. The sensitivity analysis is performed assuming the varying from 20% to 100%. In our tests, a duty-cycle implies that the data are available in the first 60 × min of each hour. The 3D orbital errors are shown in Figure 4 with varying . The observation sampling interval was 10 s. Recall that the phase-based reduced-dynamic orbit with a duty-cycle of 100% was our best-case scenario and was used as the reference for comparison as mentioned in Section 2. The RMSE increases from about 1.6 cm with a of 80% to about 3.5 cm with a of 20%. Note that the RMSE is computed from the beginning of the first power-on period to the end of the last power-on period of the test day. As shown in Figure 1 and discussed in [19], the number of satellites that can be properly tracked (at least for the P1 and P2 observations) by the CubeSat receiver could be lower than expected due to

Duty Cycling and Satellite Numbers
As mentioned before, due to the constraints of the CubeSats on the power, the data storage, and the time of active attitude control possibly during the sunlight, the GNSS observations might have breaks, and are collected in pieces of different length [20]. In this study, observation data are assumed to be collected with different duty-cycles, i.e., with the power turned on only for a part of the time. The duty-cycle D is defined here as a percentage of time data was available out of the total time considered. The sensitivity analysis is performed assuming the D varying from 20% to 100%. In our tests, a duty-cycle implies that the data are available in the first 60 × D min of each hour. The 3D orbital errors are shown in Figure 4 with varying D. The observation sampling interval was 10 s. Recall that the phase-based reduced-dynamic orbit with a duty-cycle of 100% was our best-case scenario and was used as the reference for comparison as mentioned in Section 2. The RMSE increases from about 1.6 cm with a D of 80% to about 3.5 cm with a D of 20%. Note that the RMSE is computed from the beginning of the first power-on period to the end of the last power-on period of the test day.

Duty Cycling and Satellite Numbers
As mentioned before, due to the constraints of the CubeSats on the power, the data storage, and the time of active attitude control possibly during the sunlight, the GNSS observations might have breaks, and are collected in pieces of different length [20]. In this study, observation data are assumed to be collected with different duty-cycles, i.e., with the power turned on only for a part of the time. The duty-cycle is defined here as a percentage of time data was available out of the total time considered. The sensitivity analysis is performed assuming the varying from 20% to 100%. In our tests, a duty-cycle implies that the data are available in the first 60 × min of each hour. The 3D orbital errors are shown in Figure 4 with varying . The observation sampling interval was 10 s. Recall that the phase-based reduced-dynamic orbit with a duty-cycle of 100% was our best-case scenario and was used as the reference for comparison as mentioned in Section 2. The RMSE increases from about 1.6 cm with a of 80% to about 3.5 cm with a of 20%. Note that the RMSE is computed from the beginning of the first power-on period to the end of the last power-on period of the test day. As shown in Figure 1 and discussed in [19], the number of satellites that can be properly tracked (at least for the P1 and P2 observations) by the CubeSat receiver could be lower than expected due to As shown in Figure 1 and discussed in [19], the number of satellites that can be properly tracked (at least for the P1 and P2 observations) by the CubeSat receiver could be lower than expected due to different reasons, such as the antenna stabilization problems or specific channel assignment. Therefore, to resemble this practical issue, a simplified approach is performed in this simulation to reduce the number of satellites observed under good conditions by excluding satellites with low elevation angles. The elevation mask is increased from 5 to 25 degrees in the test, resulting in a reduction of the mean satellite number from about 9 to 6 (see Table 2). The Figure 5a illustrates the distribution of the satellite numbers applying different elevation masks for GRACE FO-1 for the test period. A gradual reduction of the satellite numbers can be observed when increasing the elevation mask. The percentile of valid SPP solutions is reduced correspondingly from nearly 100% to about 83% due to the reduction in redundancy. different reasons, such as the antenna stabilization problems or specific channel assignment. Therefore, to resemble this practical issue, a simplified approach is performed in this simulation to reduce the number of satellites observed under good conditions by excluding satellites with low elevation angles. The elevation mask is increased from 5 to 25 degrees in the test, resulting in a reduction of the mean satellite number from about 9 to 6 (see Table 2). The Figure 5a illustrates the distribution of the satellite numbers applying different elevation masks for GRACE FO-1 for the test period. A gradual reduction of the satellite numbers can be observed when increasing the elevation mask. The percentile of valid SPP solutions is reduced correspondingly from nearly 100% to about 83% due to the reduction in redundancy.  Having a duty-cycle of 100%, as shown in Figure 5b, decreasing the number of the received satellites from, e.g., 7 to 6, degrades the RMSE from about 1.1 to 2.1 cm. Note that the case with the mean satellite number of 9 having an elevation mask of 5 degrees is our best-case scenario and is used as the reference for comparison, in other words, the RMSE represent the solution discrepancy between the cases of using 7 or 6 satellites and when using 9 satellites.
Although the differences observed in Figures 4 and 5 are not large, coupling a reduction in the duty-cycle with a reduction in the number of received satellites would cause larger problems. Table 3 shows the 3D RMSE of the orbits applying different duty-cycles and mean satellite numbers, assuming an observation sampling interval of 10 s still. As mentioned earlier, the solution obtained with a duty-cycle of 100% and the mean satellite number of 9 is used as the reference for assessing the solutions by computing their discrepancies. It can be observed that with a duty-cycle equal or above 60%, even when the observed satellites are limited, the RMSE is still below 3 cm. With lower Having a duty-cycle of 100%, as shown in Figure 5b, decreasing the number of the received satellites from, e.g., 7 to 6, degrades the RMSE from about 1.1 to 2.1 cm. Note that the case with the mean satellite number of 9 having an elevation mask of 5 degrees is our best-case scenario and is used as the reference for comparison, in other words, the RMSE represent the solution discrepancy between the cases of using 7 or 6 satellites and when using 9 satellites.
Although the differences observed in Figures 4 and 5 are not large, coupling a reduction in the duty-cycle with a reduction in the number of received satellites would cause larger problems. Table 3 shows the 3D RMSE of the orbits applying different duty-cycles and mean satellite numbers, assuming an observation sampling interval of 10 s still. As mentioned earlier, the solution obtained with a duty-cycle of 100% and the mean satellite number of 9 is used as the reference for assessing the solutions by computing their discrepancies. It can be observed that with a duty-cycle equal or above 60%, even when the observed satellites are limited, the RMSE is still below 3 cm. With lower duty-cycles, however, observing all or almost all visible satellites would be essential to maintain a low RMSE at this level. With a duty-cycle of 20%, decreasing the mean satellite number to 6 would increase the RMSE to about 5 cm.

Sampling Rate of the Observations
In Section 3.1, the GPS observations onboard the LEO satellite with a sampling interval of 10 s were used for producing the reduced-dynamic orbits. Depending on the capacity of the CubeSat missions, the data might need to be stored and transferred to the GPC with a lower sampling rate. A compromise thus needs to be considered between the duty-cycles, the data completeness, and the sampling rate. Figure 6 illustrates the 3D RMSE of the orbits with the sampling interval varying from 10 to 120 s and the duty-cycle varying from 100% to 20%. The cases having mean satellite numbers of 9, 7, and 6 are given between each pair of dashed lines. Recall that the number of satellites is controlled by adjusting the elevation mask (Section 3.1). In the case of a duty-cycle of 20%, the orbits generated with an observation sampling interval of 120 s (red line) are not produced due to the extremely low number of continuous observations in each data patch. Note that all the RMSE are computed based on orbits resampled to 10 s sampling intervals.
duty-cycles, however, observing all or almost all visible satellites would be essential to maintain a low RMSE at this level. With a duty-cycle of 20%, decreasing the mean satellite number to 6 would increase the RMSE to about 5 cm.

Sampling Rate of the Observations
In Section 3.1, the GPS observations onboard the LEO satellite with a sampling interval of 10 s were used for producing the reduced-dynamic orbits. Depending on the capacity of the CubeSat missions, the data might need to be stored and transferred to the GPC with a lower sampling rate. A compromise thus needs to be considered between the duty-cycles, the data completeness, and the sampling rate. Figure 6 illustrates the 3D RMSE of the orbits with the sampling interval varying from 10 to 120 s and the duty-cycle varying from 100% to 20%. The cases having mean satellite numbers of 9, 7, and 6 are given between each pair of dashed lines. Recall that the number of satellites is controlled by adjusting the elevation mask (Section 3.1). In the case of a duty-cycle of 20%, the orbits generated with an observation sampling interval of 120 s (red line) are not produced due to the extremely low number of continuous observations in each data patch. Note that all the RMSE are computed based on orbits resampled to 10 s sampling intervals. From Figure 6 it can be observed that the values in general gradually increase from the left to the right side, i.e., from the cases with good data continuity, completeness, to those with bad data conditions. With less data available for the processing in the Steps 3 and 4 (Section 2), the estimable parameters, including the six Keplerian elements, the additional dynamic parameters, and the stochastic parameters, are determined with lower precision. This results in a lower accuracy of the numerically integrated orbits. At the same time, it can also be observed that the differences between having observation sampling intervals of 10 s (green), 20 s (yellow), and 30 s (magenta) are small. From Figure 6 it can be observed that the values in general gradually increase from the left to the right side, i.e., from the cases with good data continuity, completeness, to those with bad data conditions. With less data available for the processing in the Steps 3 and 4 (Section 2), the estimable parameters, including the six Keplerian elements, the additional dynamic parameters, and the stochastic parameters, are determined with lower precision. This results in a lower accuracy of the numerically integrated orbits. At the same time, it can also be observed that the differences between having observation sampling intervals of 10 s (green), 20 s (yellow), and 30 s (magenta) are small. This can be attributed to the 30 s GPS clocks provided in the IGS final products, where interpolation into high-sampled time epochs would introduce biases within a higher percentage of the processing time (to be discussed in the next subsection). Provided that the receiver has a high duty-cycle not lower than 60%, even with a limited number of observed satellites and a high sampling interval of 120 s, the 3D RMSE is still below or around 3 cm. However, with lower duty-cycles, the degradation of the orbit accuracy could be dramatic when the observations are only available with a low sampling rate and when the observed satellites are limited.

Latency Applying Different GPS Products
In Sections 3.1 and 3.2, the IGS final products for GPS satellite orbits and clocks were used for the estimation of the reduced-dynamic orbits. The IGS final products are produced weekly and have a latency of 12 to 18 days. In case that the CubeSat orbits are required within a shorter time, the IGS (or one of its analysis centres) rapid products or products from the IGS real-time service (RTS) [45,46] can be used instead. Table 4 gives the latencies of different IGS products. In this study, the decoded IGS real-time stream IGC01, which is denoted as the IGC products and available at [47], is used for simulation of the scenario retrieving the GPS products in (near)-real-time. Products of the single IGS analysis centers, e.g., the real-time products from the National Centre for Space Studies (CNES), are also of high quality [48]. In this study, for consistency purpose, we used the IGS products of different latencies and satellite clock sampling rates.  Table 4 it can also be seen that the sampling rate of the IGR satellite clocks is lower than the other two products. In [11], the sampling interval of GPS satellite clock-offset corrections was identified as a limiting factor for the LEO orbital accuracy. This possibly explains the fact that the 3D orbital RMSE using the IGR products are not better than those using the IGC products. Figure 7 shows the differences of the between-satellite clock corrections for G02 and G32 between the IGS final and rapid clocks. The between-satellite clock differences are shown to remove possible inconsistencies between the reference time scales of different products. The blue dots indicate the differences at the given sampling interval of the rapid clocks, i.e., 300 s. The red dots illustrate the differences at a sampling interval of 30 s, where the rapid clocks were linearly interpolated using the nearest epochs to this higher sampling rate. From Figure 7 it can be observed that although the differences are generally at centimeters for the blue dots with the given sampling rate of the rapid clocks, after interpolation, the differences are increased to dm-level (see the red dots). For the test day, the average RMS of the clock differences for all satellites between the IGS and IGR final products, divided by √ 2 for the between-satellite differencing, amount to about 2 cm. Note that the large pattern of the red dots is caused by the stochastic behaviors of the clocks, and an interpolation over a longer period with a higher-order polynomial does not help to reduce the biases. As an example, interpolation with a quadratic polynomial within each half hour increases the average RMS by several millimeters. Assuming a duty-cycle of 100%, an observation sampling interval of 10 s, and a mean satellite number of 9, the last column of Table 4 shows the 3D RMSE of the final LEO orbits determined using different IGS products compared with the best-case scenario, i.e., using the IGS final products. A 3D RMSE of several centimeters can be achieved for the IGR and IGC products. This is further illustrated by Figure 8, which shows the time series of the 3D orbital errors applying the IGR and IGC products. In the case of having limited data available, i.e., with short duty-cycles, limited satellite numbers and sampling rate, the orbital accuracy using the IGR and IGC products would not be as good as shown in Figure 8 anymore. In the case of using the IGR and IGC products, Figure 9 illustrated the 3D RMSE of the difference between phase-based reduced-dynamic orbits and the best-case scenario when applying different duty-cycles, mean satellite numbers, and observation sampling rates. Note that the cases having different satellite numbers (9, 7, and 6) are given between each pair of the dashed lines. Compared to Figure 6 using the IGS final products, the 3D RMSE has generally increased due to the reduced accuracy of the GPS products and the higher sampling interval of the GPS clocks (for IGR products). As discussed for the IGS final products, when using the IGC products (Figure 9b), one can similarly observe that a high observation sampling interval of, e.g., 10 s (green) or 20 s (yellow), is not helpful to reduce the RMSE due to the GPS clock sampling of 30 s. For the cases when using Assuming a duty-cycle of 100%, an observation sampling interval of 10 s, and a mean satellite number of 9, the last column of Table 4 shows the 3D RMSE of the final LEO orbits determined using different IGS products compared with the best-case scenario, i.e., using the IGS final products. A 3D RMSE of several centimeters can be achieved for the IGR and IGC products. This is further illustrated by Figure 8, which shows the time series of the 3D orbital errors applying the IGR and IGC products. Assuming a duty-cycle of 100%, an observation sampling interval of 10 s, and a mean satellite number of 9, the last column of Table 4 shows the 3D RMSE of the final LEO orbits determined using different IGS products compared with the best-case scenario, i.e., using the IGS final products. A 3D RMSE of several centimeters can be achieved for the IGR and IGC products. This is further illustrated by Figure 8, which shows the time series of the 3D orbital errors applying the IGR and IGC products. In the case of having limited data available, i.e., with short duty-cycles, limited satellite numbers and sampling rate, the orbital accuracy using the IGR and IGC products would not be as good as shown in Figure 8 anymore. In the case of using the IGR and IGC products, Figure 9 illustrated the 3D RMSE of the difference between phase-based reduced-dynamic orbits and the best-case scenario when applying different duty-cycles, mean satellite numbers, and observation sampling rates. Note that the cases having different satellite numbers (9, 7, and 6) are given between each pair of the dashed lines. Compared to Figure 6 using the IGS final products, the 3D RMSE has generally increased due to the reduced accuracy of the GPS products and the higher sampling interval of the GPS clocks (for IGR products). As discussed for the IGS final products, when using the IGC products (Figure 9b), one can similarly observe that a high observation sampling interval of, e.g., 10 s (green) or 20 s (yellow), is not helpful to reduce the RMSE due to the GPS clock sampling of 30 s. For the cases when using In the case of having limited data available, i.e., with short duty-cycles, limited satellite numbers and sampling rate, the orbital accuracy using the IGR and IGC products would not be as good as shown in Figure 8 anymore. In the case of using the IGR and IGC products, Figure 9 illustrated the 3D RMSE of the difference between phase-based reduced-dynamic orbits and the best-case scenario when applying different duty-cycles, mean satellite numbers, and observation sampling rates. Note that the cases having different satellite numbers (9, 7, and 6) are given between each pair of the dashed lines. Compared to Figure 6 using the IGS final products, the 3D RMSE has generally increased due to the reduced accuracy of the GPS products and the higher sampling interval of the GPS clocks (for IGR products). As discussed for the IGS final products, when using the IGC products (Figure 9b), one can similarly observe that a high observation sampling interval of, e.g., 10 s (green) or 20 s (yellow), is not helpful to reduce the RMSE due to the GPS clock sampling of 30 s. For the cases when using the IGR products (Figure 9a) with a higher GPS clock sampling interval of 300 s (Table 4), no obvious degradation can be concluded even when having an observation sampling interval of 120 s.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 21 the IGR products ( Figure 9a) with a higher GPS clock sampling interval of 300 s (Table 4), no obvious degradation can be concluded even when having an observation sampling interval of 120 s.
(a) (b) Figure 9. 3D RMSE of the reduced-dynamic orbits under different duty-cycles, observation sampling intervals and mean satellite numbers using (a) the IGR and (b) the IGC products. The cases having 9, 7, and 6 satellites (on average) are given between each pair of the dashed lines.

Antenna Attitude
As shown in the example in Figure 1, it is suspected that the poor data tracking condition is caused by the lack of antenna stabilization. Good antenna stabilization is essential to keep the antenna in a favorable orientation for GNSS tracking, i.e., approximately toward the radial direction. Highaccuracy attitude control is expected to be realized also for CubeSats, as found to be feasible in the tests simulating different CubeSat scenarios in [50] and achieved for the CubeSat CanX-2 [20]. Under the condition of appropriate antenna stabilization, the antenna attitude information is often measured by extra sensors like star cameras or IMU on LEO/CubeSats satellites. For CubeSats, however, the limited size and power may limit its capacity to properly measure, store, or transfer the attitude information back to the GPC. As mentioned in Section 1, the data from the IMU onboard the CubeSat STF-1, e.g., was not collected during the GNSS data tracking because of the power allocation, where the IMU used was of a special type tested for a high-accuracy experiment that was carried out at a different time. The attitude information of the satellite (and antenna) is sometimes even ignored in CubeSat simulation studies [21]. Under such circumstances, where no attitude information is available, the antenna attitude could be estimated based on the satellite positions and velocities. Taking the GRACE FO-1 as an example, the antenna was fixed on the satellite and its up direction generally points toward the radial direction from the earth to the LEO satellite. Note that the X-axis in the satellite reference frame (SRF) is here approximately in the anti-flight direction, as GRACE FO-1 is the leading satellite in the test period. Figure 10 depicts the RSW system (green) in the radial, along-track and cross-track directions with the origin in the satellite centre-of-mass, and the antenna reference frame (ARF, red), e.g., for GRACE FO-1, in the north, east, and up (NEU) directions with the origin in the antenna reference point (ARP). To evaluate the differences between the expected and the actual orientations, the angle differences are computed with the help of the east ( ) and the opposite along-track ( ) unit vectors, the north ( ) and the cross-track ( ) unit vectors, and the up ( ) and the radial ( ) unit vectors in the inertial system: = arccos (− ), = arccos ( ), = arccos ( ).

Antenna Attitude
As shown in the example in Figure 1, it is suspected that the poor data tracking condition is caused by the lack of antenna stabilization. Good antenna stabilization is essential to keep the antenna in a favorable orientation for GNSS tracking, i.e., approximately toward the radial direction. High-accuracy attitude control is expected to be realized also for CubeSats, as found to be feasible in the tests simulating different CubeSat scenarios in [50] and achieved for the CubeSat CanX-2 [20]. Under the condition of appropriate antenna stabilization, the antenna attitude information is often measured by extra sensors like star cameras or IMU on LEO/CubeSats satellites. For CubeSats, however, the limited size and power may limit its capacity to properly measure, store, or transfer the attitude information back to the GPC. As mentioned in Section 1, the data from the IMU onboard the CubeSat STF-1, e.g., was not collected during the GNSS data tracking because of the power allocation, where the IMU used was of a special type tested for a high-accuracy experiment that was carried out at a different time. The attitude information of the satellite (and antenna) is sometimes even ignored in CubeSat simulation studies [21]. Under such circumstances, where no attitude information is available, the antenna attitude could be estimated based on the satellite positions and velocities. Taking the GRACE FO-1 as an example, the antenna was fixed on the satellite and its up direction generally points toward the radial direction from the earth to the LEO satellite. Note that the X-axis in the satellite reference frame (SRF) is here approximately in the anti-flight direction, as GRACE FO-1 is the leading satellite in the test period. Figure 10 depicts the RSW system (green) in the radial, along-track and cross-track directions with the origin in the satellite centre-of-mass, and the antenna reference frame (ARF, red), e.g., for GRACE FO-1, in the north, east, and up (NEU) directions with the origin in the antenna reference point (ARP). To evaluate the differences between the expected and the actual orientations, the angle differences are computed with the help of the east (e E ) and the opposite along-track (e S ) unit vectors, the north (e N ) and the cross-track (e W ) unit vectors, and the up (e U ) and the radial (e R ) unit vectors in the inertial system: Figure 10. The antenna reference frame (ARF, shown in red), the satellite reference frame (SRF, shown in yellow), and the radial, along-track, and cross-track (RSW, shown in green) system. 1 and 2 represent two conservative time epochs for the same satellite. The figure is scaled for a better presentation.
In Equation (11), the unit vectors and are calculated based on the satellite position and velocity vectors, and completes the left-hand system. As the antenna is fixed on the satellite, the antenna NEU unit vectors , , and are calculated based on the satellite orientations, which are tracked by extra sensors and documented in an attitude file with quaternions. The rotation matrix transforming from the SRF to the inertial system can be computed with these quaternions [51]. Based on the definition of the SRF as illustrated in Figure 10 in yellow, the , , and can be obtained with: The angle differences , , and are given in Figure 11. It can be observed that these angle differences are not large. In case that the antenna attitude file is not given, or the attitude information is not accurately documented, one could approximate the antenna orientation using the RSW directions. Figure 11. Angle differences between the east, north, and up directions in the antenna reference frame (ARF) and the opposite along-track direction, the cross-track direction, and the radial direction. The data for GRACE FO-1 was used for the plot. Note that the blue line is almost overwritten by the green dashed line. In Equation (11), the unit vectors e R and e S are calculated based on the satellite position and velocity vectors, and e W completes the left-hand system. As the antenna is fixed on the satellite, the antenna NEU unit vectors e N , e E , and e U are calculated based on the satellite orientations, which are tracked by extra sensors and documented in an attitude file with quaternions. The rotation matrix R SRF transforming from the SRF to the inertial system can be computed with these quaternions [51]. Based on the definition of the SRF as illustrated in Figure 10 in yellow, the e N , e E , and e U can be obtained with: The angle differences α E , α N , and α U are given in Figure 11. It can be observed that these angle differences are not large. In case that the antenna attitude file is not given, or the attitude information is not accurately documented, one could approximate the antenna orientation using the RSW directions. In Equation (11), the unit vectors and are calculated based on the satellite position and velocity vectors, and completes the left-hand system. As the antenna is fixed on the satellite, the antenna NEU unit vectors , , and are calculated based on the satellite orientations, which are tracked by extra sensors and documented in an attitude file with quaternions. The rotation matrix transforming from the SRF to the inertial system can be computed with these quaternions [51]. Based on the definition of the SRF as illustrated in Figure 10 in yellow, the , , and can be obtained with: The angle differences , , and are given in Figure 11. It can be observed that these angle differences are not large. In case that the antenna attitude file is not given, or the attitude information is not accurately documented, one could approximate the antenna orientation using the RSW directions. Figure 11. Angle differences between the east, north, and up directions in the antenna reference frame (ARF) and the opposite along-track direction, the cross-track direction, and the radial direction. The data for GRACE FO-1 was used for the plot. Note that the blue line is almost overwritten by the green dashed line. Figure 11. Angle differences between the east, north, and up directions in the antenna reference frame (ARF) and the opposite along-track direction, the cross-track direction, and the radial direction. The data for GRACE FO-1 was used for the plot. Note that the blue line is almost overwritten by the green dashed line.
Benefiting from the fact that the angle differences between the actual and the expected antenna orientation are small, even without having the attitude information available, differences less than 10 mm can mostly be observed in the final orbits provided that the observation data is complete and continuous (see Figure 12), i.e., with a 3D RMSE amounting to 8 mm. Applying the IGS final products, for other cases with reduced duty-cycles and mean satellite numbers, the differences in the RMSE (referenced to the best-case scenario) with and without applying the attitude file are also mostly at several millimeters. This applies also to the cases when using the IGR and IGC products. This suggests that good attitude control can maintain the POD accuracy even when the attitudes are not properly determined or transferred back to the GPC.
Benefiting from the fact that the angle differences between the actual and the expected antenna orientation are small, even without having the attitude information available, differences less than 10 mm can mostly be observed in the final orbits provided that the observation data is complete and continuous (see Figure 12), i.e., with a 3D RMSE amounting to 8 mm. Applying the IGS final products, for other cases with reduced duty-cycles and mean satellite numbers, the differences in the RMSE (referenced to the best-case scenario) with and without applying the attitude file are also mostly at several millimeters. This applies also to the cases when using the IGR and IGC products. This suggests that good attitude control can maintain the POD accuracy even when the attitudes are not properly determined or transferred back to the GPC.

Length of Arc
The orbital period of a LEO satellite is typically around 1.5 h, depending on its altitude from the Earth's surface. Taking the GRACE FO-1 as an example, based on the semi-major axis estimated in the final reference orbit, its orbital period is approximately 1 h 34 min 29 s for the test day. A 24 h arc for data analysis thus contains about roughly 15.2 orbital periods. By reducing the length of the processed arc in case of rapid processing, the amount of the data used for estimation of the orbital parameters may cover, e.g., only several orbital cycles. For an observation model with a good strength, e.g., having data with a high sampling rate, a high duty-cycle, a relatively loose latency requirement, and a high number of observed satellites, the influence of the arc length could be small. However, a long arc might be essential for scenarios with weaker model strengths. To analyze this correlation, this section presents results of a varying arc length from 24 h to 12 h and 6 h. The influences of the arc length on different scenarios are investigated in detail.
Using the best scenario with a duty-cycle of 100%, a mean satellite number of 9, and an observation sampling interval of 10 s, the 3D RMSE of the final orbits are given in Table 5 for different arc lengths and using different IGS products. For arcs shorter than 24 h, multiple arcs were processed covering the entire test day.

Length of Arc
The orbital period of a LEO satellite is typically around 1.5 h, depending on its altitude from the Earth's surface. Taking the GRACE FO-1 as an example, based on the semi-major axis estimated in the final reference orbit, its orbital period is approximately 1 h 34 min 29 s for the test day. A 24 h arc for data analysis thus contains about roughly 15.2 orbital periods. By reducing the length of the processed arc in case of rapid processing, the amount of the data used for estimation of the orbital parameters may cover, e.g., only several orbital cycles. For an observation model with a good strength, e.g., having data with a high sampling rate, a high duty-cycle, a relatively loose latency requirement, and a high number of observed satellites, the influence of the arc length could be small. However, a long arc might be essential for scenarios with weaker model strengths. To analyze this correlation, this section presents results of a varying arc length from 24 h to 12 h and 6 h. The influences of the arc length on different scenarios are investigated in detail.
Using the best scenario with a duty-cycle of 100%, a mean satellite number of 9, and an observation sampling interval of 10 s, the 3D RMSE of the final orbits are given in Table 5 for different arc lengths and using different IGS products. For arcs shorter than 24 h, multiple arcs were processed covering the entire test day. Table 5. 3D RMSE of the phase-based reduced-dynamic orbits with different arc lengths and using different IGS products. 10 s observation data with a duty-cycle of 100% and a mean satellite number of 9 were used for the processing. From Table 5 it can be observed that in case of good data condition, the accuracy differences are very small when reducing the arc length. Having a closer look at the 3D orbital errors with an arc length of 6 h and 12 h using, e.g., the IGS final products (see Figure 13a), border effects of the short arcs, i.e., the degradation of the orbital accuracy at the edges of the corresponding arcs, are the main reason of the degradation [11]. Applying the IGS rapid products (see Figure 13b), the degradation at borders is less obvious due to the generally much larger orbital errors.   Table 5 it can be observed that in case of good data condition, the accuracy differences are very small when reducing the arc length. Having a closer look at the 3D orbital errors with an arc length of 6 h and 12 h using, e.g., the IGS final products (see Figure 13a), border effects of the short arcs, i.e., the degradation of the orbital accuracy at the edges of the corresponding arcs, are the main reason of the degradation [11]. Applying the IGS rapid products (see Figure 13b), the degradation at borders is less obvious due to the generally much larger orbital errors.
(a) (b) Figure 13. 3D orbital errors computed with different arc lengths using (a) the IGS final products and (b) the IGS rapid products. A duty-cycle of 100%, an observation interval of 10 s, and a mean satellite number of 9 were assumed for the processing. The arc length experiment here refers to the processing using e.g., hours of data in each processing round, and in each hour of these hours, data are tracked in % of the time.
For cases with lower data availability, the degradation caused by shortening the arcs could be higher. For scenarios simulating different data conditions, the increase in the 3D RMSE when shortening the arc length from 24 h to 12 h is given in Table 6 as an example applying the IGS final products. While most of the increase is still at millimeters, increase of several centimeters is also observed for the cases with very limited data.  Figure 14 shows the 3D RMSE when processing with the arc lengths of 24 h, 12 h, Figure 13. 3D orbital errors computed with different arc lengths using (a) the IGS final products and (b) the IGS rapid products. A duty-cycle of 100%, an observation interval of 10 s, and a mean satellite number of 9 were assumed for the processing. The arc length experiment here refers to the processing using e.g., k hours of data in each processing round, and in each hour of these k hours, data are tracked in X% of the time.
For cases with lower data availability, the degradation caused by shortening the arcs could be higher. For scenarios simulating different data conditions, the increase in the 3D RMSE when shortening the arc length from 24 h to 12 h is given in Table 6 as an example applying the IGS final products. While most of the increase is still at millimeters, increase of several centimeters is also observed for the cases with very limited data. Table 6. Increase in the 3D RMSE (in cm) of phase-based reduced-dynamic orbits when shortening the arc length from 24 h to 12 h. The IGS final products were used. When shortening the arcs to 6 h and applying worse IGS products, a larger increase in the RMSE can be observed. Figure 14 shows the 3D RMSE when processing with the arc lengths of 24 h, 12 h, and 6 h applying different IGS products. The x-intervals between each two vertical dashed lines represent different cases having the same sampling interval, which are sorted from the highest duty-cycle of 100% to the lowest duty-cycle of 20% first for the mean satellite number of 9, and then similarly for the mean satellite numbers of 7 and 6, respectively. The peaks observed in Figure 14 are thus mostly caused by the low duty-cycles of 40% and 20%. From Figure 14 it can be observed that when shortening the arcs from 24 h (green) to 12 h (blue) or 6 h (red), the increase in the 3D RMSE could be dramatic under bad data conditions, especially when using the IGR and IGC products. While this increase is still limited to around 1 dm when shortening the arcs from 24 h to 12 h, cases with increase above 2 dm appear when shortening the arcs to 6 h and when applying the IGR or IGC products. In cases that the duty-cycle is not lower than 40% and the sampling interval is not larger than 60 s, the 3D RMSE is generally within or around 1 dm when applying the IGS final products or the IGC products, even when processing with arcs with a length of 6 h. Due to the large sampling interval of the IGR satellite clocks, reaching a 3D RMSE of this level generally requires a duty-cycle not lower than 60% when processing 6 h arcs applying the IGR products.
Remote Sens. 2020, 12, x FOR PEER REVIEW 17 of 21 and 6 h applying different IGS products. The x-intervals between each two vertical dashed lines represent different cases having the same sampling interval, which are sorted from the highest dutycycle of 100% to the lowest duty-cycle of 20% first for the mean satellite number of 9, and then similarly for the mean satellite numbers of 7 and 6, respectively. The peaks observed in Figure 14 are thus mostly caused by the low duty-cycles of 40% and 20%. From Figure 14 it can be observed that when shortening the arcs from 24 h (green) to 12 h (blue) or 6 h (red), the increase in the 3D RMSE could be dramatic under bad data conditions, especially when using the IGR and IGC products. While this increase is still limited to around 1 dm when shortening the arcs from 24 h to 12 h, cases with increase above 2 dm appear when shortening the arcs to 6 h and when applying the IGR or IGC products. In cases that the duty-cycle is not lower than 40% and the sampling interval is not larger than 60 s, the 3D RMSE is generally within or around 1 dm when applying the IGS final products or the IGC products, even when processing with arcs with a length of 6 h. Due to the large sampling interval of the IGR satellite clocks, reaching a 3D RMSE of this level generally requires a duty-cycle not lower than 60% when processing 6 h arcs applying the IGR products.

Conclusions
The CubeSat missions have attempted to determine their precise orbits combining dynamic models and dual-frequency GNSS observations. While such kind of reduced-dynamic orbits can reach high accuracy of centimeters for large LEO satellites, the limited size, power, and capacity of CubeSats always set strong constraints on the availability of data and different resources, thus limiting the achievable accuracy of the determined orbits. As an important step at the planning phase of the CubeSat missions, a proper compromise needs to be carefully considered between these constraints and the orbital accuracy, so that limited resources can be better managed and saved for other tasks of CubeSat missions under the required POD accuracy.
In this contribution, with the focus put on the dual-frequency postprocessed CubeSat POD, a detailed sensitivity analysis was performed between the orbital accuracy and different constraints, i.e., the data continuity, the data completeness, the observation sampling interval, the attitude information, the latency of the required orbital products, as well as the processed arc length. It was found that the accuracy degradation of the orbits could be large when several constraints are overlapped with each other, especially when one of the constraints is related to very low duty-cycle, e.g., below 40%. It was also found that the accuracy influences of the observation sampling interval is related to the sampling interval of the GPS satellite clock products used for the processing. In case that the data is otherwise tracked in good condition after considering the assumed constrains, the 3D orbital accuracy for arc lengths of 6 h to 24 h should generally be within or around 1 dm when applying the IGS final products or products from the IGS RTS, provided that the duty-cycle is not lower than 40% and the observation sampling interval is not larger than 60 s.