## 3 LEO-Constellation-Augmented BDS POD Strategies  
For the algorithms of the GNSS/LEO combined POD, readers can refer to [PERSON] et al. [7; 10]; the processing strategies and parameters to be estimated are listed in Table 3. Considering a regional network was used, the processing arc was three days. The dynamic models applied for POD are almost the same as those used for the constellation configuration and are listed in Table 2. Note that there are two differences between the force models for the BDS and LEO POD. The atmospheric drag model was not considered for the BDS because of the thin atmosphere at its orbital altitude. Unlike GNSS satellites, the trajectories of LEO satellites undergo perturbations from both the solar radiation pressure and the atmospheric drag; therefore, the atmospheric drag plays a dominant role in the non-gravitational forces. The air density model DTM94 was applied for the LEO satellites and the parameters were estimated every 6 h. Empirical accelerations were only considered for the LEO satellites; one set of piecewise periodic parameters per revolution were added to the along, cross, and radial components and were estimated every 2 h.  
The detailed POD procedure is as follows. 1. First, the designed BDS and LEO ephemerides are converted to a broadcast ephemeris and used as the initial values for the POD. An 18-parameter model [31] is used for the BDS broadcast ephemeris design. 2. The initial BDS positions are obtained from the generated broadcast ephemeris and then the BDS initial dynamics information is derived via orbit dynamics fitting. The positions of the LEO satellites are derived from the kinematic POD using the BDS broadcast ephemeris and LEO spaceborne data, and then the LEO orbit initial dynamics information is obtained via orbit dynamics fitting. 3. Orbit integration is applied for the BDS and LEO satellite dynamics to obtain the initial reference orbit and state transition matrix. 4. The unknown parameters are ranked, combined, and added with the constraints listed in Table 3, and the POD equations are linearized at the reference initial orbit. The observations at all epochs are stacked on a combined adjustment for parameter estimation. 5. The adjusted BDS and LEO reference orbits are obtained via orbit integration after having resolved the normal equations. 6. Considering several random cycles of cycle slips are added to the observations, posterior residual editing is conducted until new cycle slips are detected and calibrated. 7. Steps 3-5 are repeated until the residuals are smaller than the given threshold values. In terms of the carrier-phase measurements, the threshold is usually defined as three times the post-fit STD value. The carrier-phase measurement can be regarded as \"bad\" and abandoned if the corresponding residual is larger than the threshold value. An ambiguity parameter is introduced and involved in calculation when the residual is larger than one and a half times the STD value. Otherwise, the residuals are regarded as white noises with normal distribution. This step is necessary because the actual-measured code and carrier-phase residuals are included in the observation simulation.