Precise orbit determination for BDS satellites

Since the first pair of BeiDou satellites was deployed in 2000, China has made continuous efforts to establish its own independent BeiDou Navigation Satellite System (BDS) to provide the regional radio determination satellite service as well as regional and global radio navigation satellite services, which rely on the high quality of orbit and clock products. This article summarizes the achievements in the precise orbit determination (POD) of BDS satellites in the past decade with the focus on observation and orbit dynamic models. First, the disclosed metadata of BDS satellites is presented and the contribution to BDS POD is addressed. The complete optical properties of the satellite bus as well as solar panels are derived based on the absorbed parameters as well the material properties. Secondly, the status and tracking capabilities of the L-band data from accessible ground networks are presented, while some low earth orbiter satellites with onboard BDS tracking capability are listed. The topological structure and measurement scheme of BDS Inter-Satellite-Link (ISL) data are described. After highlighting the progress on observation models as well as orbit perturbations for BDS, e.g., phase center corrections, satellite attitude, and solar radiation pressure, different POD strategies used for BDS are summarized. In addition, the urgent requirement for error modeling of the ISL data is emphasized based on the analysis of the observation noises, and the incompatible characteristics of orbit and clock derived with L-band and ISL data are illuminated and discussed. The further researches on the improvement of phase center calibration and orbit dynamic models, the refinement of ISL observation models, and the potential contribution of BDS to the estimation of geodetic parameters based on L-band or ISL data are identified. With this, it is promising that BDS can achieve better performance and provides vital contributions to the geodesy and navigation.


Introduction
Global Navigation Satellite System (GNSS) is the essential technology to provide space and time information. China decided to establish its own independent GNSS in 1980 s, named as BeiDou Navigation Satellite System (BDS) (Yang et al., 2017a). Twenty years later, the BeiDou Navigation Satellite Demonstration System (BDS-1) with three GEostationary Orbit (GEO) satellites were successfully deployed to provide the regional Radio Determination Satellite Service (RDSS) in 2003. With the launch of the first Medium Earth Orbit (MEO) satellite in 2007, the global Positioning, Navigation, and Timing (PNT) services have been officially provided since 31 July, 2020.
BDS together with other GNSS systems, e.g., Global Positioning System (GPS), Global Navigation Satellite System (GLONASS), and Galileo navigation satellite system (Galileo), can offer high-accuracy positioning to meet the needs of scientific and engineering applications. Furthermore, it helps reduce the spurious signals in the GPS-derived geodetic parameters, e.g., Earth Orientation Parameters (EOP) as well as geocenter position, and improve the realization of the Terrestrial Reference Frame (TRF), which is one of the core tasks of geodesy. All of these rely on the high quality orbit and clock products, which can be determined accurately with precise modeling of well-distributed measurements and orbit dynamic perturbations. The initial access to the BDS signals was provided by the BDS Experimental Tracking Stations (BETS) established by the GNSS Research Center of Wuhan University (WU) since March 2011 (Shi et al., 2012). Besides, the networks of International GNSS Monitoring and Assessment System (iGMAS) (Jiao et al., 2011) and International GNSS Service (IGS) Multi-GNSS EXperiment (MGEX)  project are also the essential data resources for the BDS analysis. For BDS-3 POD, the ISL measurements can also be used alone or in a combination with L-band data.
With these types of data, precise BDS orbits and clocks can be determined. The absence of accurate observation models, e.g., Phase Center Corrections (PCC) for the receiving and transmitting antenna, was noticed (Dilssner et al., 2014). The special yaw attitude control mode has been revealed for BDS-2 IGSO and MEO Dai et al., 2015) as well as for BDS-3 Dilssner 2017;Wang et al., 2018), and its impacts on the orbit determination were investigated for BDS-2 IGSO and MEO satellites (Wang et al., 2013;Guo et al., 2013). Moreover, the deficiencies of the widely used Extend CODE (Center for orbit determination in Europe) Orbit Model (ECOM) of Solar Radiation Pressure (SRP) (Beutler et al., 1994) for satellites in the Orbit-Normal (ON) mode and with elongated shape were observed (Arnold et al., 2015;Dilssner et al., 2018). Although dramatic improvement of orbit quality was obtained due to the improvement of POD strategy as well as observation and orbit dynamic models, the BDS orbit consistency is still about 2-times worse than that of GPS and Galileo (Sośnica et al., 2020).
The aim of this article is to develop the methodology to improve BDS orbit quality and make the further contribution to the estimation of geodetic parameters. It starts with the current deployment status of BDS constellation ("Status of BDS constellation" section), followed by the disclosed metadata from China Satellite Navigation Office (CSNO) ("Metadata of BDS satellites" section). The evolution and status of the ground tracking networks and Low Earth Orbit (LEO) satellites with BDS tracking capability are presented in "Tracking data" section. The improvements and the remaining deficiencies of observation and orbit models are addressed separately in "Measurement models" and "Dynamic models" sections. The strategy for the orbit and clock determination with ground tracking data, LEO onboard data, and ISL data are presented in "Precise orbit and clock determination" section. The incompatible orbit and clock solutions derived from the ISL and L-band data are also discussed in "Precise orbit and clock determination" section, and the contribution of BDS to the estimation of geodetic parameters are presented in "Contributions to geodetic parameters estimation" section. Section "Summary and conclusion" summarizes the findings and the areas for further study. BDS-2 was initiated with the launch of the first MEO satellite (M01) on 13 April 2007. Up to the preparation of this article, a total of 20 satellites, including 8 GEO satellites, 7 IGSO satellites, and 5 MEO satellites, have been deployed, among which 15 satellites are fully operational. For the rest, the BDS-2 M01 was discarded due to an apparent clock problem , and other 3 GEO satellites and 1 MEO satellite are inactive. All the satellites are manufactured by CAST based on the DFH-3A satellite bus and use the Rubidium Atomic Frequency Standards (RAFS) from China and Europe as the primary and backup clock, respectively.

Status of BDS constellation
For the in-orbit validation, BDS-3 experimental (BDS-3s) constellation consists of 1 IGSO satellite (I2S) and 2 MEO satellites (M1S and M2S) made by CAST as well as 1 IGSO satellite (I1S) and 1 MEO satellite (M3S) satellites manufactured by Shanghai Engineering Center for Microsatellites (SECM) of the China Academy of Science (CAS). It was initiated on 30 March 2015 with the deployment of the first IGSO satellite. With respect to BDS-2, different satellite buses, newly designed signals (B1C, B2a, and B2b), ISL in Ka band, updated RAFSs, and new Passive Hydrogen Masers (PHMs), are used. Due to the failure of navigation signal transmitter, no signals are available for M3S . Currently, the rest BDS-3s satellites are also decommissioned for transmitting signals.
Following the success of BDS-3s, the BDS-3 constellation with the full-orbit capability started to build with the launch of the first satellite pair on 5 December 2017. Currently, 3 GEO satellites, 3 IGSO satellites, and 14 MEO satellites manufactured by CAST are active, whereas 10 MEO satellites made by SECM are also in operation to provide global PNT services. The new technology successfully validated by BDS-3s satellites is also used in BDS-3 missions. Besides the primary Ka-band, the laser ISL is also implemented (Zheng, 2020). For SECM MEO satellites, two PHMs are used as the primary frequency standard, while two improved Chinese RAFS clocks are used as the backup. For CAST MEO satellites launched before 2019, four RAFS clocks are equipped, whereas two PHM clocks are used as the primary frequency standard for CAST GEO satellites and IGSO satellites. However, for the rest CAST MEO satellites, only one PHM clock is carried . Table 1 lists the in-orbit status of BDS constellation as 9 September 2021, and the readers can refer the website of Test and Assessment Research Center (TARC) of CSNO for the latest status of BDS system (http:// www. csno-tarc. cn/ system/ const ellat ion).

Metadata of BDS satellites
Satellite metadata, the definition and description of satellite properties, includes unique identifiers like Satellite Vehicle Number (SVN), time of life, the geometry properties (e.g., attitude, the position of mass center, transmit antenna and laser retroreflector array, the shape as well as the dimensions of satellite panels or bus), and the physical properties (e.g., materials, mass, transmit power of signals, optical and thermal properties). They are vital for the accurate modeling of GNSS satellites, particularly for the non-conservative perturbation modeling.
At the end of 2019, CSNO disclosed the metadata for BDS-2 and BDS-3 satellites (CSNO 2019a). The Phase Center Offsets (PCO) of transmitter antenna in B1, B2, and B3 frequency are archived in IGS Antenna Exchange Format (ATX) and will be discussed in "Phase center correction" subsection. While the rests are published in another file containing the satellite identifiers (SVN, Pseudorandom Noise/PRN, COSPARID), block type, active period, mass, laser retroreflector eccentricities, satellite effective surfaces, and the absorption coefficient of satellite bus as well as solar panels. Besides, the attitude law is described in the document 'Definitions and descriptions of BDS/GNSS satellite parameters for high precision application' (CSNO, 2019b) and will be described in "Yaw attitude" subsection. These disclosed metadata can be used for the precise analysis of BDS data. However, it is not enough for orbit dynamic modeling, particularly the SRP, Earth Radiation Pressure (ERP) and antenna thrust, as the specular and diffuse reflection coefficients as well as the signal transmit power are missing.
For illustration purpose, Fig. 1 shows the artist's impression of BDS-3 GEO, IGSO, CAST MEO and SECM MEO satellites from TARC of CSNO and SECM. Those can be found in Yang et al. (2017a) for BDS-2 GEO, IGSO and MEO satellites. Table 2 lists their dimensions, mass, area of solar panels, as well as the ratio of Z to X panel area. Here, the directions of satellite body reference frame follow the IGS conventions (Montenbruck et al., 2015b). As mentioned earlier, BDS-1 satellites are manufactured based on the DFH-3 satellite bus with a size of 2.2 m × 1.72 m × 2.0 m, while the DFH-3A satellite bus is used for all BDS-2 satellites with the similar dimension as its predecessor. Additionally, BDS-2 GEO satellite carries a large Communication Antenna (CA) folded on the +X surface connecting to the +Z surface for time synchronization and RDSS telecommunications with an area approximately 3.14 m 2 . BDS-3 satellites are produced by different manufactures based on their own satellite buses. All GEO and IGSO satellites are manufactured by CAST based on the DFH-3B platform . Moreover, one can see from Fig. 1 that BDS-3 GEO and IGSO are equipped with two and one hexagonal antennas folded on the ±X and +X surfaces, respectively, while there are two smaller circle antennas on IGSO satellite for the in-orbit validation of high-speed data transformation as well as procession . For CAST MEO satellites, a dedicated T-shaped platform is developed to facilitate the lunch of multiple satellites at same time , and there could be an additional panel on the −X surface flipped up for Search and Rescue (SAR) service. However, the satellite body of BDS-3 CAST MEO satellites published by CSNO is a standard cuboid, but it is a bit far from the real. The more precise dimensions are reported by Duan et al. (2021b) and listed in Table 2.
As a new provider of BDS satellites, SECM produces their own satellite platform. Most of SECM satellites are manufactured based on the same bus (SECM-A) expect for M11 and M12 (SECM-B), which have slightly different shapes, as shown in the CSNO released metadata.
The ratio of Z to X panel areas listed in Table 2 indicates stretched degree of the satellite bus. The closer to 1 this value is, the closer to a cubic the satellite bus will be. Hence, BDS-3 GEO, IGSO, and MEO satellites are much more elongated than their counterparts of BDS-2 satellites. However, the mass of BDS-3 GEO and IGSO is about 1.5 times the weight of BDS-2 satellites, whereas BDS-3 MEO satellites are lighter than that of BDS-2. Hence, the deficiencies of non-gravitational perturbations will be more noticeable for BDS-3 due to the large area-to-mass ratio compared to BDS-2. Moreover, BDS-3 SECM MEO satellites are stretched along different direction as BDS-3 CAST MEO satellites, hence, the characteristics of orbit errors will be different and will be presented and discussed in "Solar radiation pressure (SRP)" subsection. Lack of the reflection and diffusion properties limits the accurate modeling of SRP and ERP. For BDS-2 satellites, the optical properties including materials are reported by Chen et al. (2019) from CAST. In general, the satellite surfaces are covered by three kinds of materials, i.e., Multi-Layer Insulation (MLI) blankets with an absorbed coefficient ( α ) of 0.36 and a specularly reflected coefficient ( ρ ) of 0.0, Optical Solar Reflectors (OSRs) with α = 0.135 and ρ =0.865, and anti-static white paint with α = 0.23 and ρ =0.0. The +X, −X, and −Z surfaces are covered by MLI entirely. The area of OSRs on +Y and −Y surfaces is about 2.223 and 2.18 m 2 , respectively, while the rest are MLIs. The +Z surface consists of anti-static white paint (about 1.039 m 2 ) and MLIs. Table 3 lists the corresponding optical properties for all surfaces, which are slightly different from the CSNO released values for the +Y, −Y, and +Z surfaces with two kinds of materials. Solar panels use the silicon photovoltaic solar cells, and the absorbed, reflected, and diffused properties are 0.72, 0.238, and 0.042, respectively . While the optical values for the CA antenna of BDS-2 GEO are assumed the same as that of −Z to facilitate the SRP modeling.  For BDS-3, the optical properties can be roughly derived by taking BDS-2 satellites as the reference. All BDS-3 satellites use triple junction GaAs (gallium arsenide) solar cells for solar panels, and their absorption property is 0.92 and no diffusion. The satellites are also wrapped up with MLIs, as shown in Fig. 1. It seems that +X and −X surfaces of GEO and IGSO satellites are covered by this material completely, while the +Y and −Y surfaces are made of OSRs. For +Z and −Z surfaces, the CSNO disclosed absorbed property is 0.87, which is larger than that of all materials used by BDS-2 satellites. As no further information is available, we simply assume the specularly reflected coefficient of 0.13. More details on the additional antenna mounted on IGSO and GEO satellites are also missing, hence, their optical properties and dimensions are not presented in Table 3. For CAST MEO satellites, +X and −Z surfaces are covered by MLIs totally, while the +Y and −Y surfaces are made of OSRs. The −X surface in fact consists of two parts: one with an area of 1.11 m 2 is OSR, and the other with an area of 1.76 m 2 has high absorbed property of 0.92. The triple junction GaAs used by solar panels also has the same absorbed coefficient. In addition, +Z surface also has the same absorbed property. For these surfaces, no diffusion is assumed. Table 3 lists the average optical properties. Hence, the absorbed coefficient is different from CSNO released values for −X surface. Besides, the details on the SAR antenna are also missing. For SECM satellites, it is reported that all the six surfaces have the same absorption coefficient of 0.20 by CSNO. We assume the OSRs is used, hence, a specular reflection coefficient is 0.80.
It is worth to mention that the optical properties listed in Table 3 are quite coarse, particularly for BDS-3 SECM MEO satellites. However, with the real tracking data, they can be calibrated as shown by Duan et al. (2019Duan et al. ( , 2021b.  ) and Australian Regional GNSS Network (ARGN) of Geoscience Australia also provide BDS observations. As the both contribute to IGS network, we do not discuss them here. The BETS has been operated since March 2011 by GNSS Research Center of WU to monitor the PNT performance of BDS-2 system. This network comprises 15 stations, among which 6 stations are located outside China (Shi et al., 2012). All the stations are equipped initially with UB240-CORS receiver and UA240 high gain antenna manufactured by Beijing Unicore Communications Incorporation, and can track the GPS L1 and L2 as well as BDS B1I and B2I signals. Later some receivers were replaced by Trimble NetR9 to support Multi-GNSS research . Currently, this network is decommissioned. However, it provided the initial tests of BDS-2 signals, and made BDS-2 constellation recognized in the GNSS community, particularly in the early age. The data can be provided by GNSS Research Center on request.
The iGMAS network is built by the iGMAS project aiming at monitoring and assessing the performance of Multi-GNSS system (Jiao et al., 2011). It consists of 31 globally distributed stations equipped with six types of receivers from three manufactures. All receivers can track BDS-2 B1I, B2I, and B3I signals. Currently all of them are also able to track the B1I, B3I, B1C, and B2a open service signals of BDS-3 satellites. Furthermore, B2b signal can be tracked by two kinds of receivers. These data can be downloaded by authorized parties from three data archive centers.
IGS initiated the MGEX project in the middle of 2011, and simultaneously built a network with Multi-GNSS tracking capability as its backbone. The IGS MGEX stations utilize diverse types of receivers mainly from Javad, Leica, NovAtel, Septentrio, and Trimble. In 2016, all MGEX stations were fully incorporated into the official IGS network . As September 2021, there are more than 250 Multi-GNSS stations tracking BDS signals, as shown in Fig. 2. The observations from the IGS Multi-GNSS network are freely available to users from IGS data centers. Figure 3 shows the number of the iGMAS and MGEX stations which track BDS-2 and BDS-3 satellites since 2014. In general, the tracking capability of iGMAS stations is promoted gradually. For MGEX, the tracking capability for BDS-2 is also improved gradually, while a rapid change around the middle of 2018 is observed for BDS-3 due to the deployment of Trimble Alloy and the update of Trimble NetR9, Septentrio POLARX5/ POLARX5TR.
However, not all deployed BDS-3 satellites can be tracked by each receiver, possibly due to the limitation of tracking channels or firmware. Figure 4 shows the daily number of IGS and iGMAS stations tracking the representative satellites, i.e., C19, C23, C32, C35, C38, C41, C56 and C59 from 2019 to 2021. The number of stations with C19 tracked represents those for tracking C19-C22 satellites, whereas that for C23 indicates those for satellites C23-C31, and the similar for the rest. In general, the early launched satellites can be tracked by more receivers. At the beginning of 2019, except for Javad TRE-3 Delta receiver and iGMAS receivers capable to track all deployed BDS-3 satellites, Trimble receivers can only track the C19-C30 satellites, and Septentrio receivers with the latest firmware can only track C19-C22 and C32-C34 satellites. Hence, the number of the stations with C19 tracked is the greatest, then followed by C32.
With the update of Trimble NetR9 as well as deployment of Trimble Alloy, the number of the stations tracking C23 increase gradually, and even exceeds that for C32. Thanks to the update of Septentrio firmware, C35 is tracked by more stations, and the number is over 100 by the end of 2019. Moreover, insignificant difference in the number of tracking stations can be observed for IGSO (C38-C40) and MEO satellites with PRN beyond 40 in 2019, as only iGMAS stations can track them. However, the situation changed at the beginning of 2020. Significant improvement can be observed during January 2021 due to the firmware updates of Septentrio receivers. The different tracking capability will impact the derived orbit quality.

LEO satellites
Besides ground tracking stations, there are a few LEO satellites equipped with GNSS receivers with BDS tracking capability, e.g., LING QIAO (Chen et al., 2016), FengYun-3C/3D (Li et al., 2020c), Luojia-1A (Wang et al., 2020d), and Tianping-1B . The LING QIAO satellite was launched on 4 September 2014, which is the China's first LEO mobile communication experimental satellite and carries two single-frequency GPS receivers for navigation and time synchronization. Besides, one independent single-frequency BDS receiver is also equipped for the in-orbit validation of its performance. However, the actual number of tracked BDS satellites is lower than the expected due to lower acquisition efficiency. Lack of dual-frequency observations limits the enhanced BDS POD with its onboard data.
The above problem is overcome by FengYun-3C satellite, which was launched on 23 September 2013, and developed by the Meteorological Administration/ National Satellite Meteorological Center (CMA/NSMC) of China. This satellite carries a GNSS Occultation Sounder (GNOS) developed by the Center for Space Science and Applied Research (CSSAR) of CAS for positioning, and GPS L1/L2 and BDS B1I/B2I signals from up to six BDS satellites and more than eight GPS satellites can be tracked simultaneously. Zhao et al. (2017) reported that the accuracy is about 0.69, 0.62, 0.38 and 0.820 m for BDS B1I, B2I, GPS C/A, and P2 code measurements, respectively, and is around 2 mm for phase observations. However, the quality of observation drops gradually from 2013 to 2017 due to the aging of the GNOS receiver (Li et al., 2020d). With this data, the elevation-dependent group delay variations of BDS-2 satellites were analyzed, and the BDS orbits, particularly for GEO, are improved due to the rapid variation of the observation geometry . As the successor of FengYun-3C, FengYun-3D carries an enhanced GNOS instrument with up to 17 and 12 tracking channels for GPS and BDS, respectively, resulting in much more observations than that of FengYun-3C. The same dual-frequency signals are tracked, and the loss of tracking GPS L2 and BDS B2I signals, especially at low elevations, has been improved. Hence, FengYun-3D is more stable and has continuous tracking capability than that of FengYun-3C. Thanks to the good performance, its onboard data combined with the FengYun-3C data are also used to enhance BDS POD (Li et al., 2020c).
Luojia-1A is a nanosatellite launched in June 2018 for the technology demonstration of LEO-based navigation signal augmentation as well as nigh-light remote sensing. This satellite is operated by WU and carries an onboard receiver with GPS/BDS dual-frequency tracking capability. However, due to the distribution and the number of BDS-2 satellites, the percentage of the epochs with more than four satellites tracked are around 57%. However, Wang et al. (2020b) shows that the accuracy of code and phase measurements is a few decimeters and millimeters, respectively. Hence, it has potential for BDS enhanced POD.
Tianping-1B is also a nanosatellite equipped with a Multi-GNSS receiver manufactured by Beijing Hexie Avionics Technology Co., Ltd. with 12 tracking channels for GPS and BDS separately . Different from the receivers in the above LEO satellites, its receiver supports BDS-3 tracking and collects code and phase measurements on B1I and B3I frequencies. However, the signals from the BDS-3 satellites with PRN beyond 32 cannot be tracked. In average 8 and 10 satellites can be tracked for BDS and GPS, respectively. The accuracy of a few decimeters is achieved for code measurements. It is vital to analyze BDS-3 observations, but unfortunately the data cannot be freely accessed. In "With ground and onboard L-band data" subsection, the contribution of LEO onboard data to BDS POD will be presented and discussed.

Inter-satellite link data
As one of the significant developments of BDS-3, ISL technology is employed for satellite communication as well as ranging to achieve global SAR service, POD, and time synchronization without the limitation of regional distribution of monitoring stations. Different from GPS BLOCK IIR satellites using Ultra-High Frequency band, the Ka band and laser ISL are used by BDS-3. The Concurrent Spatial Time Division (CSTD) system is employed by BDS-3 ISL to ensure the safe and complete connection among the satellites (Yang et al., 2017b). The word 'concurrent' means that many links are formed simultaneously between the satellites within the constellation, however, there is only one for each satellite. Time-division technology allows each satellite to be connected to other satellites in deferent timeslots (3 s), while the satellite within the view of nadir region (− 60°, 60°) can be targeted using a phased-array antenna, which can switch the beam direction by changing the feed phase of different array elements to implement spatial division. Moreover, the quick orientation of the phased-array antenna ensures the implement of dual one-way observation within 3 s. For ISL ranging, in each timeslot a pair of satellites send signals to each other in turn based on a predefined timeslot schedule, which arranges the order of connection among satellites and is updated weekly by the ground control segment. Within 3 s, the forward and backward links are established in each 1.5 s. Hence, a satellite can establish the links with up to 20 satellites or anchor stations within 60 s, and the links with the whole constellation can be completed in a relative short period to support POD and autonomous navigation (Zheng, 2020). Figure 5 demonstrates the connection for C32 (M13, Slot-B1) with other BDS-3 MEO satellites. The red solid lines represent the continuous satellite links, while the blue dash lines are the discontinuous links. All MEO satellites are evenly distributed on three orbital planes. The satellite cannot connect with its adjacent satellite and the farthest satellite in the same plane due to the nadir angle above 60° as well as earth occlusion, while continuous links with the second and third nearest satellites in the same orbit plane can be formed. For the cross-orbitplane links, there are four continuous and 12 discontinuous links. In addition, each IGSO and GEO satellite can establish discontinuous link with MEO satellites, while the IGSO satellites can be connected with each other.

Yaw attitude
The attitude of a GNSS satellite determines its orientation in space. It is essential for POD, as it has impact on the correction of observation errors and the modeling of non-gravitational perturbations. To ensure that users on Earth can effectively receive satellite signals, the transmitting antenna needs to point to the center of the Earth. Besides, the solar panels should be perpendicular to the direction of the Sun to ensure that the satellite has enough power supply. Under these two constrains, the nominal GNSS satellite attitude is determined. Usually, the Z-axis points to the Earth center, the Y-axis is the rotation axis of the solar panels and perpendicular to the direction of satellite to the Sun, and the X-axis completes the right-hand coordinate system and points to or away the velocity direction (Bar-Sever, 1996;Montenbruck et al., 2015b). However, the satellite does not always follow the nominal attitude law. The yaw maneuvers occur near the midnight and noon points of the orbital plane in low β-angle (the elevation angle of the Sun above the orbital plane) regime, as the required yaw rate exceeds the maximum value provided by the momentum wheel of satellite attitude control subsystem. This is termed as the midnight-/noon-turn maneuvers. Moreover, the eclipse maneuver happens when satellites cross the eclipse seasons, but only for GPS BLOCK IIA, due to the zero output of Sun sensors in the shadow of the Earth ( Bar-Sever, 1996).
Usually, the attitude of different satellite blocks behaviors differently in the maneuver periods. The GEO satellites of BDS constellation adopts ON mode. In this case, the Z-axis points to the radial direction and the X-axis ig. 5 The connection for C32 (M13, Slot-B1) between BDS-3 MEO satellites. The red solid line represents the continuous satellite links, while the blue dash line is the discontinuous links points toward the along-track direction, which thus yields a zero-yaw angle. The Y-axis is perpendicular to the orbital plane and completes the right-hand frame. For BDS-2 IGSO and MEO satellites, the Yaw-Steering (YS) and ON modes are adopted. The transformation from YS to ON occurs when the absolute value of β angle is less than 4° and the absolute value of yaw angle ( ψ angle) is no greater than 5° , while the conditions for switching from ON to YS are |β| > 4 • and |ψ| < 5 • . They present the yaw angles during the attitude switch period from the telemetry data for C08, C11 and C12 satellites, and the sudden rump-up and rump-down jumps instead of smooth switch can be clearly observed at the attitude switch points. With the Revise Kinematic Precise Point Positioning (RKPPP) approach proposed by Dilssner et al. (2011), the yaw angles were estimated with ground L-band measurements Dai et al., 2015). Almost the same conditions as  are derived by Dai et al. (2015). However, by assuming that the attitude mode will switch when the yaw angle is closest to the required orientation, Guo et al. (2017) derived the conditions for yaw attitude switch when µ = 90 • and |β| closest to 4°, where µ is the orbital angle. The above models are quite consistent because the yaw angles at orbit angle of 90° are about ± 4° for the low β angle. However, the analysis of all historical attitude switch events of BDS-2 IGSO and MEO satellites from 2016 to 2019 demonstrates that there are a few events occurring with |β| > 4° for YS to ON and |β| < 4° for ON to YS switch, respectively (Xia et al., 2019).
As the deficiency of widely used ECOM model for SRP modeling in the ON regime, the orbit accuracy degenerates significantly in this period. Hence, the newly launched BDS-2 as well as all BDS-3 IGSO and MEO satellites abandon the ON mode in low β regime, and the continuous YS model is used. In the estimation of the yaw angles, the behavior was first identified for BDS-2 I06, and the corresponding attitude model was proposed (Dilssner 2017;Wang et al., 2018). It is also confirmed that BDS-2 C14 (M06) has abandoned ON mode since September 2017. Generally, the model can predict the yaw behaviors quite well, however, a reverse midnightturn maneuvers was first observed for I06 on DOY 171, 2017, when the β angle was less than 0.1° . Later more evens were identified for BDS-2 C13 and C14 when β angles falls into the range (0°, 0.14°) (Xia et al., 2019). Hence, a yaw bias of 0.14° is introduced to predict the features of the reverse midnight-turn yaw maneuvers (Xia et al., 2019).
The BDS-3 CAST MEO satellites obey the continuous yaw-steering law, while the yaw attitude for SECM MEO satellites was presented by Lin et al. (2018) and is confirmed by CSNO disclosed metadata (CSNO et al. 2019b). Generally, the yaw attitude with β = ± 3° is applied for the satellite when |β| is less than 3°. When the Sun crosses the orbital plane, the sign of the yaw attitude will change in consistency with that of β angle. However, these behaviors could not be validated with L-band data due to near zero values of horizontal PCOs of the satellites. For BDS-3 IGSO satellites manufactured by CAST, the model proposed by Wang et al. (2018) can be used, but still needs further validation. High-rate Satellite Laser Ranging (SLR) and ISL measurements provide the possibility of assessment. With ISL data, the yaw attitudes of BDS-3 IGSO and MEO satellites are derived and analyzed by Yang et al. (2021), and the results clearly confirm the yaw attitude of BDS-3 CAST IGSO follows the model proposed by Wang et al. (2018), however, the yaw behaviors of SECM MEO satellites are occasionally different with the prediction of CSNO disclosed model.

Phase center correction
GNSS observations measure the geometric distance between the antenna phase center of a navigation satellite at signal emission time and the receiving antenna at signal reception time. The phase center is usually different from the Antenna mechanical Reference Points (ARP). Hence, high-precision GNSS applications require a set of PCCs to tie the GNSS measurements consistently to the ARP of the satellite and receiver antenna. Normally, the set of PCCs for the antenna phase center consists of PCO and Phase Center Variation (PCV).
For receiver antenna, the lack of consistent Multi-frequency Multi-GNSS PCCs limits the accurate modeling of GNSS data, particularly for the new systems. The official IGS ATX file contains dual-frequency PCCs of the GPS and GLONASS receiver antennas only. Hence, the values of GPS L1 and L2 were used for BDS data analysis. To overcome the dilemma, the independent PCO/PCV calibrations for BDS and Multi-GNSS signals were conducted and reported (Su et al., 2018;Wang et al., 2020a;Kröger et al., 2021;Krzan et al., 2020;Willi et al., 2020). The PCCs for GPS and BDS signals using the robotic method were calibrated (Su et al. 2018). The similar work was also carried out at CAS (Wang et al., 2020a). However, only two antennas, i.e., TRM5791.00 NONE and TRM59800.00 NONE, were calibrated. Since April 2019, the IGS ACs have carried out the third run of reprocessing the full history of GNSS data in a fully consistent way using the latest models and methodology. The Multi-frequency Multi-GNSS calibration results for a set of antennas were provided by Geo++ for incorporating Galileo into analysis (Wübbena et al., 2019), while the PCCs for BDS open-service signals are also provided for some types of antennas, and they can be used for PCC calibration of the transmitter antenna.
For the satellite antenna, conventional PCOs were initially adopted and recommended by IGS MGEX for BDS-2 based on the approximate spacecraft body dimensions. Afterwards, the PCCs of satellite antennas were calibrated in the IGS08 frame (Dilssner et al., 2014;Guo et al., 2016a;Huang et al., 2018). However, their estimates show significant differences in the Z-offset, especially for IGSO satellites (at meter level), possibly caused by data processing methods and software. Furthermore, the FengYun-3C onboard data were used to extend the PCV estimation of BDS-2 IGSO and MEO (Qu 2021). For the BDS-3 satellites, although the ground calibrated PCO values have already been provided by the CSNO, the preliminary estimates of satellite antenna PCCs in the ITRF14 frame are carried out by Yan et al. (2019a, b), Springer et al. (2020), and Xia et al. (2020). The results show that the discrepancy between the estimated and the ground calibrated PCOs is within a decimeter, and the estimation is affected by the receiver antenna models as well as the SRP model. For the first two estimates, the antenna corrections for BDS B1I and B3I signals are replaced by that of GPS L1 and L2. For the last one, although the receiver antenna PCCs are modeled, only PCOs are calibrated within ITRF14 frame without regards of PCCs. Qu et al. (2021) estimated PCOs and PCVs for BDS-3 IGSO and MEO satellites in ITRF14 and IGS Repro3 frame. Villiger et al. (2021) also calibrated the PCCs for BDS in the IGS Repro3 frame. The two estimates show a quite good consistency and confirm that BDS-3 MEO satellites have an offset of about 10 cm for Z-PCO in IGS Repro3 frame with respect to the CSNC disclosed values. Excellent consistency results are also obtained for the estimates in ITRF14 frame. The PCCs calibration of BDS transmitter antennas facilitates the consistency analysis of Multi-GNSS data for coordinates estimation as well as TRF establishment. However, the current studies only provide the estimation for satellites with PRN less than 37, hence, the PCCs should be further calibrated for the rest BDS-3 IGSO and MEO satellites. Besides, precise point positioning with multi-frequency signals show noticeable improvement in convergence time, precise modeling the measurements also require accurate multi-frequency PCC for satellites antenna, and the initial work has been done by Qu (2021) to calibrate the PCC for multi-frequency signals of BDS.

Biases of L-band data
For BDS, the well-known bias is the elevation-dependent group delay variations. It was diagnosed in BDS-2 satellites first by Wanninger and Beer (2015), and the elevation-dependent piece-wise linear models are proposed to reduce its impacts on BDS data analysis. Because of almost static observation geometry for GEO, the elevation-dependent group delay variations cannot be established with ground data. With FengYun-3C onboard data, the similar biases are identified for IGSO and MEO with the elevation angles above 40°, while larger discrepancy is observed for the elevation angles below 40° possibly due to that the measurements are smoothed . This kind of bias is significantly reduced for BDS-3s and BDS-3 satellites , but still show a variation of 0.1 m revealed by the measurements from a 40-m dish antenna .
Another bias is the inter-system bias between BDS-2 and BDS-3 signals. By analyzing the measurements of a baseline with non-identical receiver pairs (Trimble Alloy and Septentrio POLARX5), the noticeable inter-system biases up to meter level are identified for B1I and B2b/ B2I between BDS-2 and BDS-3 (Mi et al., 2021). This type of bias will cause integer wide-lane cycle difference in the MW linear combination if it is not considered. One wide-lane cycle difference has an impact of about 5 cm on the fractional part of the narrow-lane ambiguity. Therefore, narrow-lane ambiguities can be easily fixed to wrong integers, which will impact the derived orbits and clocks. Figure 6 shows the inter-system biases for the B1I and B3I ionosphere-free combination between BDS-2 and BDS-3 constellation for some stations with different types of receivers. It is clear that receiver dependent biases can be identified, and it should be noted that these values are relative instead of the absolute. For Javad receivers, the bias can be up to 5 m, while up to 3 m for Leica and Trimble receivers. Among them, the Septentrio receivers have the smallest one. Hence, we should treat BDS-2 and BDS-3 as two constellations in data processing, and ambiguity resolution as well as clock estimation must be carefully considered.

Systematic errors in ISL data
For the one-way ISL measurement, the observation corrections mainly include the orbit errors, clock offsets, the eccentrics of ISL antenna, and the hardware delays of the receiving and sending satellites. Besides, there are relativistic effect, ionosphere delay, gravitational time delay, as well as the troposphere delay in the case of the connection with an anchor station. Usually, the clock-and orbit-free observables are formed by the dual one-way measurements for orbit and clock determination .
The orbit-free observables mainly contain the clock offsets as well as hardware delays of the two satellites. Usually, the satellite' hardware delays are assumed constant, and the atomic clock varies linearly, especially in a short period. Hence, the residuals of orbit-free observables with the bias and linear trend removed are used as an indicator of ISL measurement noises. By linearly fitting the 1-h orbit-free observables, the results show that the ISL measurement noises are mainly dependent on the types of ISL terminals (Xie et al. 2019). The ISL noise for CAST satellites is about 1 cm, and around 3 cm for SECM satellites . However, once the orbit-free observables are detrended in long period, e.g., 24 h, the different characteristics of the in-plane links and out-of-plane links can be observed due to the appearance of other slow-varying errors .
As the ISLs measured in different epochs are reduced to the observables at same epoch, the closed-loop observables can be formed by more than three satellites connected with each other. For example, by summing up the three orbit-free observables, the clock offset and hardware delays can be canceled, leaving observation noises only . Figure 7 shows the misclosures for four selected groups. For group of C20, C21, and C41, there are only noises, and no biases and other periodic signals are observed. However, other three groups show noticeable biases and periodic variations. The sources of these errors are still unclear, and they may be caused by the time-dependent ISL hardware delays.
Furthermore, Xie et al. (2020) identifies that the post-fit residuals of ISL clock-free observables for some satellite pairs have constant biases, as shown in Fig. 8  for ISL links. With these biases corrected, the post-fit residuals are noticeably reduced, and the orbit quality is improved further. However, the reason is still unknown, possibly due to the ISL antenna, whose link or channel dependent biases exist. It is also worth to note that these biases can be compensated if the link dependent biases are estimated. Hence, the characteristics of ISL measurements require further investigation to find the sources of errors and the strategy to reduce them. Without doubt, the orbit and clock accuracy will be improved further, particularly for the clock estimated with orbit-free observables as its estimation completely depends on the accuracy of the observations. The reduction or elimination of the systematic errors in ISL observables with better observation models can further refine the stochastic model for orbit and clock determination with ISL only or combined with other measurements, e.g., L-band and SLR data.

Solar radiation pressure (SRP)
As the largest non-gravitational perturbation acting on the GNSS satellite, SRP must be accurately modelled to determine the quality GNSS orbits. Various methods for SRP modeling have been proposed since 1990s, and they can be classified into three types: (1) Empirical models, e.g., the empirical CODE orbit model and its reduced or extend version (Beutler et al., 1994;Springer et al., 1999;Arnold et al., 2015). These models best fit the real GNSS tracking data, though they do not consider the actual physical forces acting on the satellite. (2) Analytical models based on the optical and geometrical properties of the satellite, e.g., Rock models (Fliegel et al., 1992;Fliegel and Gallini 1996) and University College London model (Ziebart & Dare, 2001). The main disadvantage of these models is that they cannot compensate accurately for the real on-orbit behavior of the satellites, e.g., due to the change or uncertainty in the a priori properties of the satellite surface or deviations from nominal attitude (Rodríguez-Solano et al., 2012a). (3) Semi-analytical and semi-empirical models, e.g., the adjustable box-wing model (ABW; Rodríguez-Solano et al., 2012a) and GPS solar pressure model (GSPM; Bar-Sever & Kuang, 2004, 2005. Such models represent intermediate approaches between analytical SRP models and empirical ones and combine a good fit to real tracking data with a clear physical understanding of SRP. Among these SRP models, the reduced version of ECOM (Springer et al., 1999; labeled as ECOM1 in this study) is widely used in the GNSS community. However, it has deficiencies for SRP modeling in ON mode as well as for non-cubic satellites. Hence, the extended ECOM model (ECOM2) are proposed, or it has been augmented with a prior model.
For modeling the SRP for BDS satellites, many efforts have been made. For GEO, POD is particularly challenging due to the almost static geometry as well as the ON attitude mode employed. As a result, the strong correlations among orbital elements, SRP parameters, and ambiguities are present. To cope with these correlations, Steigenberger et al. (2013) proposed the estimation of only the D 0 of ECOM1. However, Liu et al. (2016) demonstrates that the GEO orbit can be improved by estimating six parameters (D 0 , Y 0 , B 0 , D s , B s , and Y c ) instead of the typical five parameters (D 0 , Y 0 , B 0 , B c , and B s ) of ECOM1 model. Tan et al. (2016) developed an analytical model of SRP based on a ray-tracing approach to use as the a priori model for the ECOM1. Wang et al. (2019a) identified that the perturbation caused by the C-band communication antenna can generate the Sun-elongation-angle-dependent variation and the bias of about 14.9 cm in BDS-2 G01 SLR residuals. Besides, the ON attitude mode used by BDS GEO satellites as well as an orbital inclination of nearly 0° result in strong linear correlations between the satellite's initial position in the Z-axis and the Y 0 parameter. An empirical SRP model is established for BDS-2 GEO satellites to enhance the ECOM1 using an empirical fitting approach. Better than 10 cm Root-Mean-Square (RMS) of SLR residuals is achieved, and an improvement by 4-5 times over the ECOM1 model is obtained. However, the reduced performance is observed for WU's MGEX C01 solution after the replacement of G03 by G08, hence, the a priori model should be re-estimated.
For BDS-2 IGSO and MEO satellites, the challenge is to determine the accurate orbits in ON mode. Dramatic degeneration of orbit accuracy is first observed for BDS-2 IGSO and MEO in ON mode (Wang et al., 2013), and the  Guo et al. (2013), Prange et al. (2016) also confirm that ECOM model are strictly designed for the satellites in YS mode. This attitude-related POD issue was analyzed by , and the shortcomings of the ECOM1 can be compensated by adding constrained empirical accelerations in the along-track direction, resulting in significant improvement of the orbit accuracy Guo et al., 2016a). Prange et al. (2020a) also developed a series of ECOM-like empirical SRP models for the satellites in ON mode, which are slightly different for BDS-2 IGSO and MEO. The SLR validation confirms that the orbit accuracy of BDS-2 IGSO and MEO satellites can be improved by a factor of 2 in ON mode with all above models, but still lower than those in YS mode. Furthermore, semi-analytical SRP modeling for QZS-1 satellite in ON attitude is developed by Montenbruck et al. (2017b) based on a generic box-wing model, and it can also be used for modeling SRP of BDS-2 IGSO and MEO satellites in ON mode. However, few studies were conducted to demonstrate its performance for BDS satellite, and further validation is needed. To ensure a proper use of the a priori box-wing model, the optical parameters of BDS-2 satellites are estimated (Duan et al., 2019), and improvement by more than 60% over the initial values was obtained.
The ECOM1 has also been used with the emergence of BDS-3 MEO satellites. However, the SLR validations show obvious linear systematic errors with respect to the Sunelongation (ε) angle for MEO satellites (shown in Fig. 9), and its patterns are similar to that of Galileo IOV and FOC satellites (Montenbruck et al., 2015a;Sośnica et al., 2020). CAST and SECM MEO satellites have opposite patterns of SLR residuals with the similar absolute value of slope due to that the Z-bus area with transmitter antennas is much larger than the X-bus surface area for SECM satellites and opposite for CAST. The systematic patterns of the SLR residuals can be reduced by using the a priori box-wing or cuboid box model based on CSNO released or calibrated metadata of BDS-3 satellites (Li et al., 2020a;Duan et al., 2021a;Yan et al., 2019b), as well as the purely empirical ECOM2 model (Yan et al. 2019b). In general, better performance can be obtained for those models with calibrated optical coefficients. Besides, analyzing the estimated B 0 parameters of ECOM1 model, Li et al. (2021a) identifies that C45 (M23) and C46 (M24) satellites show different variations compared with other CAST satellites, possible due to different structure. Besides, two SECM satellites (C34/M15 and C35/M16) have relative larger D 0 with respect to other SECM satellites possible due to higher thermal radiation of solar panels or much smaller mass than the published value, as shown in Fig. 10. These mismodeling issues need further investigation.
For BDS-3 GEO and IGSO satellites, it is a challenge to model the SRP perturbations, as the additional antenna and elongated shape should be considered. Lack of accurate SRP model for GEO and IGSO limits the determined orbit accuracy to the decimeter and meter level, respectively, and needs further investigation. However, the approach for the establishment of empirical SRP for BDS-2 GEO can be referenced .

Antenna thrust
Antenna thrust is a perturbation along the radial direction generated by the signal transmission of satellites. The total transmitting power and the in-orbit mass of the satellite are the prerequisites for modeling this perturbation. For BDS, the mass has been disclosed by CSNO, however, the lack of total transmitting power limits the accurate modeling of antenna thrust. Fortunately, the values are measured with a 30-m highgain dash antenna of the German Aerospace Center (DLR) with considering the power loses along the propagation path between satellite and ground receiver by Steigenberger et al. (2018) and Steigenberger and Thoelert (2020) for BDS-2 and BDS-3 satellites, respectively. Generally, 180 W and 130 W are recommended for BDS-2 IGSO and MEO, while 310 W and 280 W are assumed for BDS-3 CAST and SECM MEO satellites. With the modelled antenna thrust, a bias in the orbit radial direction can be observed by SLR. In general, a radial effect is about 28, 5, 16, and 19 mm for BDS-2 IGSO, MEO, BDS-3 SECM MEO, and CAST MEO satellites, respectively . As expected, a larger bias is observed for BDS-2 IGSO than that of MEO, as the former has greater signal transmission power with the similar mass as that of MEO (see Table 1). While the orbits of BDS-3 CAST MEO satellites shift much more than that of SECM ones due to greater transmit power and lighter mass compared with SECM MEO satellites.

Earth radiation pressure
ERP is a non-gravitational perturbation caused by the reflected or re-emitted solar radiation by the Earth, which is mainly along the radial direction. Similar to the antenna thrust, the ERP mainly affects the radial component of the orbit. The detailed metadata of the satellite, including satellite dimensions, optical and infrared coefficient, are useful for modeling ERP. Rodriguez-Solano (2009) presented the model for ERP computation, and it has been applied by IGS ACs for GPS and GLONASS.
With the values listed in Table 3, the impacts of ERP on BDS orbits are evaluated. SLR validation demonstrates that the radial components of orbits are biased by about 25, 20, 15, and 12 mm for BDS-2 IGSO, MEO, BDS-3 CAST MEO, and SECM MEO satellites, respectively .

Thermal radiation pressure
The thermal radiation of a satellite is the electromagnetic radiation emitted from the surface-covered materials of the satellite as well as its entrails and generates a perturbation, called thermal radiation pressure, to the satellite. In this study, the thermal force caused by the re-radiation of the absorbed energy from the Sun in the form of heat is indicated as Thermal Radiation Pressure (TRP), whereas the force due to the emitted heat from the thermal radiator is called Thermal Radiation (TR). Since TRP and SRP show similar characteristics as both are originated from the Sun, the former is often incorporated with SRP in a combined form. However, if the satellite is inside eclipse seasons, thermal radiation should be carefully modeled. It is believed that mismodeling its impacts limits the POD accuracy in eclipse seasons. On the other hand, the heat generated by the satellite internal subsystems is usually Beta ( vented to space via radiators or louvers on the surface of the bus. Hence, the details of the thermal radiators for a satellite are required for TR modeling. However, this kind of information is not disclosed by the spacecraft manufacturer. Usually, the radiators are assumed to be monumented on the unilluminated surfaces, i.e., +Y, −Y, and −X. Recently, Sidorov et al. (2020) identified that TR force due to the emitted heat from the thermal radiator caused by the onboard atomic clock generated the orbit errors in eclipse seasons for Galileo satellites. Moreover, based on the analysis of Y 0 of ECOM1 model, the potential radiators on the −X surface were identified for GLONASS satellites (Duan et al., 2020). Wang et al. (2019c) also introduced a periodic acceleration along +X direction to compensate TRP for BDS satellites, particularly for C13. The systematic errors dependent of Sun elongation angle can be significantly reduced, and SLR validation indicates that the orbit accuracy is improved by 10-30% to approximately 4.5 cm and 3.0 cm for BDS-2 IGSO and MEO satellites, respectively (Wang et al., 2019a, b). Besides, some BDS-3 MEO satellites have clear radiator emissions, as their Y 0 estimates are not zero, as shown in Fig. 11. Comprehensive analysis and modeling of the radiator emission and the thermal radiation of solar panels are presented by Duan et al. (2021b) for all BDS-2 and BDS-3 satellites, and the orbit misclosures of BDS-3 satellites are reduced by a factor of two for the ECOM1 model during eclipse seasons.

Precise orbit and clock determination
In the above two sections, the observation and dynamical models for BDS POD were discussed. This section will focus on the POD strategy, and present the achievable accuracy of determined BDS orbits.

With ground and onboard L-band data.
For BDS, the ground L-band measurements are usually used for orbit determination. At the beginning, due to the limitation of tracking data a three-day solution is usually employed for orbit determination to improve the solution strength. Furthermore, the two-step approach is applied, which uses the GPS data to derive the common parameters, i.e., the station coordinates, epoch-wise receiver clock offsets, and Zenith Troposphere Delay (ZTD) parameter. They are fixed later for BDS POD. This approach benefits from the application of GPS constellation and results in a better quality for the common parameters. With this strategy, the initial POD results for BDS-2 satellites are reported by Shi et al. (2012) and Zhao et al. (2013) as well as Steigenberger et al. (2013) based on 15 BETS stations and 6 CONGO stations, respectively. Furthermore, the BDS-3s orbits are determined (Tan et al., 2017;Zhao et al., 2018). On the other hand, all parameters are estimated simultaneously from BDS only or its combination with GPS measurements, as done by Ge et al. (2012) for BDS-2 and Li et al. (2019) for BDS-3s. With the increase in tracking stations as well as their better distribution, the one-step is widely used for BDS data analysis, as used by most of IGS MGEX AC (Loyer et al., 2018;Selmke et al., 2018;Prange et al., 2020a, b), and better orbit consistency can be obtained. Currently, the BDS-2 GEO orbit consistency between

GFZ (Deutsches GeoForschungsZentrum Potsdam) and
WU is at the 2-4 m level, and the best agreement of about 15 cm is achieved between GFZ and WU among all MGEX AC for IGSO and MEO satellites (Steigenberger & Mentenbruck, 2020). While the consistency in BDS-3 orbits between European Space Operations Centre of the European Space Agency (ESOC/ESA) and WU products is about 12.3, 9.2 and 4.1 cm in along-track, cross-track, and radial direction for CAST MEO satellites, while for the SECM MEO satellites there are larger mean RMS values of 25.0, 15.2 and 5.7 cm in the respective direction . The SLR validation demonstrates that the achieved accuracy is below 20 cm, 5-7 cm, 3.5 cm for BDS-2 GEO, IGSO, and MEO satellites, while the RMS values of 3-4 cm are achieved for BDS-3 MEO satellites using the improved SRP model . To achieve the better quality of orbit solution, integer ambiguity resolution is essential. Usually, the Double-Difference (DD) ambiguities are fixed to the integers. With the development of the integer ambiguity resolution for Un-Difference (UD) ambiguities, it has been employed for GNSS POD. GFZ announced that their MGEX rapid satellite orbit and clock are estimated with the UD Ambiguity Resolution (UD-AR) using daily wide/ narrow-lane Uncalibrated-Phase-Delay (UPD) method. Compared to DD, UD-AR for POD has two advantages: more phase observations can be fixed, and the selection of independent DD ambiguities can be avoided for a set of independent baselines. The 6 h orbit overlapping differences can be improved by about 30% and 18% for BDS-2 and BDS-3, respectively (Deng 2021). For BDS-3, the backup B1I and B3I signals are usually used for POD, as they are the overlapping frequencies supported by both BDS-2 and BDS-3. However, the designed signals for compatibility and interpolation are B1C and B2a, with which Li et al. (2021b) demonstrates better POD results can be obtained due to their low noises.
Besides, the LEO onboard data can also be applied for the BDS POD. The fast movement of LEO can result in a rapid change in observation geometry, and better geographical distribution of observations. Hence, it can be expected that the orbit accuracy can be improved, and the number of ground stations for achieving the required accuracy can be reduced. Zhao et al. (2017) showed that the orbit solution was enhanced by FengYun-3C onboard data. Li et al. (2020c) also presented the integrated POD for GPS, BDS, and FengYun-3C. The most pronounced benefit is observed in BDS GEO orbits, which is improved by 44% for the regional solution and 41% for the global solution. With more LEO onboard data used, the better accuracy is expected.
As the orbit and clock are determined simultaneously using the L-band data, these two types of the estimated parameters are highly correlated. To separate them, the DD approach is used for orbit determination, as CODE does in the IGS MGEX ACs. Besides, L-band Two-Way Satellite Time and Frequency Transfer (TWSTFT) and ISL are applied for time synchronization between BDS satellites and the ground master. Hence, the clocks derived from TSSTFT and ISL can be fixed for BDS POD. This approach has been used to generate the ephemeris of BDS-2 (Zhou et al., 2011;Guo et al., 2015) and BDS-3 Pan et al., 2021) which are detailed by Zhou et al. (2020). Moreover, the satellite clock parameters within a relative short period can be constrained to a linear model, as GNSS onboard frequency standards are quite stable, particularly for the improved RAFSs and PHMs carried by BDS. This idea was assessed by using Galileo and BDS-2 data (Hackel et al., 2015;Qing et al., 2017). The reduction of orbit boundary discontinuities and SLR residuals is clearly observed.
Besides L-band data, the SLR measurements can also be used for the BDS POD, especially for the case of few L-band data available. For example, Hauschild et al. (2012) determined the orbit of BDS-2 M01 using the SLR tracking data. For BDS-2 IGSO and MEO, Bury et al. (2018) analyzed the impacts of orbit arc length and the number of SLR measurements on the orbits. The consistency of MEO orbits derived with SLR-only is 22.9, 13.4, and 4.4 cm in along-tack, cross-track, and radial direction, respectively, while that of BDS-2 IGSO is much worse due to the poor SLR observation geometry. As no clock to be estimated, it is expected that SLR can contribute to the improvement of the GEO orbits with a combination with L-band data (Sun et al., 2016). The bias in the microwave orbit solutions can be reduced by 20 cm using SLR measurements. However, generally, SLR's contribution to the BDS POD is marginal due to the poor distribution as well as fewer measurements.

With inter-satellite-link data.
For BDS-3 s and BDS-3 satellites, ISL technology is used to implement autonomous and global navigation with less or no ground segment supports. In this case, the orbit and clock can be precisely determined even when the satellite moves out of the coverage of ground monitoring stations. Tang et al. (2018) and Yang et al. (2017b) reported the initial POD results for BDS-3 s satellites with ISL observables only. With the deployment of BDS-3 satellites, the POD results with ISL data were presented by many researches, e.g., Wang et al. (2019b), Xie et al. (2019), and Yang et al. (2019). For the ISL-only solution, the 3D orbit consistency is less than 20 and 30 cm for MEO and GEO/IGSO, respectively. Once the ISL data is combined with the ground tracking data, the 3D accuracy is improved to 7 cm, 13 and 20 cm for MEO, IGSO, and GEO, respectively . For the ISL-only clock estimation, although the accuracy of the clock estimated by ISL observations is about 0.2-0.3 ns due to larger noises of ISL, the estimated clocks are free of orbit errors  and show better stability in a long period. The analysis performed by Cai et al. (2020) indicated that the ISL data from different orbit planes contribute to POD more than that from the same plane, as they have better observation geometry.
Instead of processing the clock-free and orbit-free ISL measurements for orbit and clock determination, Ruan et al. (2020) proposed a method to process the raw ISL measurements directly to obtain satellite orbits, clocks, and hardware delays of ISL equipment simultaneously. Compared with the dual-way model, it can process ISL observables as much as possible, and obtain consistent orbit and clock solutions. However, the clock will be potentially contaminated by the orbit errors.

Orbit and clock products
Currently, iGMAS and IGS MGEX provide the BDS orbits. Five ACs contribute BDS-2 and BDS-3 products to IGS MGEX routinely, i.e., CODE, GFZ, Shanghai Astronomical Observatory (SHAO), Information and Analysis Center (IAC), and WU, while 13 ACs to iGMAS. Although WU and SHAO are the ACs for both MGEX and iGMAS, different solutions are submitted to each project. Among MGEX ACs, CODE does not provide BDS GEO solutions. Besides, ESOC/ESA has also been providing BDS-3 orbit and clock solutions since DOY 1, 2019 (http:// navig ation-office. esa. int/ produ cts/ gnssprodu cts/). Within iGMAS, the ultra-rapid, rapid, and final products for BDS are generated from those submitted by all ACs, while only WU releases its hourly ultrarapid Multi-GNSS solutions to support (near) real-time applications. The consistency in the orbit and clock products for each AC with respect to the combined solutions is presented on the iGMAS website (http:// www. igmas. org). For IGS MGEX, the consistency in the MGEX orbit products is assessed by comparing the results from different ACs, whereas their accuracy is evaluated by SLR residuals. Complementary plots of such quality assessments are presented on the MGEX website with weekly updates (see https:// igs. org/ mgex/ analy sis), and summarized by Steigenberger and Montenbruck (2020). Guo et al. (2016b) made a comprehensive assessment of the MGEX BDS orbit products from January 1, 2013, to May 1, 2015. As the successor, Li et al. (2020b) have assessed the products quality since 2018.
As the orbit and clock solutions have been provided by many ACs, it is possible to implement the Multi-GNSS combination. This has been done within iGMAS (Chen et al., 2015). IGS Analysis Center Coordinator initiated an experimental Multi-GNSS orbit combination service in 2019 by adapting the combination software used for many years for IGS GPS and GLONASS combinations. The Multi-GNSS orbits are combined based on the individual products generated by IGS and MGEX ACs. The validation performed by Sośnica et al. (2020) demonstrated the good quality achieved for the combined solutions. However, the BDS-3 MEO satellites suffer the orbit errors for mismodeling SRP. This is expected, as ECOM1 model is used to generate the WU BDS-3 orbits in the selected period. The mean standard deviations of SLR residuals are 87, 51, 40 mm for BDS-2 GEO, IGSO, and MEO satellites, respectively. Besides, WU also release the combined orbits for Multi-GNSS satellites based on the MGEX AC products .

Consistent orbits and clocks from ISL and L band data
Previous studies demonstrated that the SLR residuals of BDS-3 CAST and SECM MEO satellite orbits determined with the ground L-band data had linear trend but opposite dependency with respect to the position of Sun. Usually, this behavior is believed to be caused by the large area-to-mass ratio as well as the elongated shape. With the ISL data only in 2019, the orbits for those satellites are determined based on ECOM1 model. However, the linear variations in the SLR residuals with respect to Sunelongation angle are significantly reduced, as shown in Fig. 12 for C20 and C30. The possible explanation is that the observations have different sensitivity to the satellite position in each direction. The ground L-band data can only capture the variations in radial direction, while ISL data observe the orbit errors in along-track more precisely, resulting in the precise estimation of SRP along D direction. This explanation needs further investigation.
Besides, the different characteristics are also observed for the clocks determined with L-band and ISL data separately, as shown in Xie et al. (2020). As the noise level of Ka-band ISL measurements is about 10 times higher than that of the L-band carrier phase observations, the clocks derived with ISL data display a superposition of random walk and white noises, and the corresponding Allan Derivation (ADEV) with the sampling interval of 2000 s is generally worse than that of the L-band clocks. However, a pronounced "bump" appears at around second 20,000 for ADEV of the L-band clocks, as the clocks are contaminated by the orbit errors, while ISL derived clock do not show such behaviors. The spectra analysis also demonstrates that ISL and L-band clocks have different characteristics. For L-band clocks, there are pronounced 1-CPR (cycle-per-revolution), 2-CPR, and 3-CPR harmonics. However, only 1-CPR and 2-CPR signals are observed for ISL clocks. And the amplitude of 1-CPR is only a half of that of L-band clocks. These indicate the orbit errors are significantly reduced within ISL clocks. Hence, for the clock estimation with the L-band and ISL data, it is not just a simple combination of the measurements, but the different characteristics of both types of data needs considering to get clock products with good stability in both short and long periods.

Contributions to geodetic parameters estimation
In addition to the PNT services, GNSS also play a fundamental role in the establishment of the TRF with the estimation of the EOP, geocenter, and terrestrial station coordinates with orbit and clock simultaneously. These geodetic parameters are time-variable with a geophysical origin. Theoretically, their estimates should be independent of the GNSS constellation. However, the spurious draconitic signals are identified in the GNSS-based estimates originated from orbit modeling issues (especially the SRP) and GNSS constellation characteristics. As shown by Scaramuzza et al. (2018), 3 instead of 6 orbital planes in the constellation may lead to spurious signals in polar motion, especially at the harmonic of 3 cycles per draconitic year (cpy). Hence, the spurious signals in polar motion can be expected for BDS as 3 orbital planes used. And serious collinearity issues due to the simultaneous estimation of epoch-wise station positions and satellite clock offsets and tropospheric parameters in global GNSS data analyses contaminate the determination of all three components of geocenter . To improve the estimation, the collinearity issues can be reduced by the simultaneous analysis of Multi-GNSS data collected by ground stations and LEO satellites, the modelling of ultra-stable satellite docks, and the mitigation of orbit modelling errors . Zajdel et al. (2020Zajdel et al. ( , 2021 investigated the contribution of Multi-GNSS to the estimation of EOP and geocenter with the focus on the SRP modeling and length of orbit arc. However, the contribution of BDS, particularly BDS-3, is less investigated. Xu et al. (2014) investigated the impacts of BDS-2 on the EOP estimation with a combination of GPS. Recently, Duan et al. (2021) analyzed the impacts of BDS SRP models on EOP and geocenter estimation and show that the use of the a priori box-wing model can mitigate a large portion of the spurious signals in the geodetic parameters.
Another issue for the GNSS contribution to TRF is its weak ability to determine the terrestrial scale, as it cannot be separated from conventional satellite PCOs. Two approaches can be used to overcome the issue. One is to fix the transmitter antenna patterns of at least one satellite, and the other is to use LEO onboard data for the PCO estimation because the orbits of the LEOs are scale independent. Although the PCO values are not available for GPS and GLONASS satellites, the ground calibrated PCO values for Galileo, BDS, and QZSS have been released, making a reliable estimation of the terrestrial scale with GNSS possible. The preliminary estimation by Qu et al. (2021) demonstrates that the scale inconsistency derived from BDS and Galileo released PCOs reaches about +1.854 ± 0.191 ppb (part-per-billion). Hence, more efforts should be made to bridge the gap to obtain a consistent scale determined by different GNSS constellations. One of the possible ways is to use the onboard BDS, Galileo, and GPS data from LEO satellites, as done by Huang et al. (2021).
Besides ground and onboard L-band data, BDS also provide ISL measurements, which has potential for geodetic estimation. The simulation performed by Glaser et al. (2020) assess the potential improvements of the Kepler constellation on the TRF origin and scale. The results demonstrate that the fully developed Kepler system significantly improves the geocenter estimates, as it increases the reliability due to a complete de-correlation of the geocenter coordinates and the orbit parameters related to the SRP modeling. The inclusion of ISL can also make the estimation of the Z component of geocenter better by a factor of 13. However, the preliminary validation demonstrates the BDS ISL has marginal impacts on the estimation of geodetic parameters . This is possibly due to larger noises (2-3 cm) of BDS Kaband ISL than that of Kapler laser ISL (1 mm) as well as unmodelled systematic errors in BDS ISL data as shown above. Hence, the ability of BDS for geodetic parameter estimation with ISL should be explored further.

Summary and conclusion
Since the development of BDS, substantial process in the model and strategy for BDS POD have been made. The expanding tracking network and the disclosed metadata lay the vital foundation for the refinement of observation model and dynamic model for BDS POD. Initial orbit and clock solutions have been made publicly available to support precise positioning and to identify the weaknesses of POD strategy. Noticeable improvements in the calibration of the PCCs for receiver and transmitter antennas, yaw attitude, biases, SRP modeling, and antenna thrust have been made. POD with ISL measurements also achieves a great success. The SLR validation shows the accuracy of BDS-3 orbits is at similar level as GPS, however, the consistency of the orbit is still lower than that of GPS. Besides, the further study of new ISL measurements are still needed regarding its characteristics and systematic errors to explore the potential not only in orbit and clock determination, but also in geodetic parameter estimation. Hence, the following works remain to be done to improve BDS POD accuracy.
• Further disclosure of the completed and detailed spacecraft parameters, e.g., the completed dimensions, optical coefficients as well as transmitting power, by the manufactures or system provider. • Calibration of the PCOs and PCVs for all BDS satellites based on the IGS recommended PCCs of receiver antennas using the ground and LEO onboard data. • Further refinement of the SRP perturbation for BDS satellites, particularly for BDS-2 IGSO and MEO in ON mode as well as BDS-3 IGSO and MEO satellites with the consideration of the additional antennas and different characteristics of satellite groups. • Modeling the thermal radiation for better orbit determination and predication in eclipse seasons. • Investigation of the characteristics of ISL measurements for mitigating the observation errors. • Analysis of the contribution of the L-band and ISL data to orbit and clock products, and investigation of the sources of SRP-deduced orbit errors.
• Orbit determination with short-term clock modeling for ultra-stable clock or clock-free measurements. • Full exploration of the BDS contribution to geodetic estimation as well as the establishment of terrestrial reference frame with the focus on the scale and origin.
No doubt, these tasks should be investigated in a close cooperation with the manufacturers and system provider. Better orbit solution for BDS will be obtained with the combination of various kinds of observations.