A Priori Solar Radiation Pressure Model for BeiDou-3 MEO Satellites

Due to the cuboid satellite body of BeiDou-3 satellites, the accuracy of their orbit showed a trend of systematic variation with the sun-satellite-earth angle (ε) using the Extend CODE Orbit Model (ECOM1). Therefore, an a priori cuboid box-wing model (named the cuboid model) is necessary to compensate ECOM1. Considering that the body-dimensions and optical properties of the BeiDou-3 satellites used to construct the box-wing model have not yet been fully released, the adjustable box-wing model (ABW) was used for precise orbit determination (POD). The a priori cuboid box-wing model was directly estimated by the precision radiation accelerations, obtained from ABW POD. When using ECOM1 model, for 14 < β < 40◦, a linear systematic variation of D0 related to the elevation of the sun above the orbital plane (β-angle) with a slope of 0.048 nm/s2/◦, was found for C30. After adding the cuboid model to assist ECOM1 (named Cuboid + ECOM1), the slope was reduced to 0.005 nm/s2/◦, and for C20 satellite, the standard deviation (STD) of D0 was improved, from 1.28 to 0.85 nm/s2 (34%). For satellite laser ranging (SLR) validation, when using the ECOM1 model, the systematic variation with the ε angle was about 14 cm for C20 and C30. After using the Cuboid + ECOM1 model, the variation was significantly reduced to about 5 cm. For C20 and C21, compared with the ECOM1 model, the root mean square (RMS) of the ECOM2 and Cuboid + ECOM1 model was improved by about 0.54 (10.3%) and 0.43 cm (8.7%). For C29 and C30, the RMS of ECOM2 and Cuboid + ECOM1 model was improved for about 0.7 (10.9%) and 1.6 cm (25.6%). Finally, the RMS of the SLR residuals of 4.37 to 4.88 cm was achieved for BeiDou-3 POD.


Introduction
The construction of the BeiDou-3 global navigation system began in November 2017. As of April 2019, 18 Medium Earth Orbit (MEO) satellites, one Geostationary Orbit (GEO) satellite, and one Inclined Geosynchronous Orbit (IGSO) satellite were launched into orbit. By the end of 2020, the full operational capability of the global constellation will be completed, and there will be three GEO satellites, three IGSO satellites and 24 MEO satellites [1,2].
Satellite precise orbit determination (POD) is an indispensable task in the global navigation satellite system (GNSS) service. The main factors affecting the accuracy of satellite orbit determination are the distribution of ground tracking stations, data quality of tracking stations, and accuracy of satellite-end models. Satellite attitude models, the phase center model, and the solar pressure model (SRP) are the main models of the satellite-end, and the SRP is the most important. Extending the CODE Orbit model (ECOM1) with five parameters (D0, Y0, B0, Bc, and Bs) produces an empirical SRP model, which is simple and does not require an a priori model. It was widely adopted for the SRP of GPS and GLONASS for the most of Multi-GNSS Experiment (MGEX) analysis centers [3][4][5]. However, for the

ECOM1 SRP Model
Solar radiation pressure, which is the main source of non-gravitational perturbation of the GNSS satellites [22], and has significant influence on precision orbit determination, depends on the optical properties of the reflecting surface and its orientation with respect to the incident radiation.
In the absence of optical properties and dimensions, as an empirical SRP model, the ECOM1 is widely used among most International GNSS Service (IGS) [23] analysis centers for the processing of GPS and GLONASS, as well as the BeiDou-2 IGSO and MEO satellites [24]. This model introduced the right-handed DYB frame [3].
where a D , a Y , a B are accelerations at 1AU in the DYB frame. This model consists of three constant accelerations (D 0 , Y 0 , B 0 ) and two items B c and B s , and u is the argument angle of the latitude of the satellite.

ECOM2 SRP Model
Arnold et al. implemented the ECOM2 model [7], to which additional estimation parameters are added in the model of Beutler et al. [3]. Second-and fourth-order harmonic terms are added for the D-component. a D = D 0 + D c2 · cos 2∆u + D s2 · sin 2∆u + D c4 · cos 4∆u + D s4 · sin 4∆u a Y = Y 0 a B = B 0 + B c · cos ∆u + B s · sin ∆u. (2) Instead of the argument of latitude u, the angular argument ∆u = u − u s is used, where u s stands for the argument of latitude of the sun projected on the satellite's orbital plane.

Precise Orbit Determination Strategies
Other models involved in precision orbit determination, such as observation models, error models and parameter estimation models are described in Table 1. The distribution map of the stations of POD is shown in Figure 1. Among them, 56 MGEX/IGS stations [5] and 17 international GNSS Monitoring and Assessment System (iGMAS) stations provide BeiDou-2, BeiDou-3, and GPS observations [18]. The distribution map of the stations of POD is shown in Figure 1. Among them, 56 MGEX/IGS stations [5] and 17 international GNSS Monitoring and Assessment System (iGMAS) stations provide BeiDou-2, BeiDou-3, and GPS observations [18]. We selected one day (4 April 2019) to analyze the number of observations of BeiDou satellites. From Figure 2, it can be seen that: For the BeiDou-3 MEO, the average number of observations of C20, C21, and C22 was 10487, which had a factor that was of 1.5 times larger than that of other BeiDou-3 MEO satellites. For the BeiDou-2 MEO, the average number of observations was about 11758. In general, the number observations of BeiDou-3 MEO were smaller than that of the BeiDou-2 MEO satellite, especially for the C35 satellite.

C06
C07 C08  C09  C10  C11  C12  C13  C14  C16  C20  C21  C22  C23  C24  C25  C26  C27  C29  C30  C32  C33  C34  C35  --0   We selected one day (4 April 2019) to analyze the number of observations of BeiDou satellites. From Figure 2, it can be seen that: For the BeiDou-3 MEO, the average number of observations of C20, C21, and C22 was 10487, which had a factor that was of 1.5 times larger than that of other BeiDou-3 MEO satellites. For the BeiDou-2 MEO, the average number of observations was about 11758. In general, the number observations of BeiDou-3 MEO were smaller than that of the BeiDou-2 MEO satellite, especially for the C35 satellite. The distribution map of the stations of POD is shown in Figure 1. Among them, 56 MGEX/IGS stations [5] and 17 international GNSS Monitoring and Assessment System (iGMAS) stations provide BeiDou-2, BeiDou-3, and GPS observations [18]. We selected one day (4 April 2019) to analyze the number of observations of BeiDou satellites. From Figure 2, it can be seen that: For the BeiDou-3 MEO, the average number of observations of C20, C21, and C22 was 10487, which had a factor that was of 1.5 times larger than that of other BeiDou-3 MEO satellites. For the BeiDou-2 MEO, the average number of observations was about 11758. In general, the number observations of BeiDou-3 MEO were smaller than that of the BeiDou-2 MEO satellite, especially for the C35 satellite.

D-Acceleration Analysis
In yaw-steering mode, the solar panel of the satellite is always perpendicular to the sun. The solar radiation acceleration generated by the solar panel is constant regardless of the multiplication distance-factor of the satellite to the sun. Thus, the D-component variation can reflect the solar radiation acceleration variation generated by the satellite body with the β (beta) angle. Therefore, the difference of SRP between the BeiDou-2 and BeiDou-3 satellites can be reflected on the D0 (estimation of D-component of ECOM) variations, where C09 and C12 were selected as representatives for the BeiDou-2 satellite, and C20 and C30 were selected for the BeiDou-3 satellites from the CAST and SECM manufacturers, respectively. The D0 variations are shown in Figure 3. When the POD is carried out, the phase center offsets (PCO) of the BeiDou-3 MEO satellites are of crucial importance. However, the accurate PCOs of the BeiDou-3 satellites are still un-released. Therefore, the coarse parameters, recommended by iGMAS in 2018, were used and are shown in Table 2.

D-Acceleration Analysis
In yaw-steering mode, the solar panel of the satellite is always perpendicular to the sun. The solar radiation acceleration generated by the solar panel is constant regardless of the multiplication distance-factor of the satellite to the sun. Thus, the D-component variation can reflect the solar radiation acceleration variation generated by the satellite body with the β (beta) angle. Therefore, the difference of SRP between the BeiDou-2 and BeiDou-3 satellites can be reflected on the D0 (estimation of D-component of ECOM) variations, where C09 and C12 were selected as representatives for the BeiDou-2 satellite, and C20 and C30 were selected for the BeiDou-3 satellites from the CAST and SECM manufacturers, respectively. The D0 variations are shown in Figure 3.  From Figure 3, it can be seen that the BeiDou-2 satellites had a smaller D0 variation than the BeiDou-3 satellites, except for the larger scatters near the maximum absolute β angle and near the zero, and the range of the D0 variations was smaller than 0.5 nm/s 2 for C09 and C12. However, the systematic variation of D0 was about 5 nm/s 2 when |β| > 40° for C20.
When |β| < 20°, the variation of D0 was about 1.5 nm/s 2 for C30; for 20 < β < 40°, there was a linear systematic of D0, with a range of 1.5 nm/s 2 . Combined with the function of the D0 variations related to the β-angle, derived by Montenbruck et al. [10] from the Box-wing model, there are four parameters affecting the systematic variations of D0, which are C a αδ , C a ρ , S a αδ , and S a ρ (section 3.2).
Only the partial derivatives of S a αδ and S a ρ were a near linear-slope with respect to the β-angle. From Figure 3, it can be seen that the BeiDou-2 satellites had a smaller D0 variation than the BeiDou-3 satellites, except for the larger scatters near the maximum absolute β angle and near the zero, and the range of the D0 variations was smaller than 0.5 nm/s 2 for C09 and C12. However, the systematic variation of D0 was about 5 nm/s 2 when |β| > 40 • for C20.
When |β| < 20 • , the variation of D0 was about 1.5 nm/s 2 for C30; for 20 < β < 40 • , there was a linear systematic of D0, with a range of 1.5 nm/s 2 . Combined with the function of the D0 variations related to the β-angle, derived by Montenbruck et al. [10] from the Box-wing model, there are four parameters affecting the systematic variations of D0, which are a αδ C , a ρ C , a αδ S , and a ρ S (Section 3.2). Only the partial derivatives of a αδ S and a ρ S were a near linear-slope with respect to the β-angle. Therefore, it can be inferred that the C30 satellite has significant a αδ S and a ρ S parameters, resulting from the cuboid of the satellite bus, like the Galileo satellites.
A similar phenomenon was also found for the other satellites of CAST and SECM. Therefore, an a priori box-wing model is necessary to decrease the D0 variation of the BeiDou-3 satellites.

A Priori SRP for the BeiDou-3 Satellites
The box-wing model is an analytical model, with physical properties that reflect the actual solar radiation pressure. Therefore, the box-wing model can be an a priori model to compensate for the inadequacy of the empirical model, like ECOM.

SRP Acceleration Obtained by the Semi-Analytical Model
The box-wing model needs a priori mass, dimensions, and optical properties. For an individual surface, the acceleration vector of the inertial system can be expressed as [22]: where α, ρ, and δ are the physical optical properties of a satellite surface, named the fraction of the absorbed, specular, and reflected photons; A is the surface area; M is the total mass of the satellite; S = 1367 W/m 2 denotes the solar irradiance at the Sun-Earth distance; and c is the vacuum velocity of light.
→ e D is the vector from the satellite to the sun, and → e N is the normal vector of the illuminated surface. θ is the angle between → e D and → e N . The satellite body surfaces are covered by multi-layer insulation (MLI), the energy absorbed by the satellite surface is instantaneously re-radiated back into the space, according to Lambert's law, and the acceleration vector can be expressed as: For both Equations (3) and (4), a condition of α + δ + ρ = 1 is implicit. The performance of the box-wing model mainly depends on the accuracy of the surface areas and their optical properties, which is presently only available for GPS, Galileo and, partly, GLONASS. These problems limit the accuracy of the box-wing model. Rodriguez-Solano et al. used the box-wing model to improve the performance of GPS POD by adjusting the optical properties using the ground tracking data. The estimated parameters are absorption plus diffuse reflection (AD) as well as the specular reflection (R) for the illuminated surfaces (+X, +Z, −Z), the optical properties of the solar panel (SP), the constant acceleration (Y0), and the solar panel rotation lag (SB). Thus, this model was named the adjustable box-wing (ABW) [14]. Since the ABW model directly estimates the optical properties, the accuracy of the a priori optical properties is not very critical, unlike for the box-wing model. The metadata of the BeiDou satellites have so far not yet been published by an official organizer. Therefore, the ABW model is an alternative to replace the box-wing model for the BeiDou-3 satellite POD. To achieve the precise solar radiation accelerations of the box-wing, firstly, the precise accelerations of satellite were integrated from the illuminated panels and body-surfaces with the adjustable optical properties of ABW model. The precision of the accelerations was indirectly assessed by the accuracy of POD.
The strategies of data processing are similar to those presented in Table 1, and only the ABW model, with the nine parameters of +XAD, +XR, +ZAD, +ZR, −ZAD, −ZR, SP, Y0, and SB, was used. Considering the correlations between the nine parameters, the constraint of 0.1 was added to +XAD, +XR, +ZAD, +ZR, −ZAD, and −ZR, and a constraint of 1000 was added to SP, Y0, and SB. For the SECM's satellites, the areas of the satellites were published in the ION GNSS+ conference in 2018 [19] ( Table 3), and only the optical properties were used for BeiDou-2 values as an approximation, which are shown in Table 3. For the CAST's satellites, the initial optical properties and areas were used for the BeiDou-2 values instead (Table 4, [20]).

Estimation of the A Priori Cuboid SRP Model
The box-wing model can be reformulated as a function of the sun-satellite-earth angle (ε), where the attitude obeys the nominal yaw-steering law [25]. The contributions of the satellite bus to the acceleration can be isolated as the cube, cuboidness, and z-asymmetry, which was described as the parameters of a αδ , and superscript C, S, and A represent the cubic, stretched body shape, and asymmetric [6,10]. The detailed expression can be written as: where The acceleration of box-wing model is converted from the inertial system (Equation (4)) to the DYB frame (ECOM1 model). The box-wing acceleration can be expressed as: Remote Sens. 2019, 11, 1605 where a sp is the main part of the D-acceleration, contributed and engendered by the solar panels, which was nearly a constant for an arc of the POD in the yaw-steering period. So far, the function expressions for the box-wing accelerations, relative to the ε(eps) angle, are all completed.
There  To cover a whole range of the β-angle, long range data processing is necessary.
Considering the efficiency of the data processing, we used a step of 10 days to cover the β-angle range. The period is selected from the DOY of 218 in 2018 to the DOY of 058 in 2019. The estimation of the cuboid parameters of C20 and C30 are shown in Figure 5. 16 12 Based on the function model of Equations (6) and (7), the parameters of a αδ C , a αδ S , a αδ A , a ρ C , a ρ S , a ρ A ,a sp were estimated by a separated least-squares adjustment, where the D-and B-accelerations derived from the ABW model were used as the observations. Therefore, the total parameters are a αδ C , a αδ S , a αδ A , a D-and B-accelerations derived from the ABW model were used as the observations. Therefore, the total parameters are , , a a a a a a αδ αδ αδ ρ ρ ρ , , , (named cuboid parameters) and sp a (named constant parameter). To cover a whole range of the β-angle, long range data processing is necessary. Considering the efficiency of the data processing, we used a step of 10 days to cover the β-angle range. The period is selected from the DOY of 218 in 2018 to the DOY of 058 in 2019. The estimation of the cuboid parameters of C20 and C30 are shown in Figure 5. For the BeiDou-3 satellites, a clear separation between the scatters of a αδ C and a ρ C , which are the main parameters for the cuboid model, can be seen from Figure 5. For the a αδ S and a ρ S , it can also be seen that the coefficients were significant for the C30 satellites. Nevertheless, the a αδ S and a ρ S were close to zero for C20, and the scatters were still stable. For the a ρ A , the estimations of both C20 and C30 were near to zero. However, for the a αδ A , the scatters were larger than the other cuboid parameters, for both C20 and C30. Considering the poorer accuracy and small magnitude, we ignored a αδ A in the a priori cuboid box-wing model (named the cuboid model). The standard deviation (STD) of the corresponding time series was used to evaluate the precision of the cuboid parameter. For C20, the STD values of the corresponding time series of a αδ C ,a αδ S , a αδ A , a ρ C , a ρ S , and a ρ A were 0.08, 0.14, 0.60, 0.11, 0.17, and 0.32 nm/s 2 . For C30, the STD values of the corresponding time series of a αδ C , a αδ S , a αδ A , a ρ C , a ρ S , and a ρ A were 0.08, 0.14, 0.60, 0.11, 0.17, and 0.32 nm/s 2 , respectively. The cuboid estimations of the other BeiDou-3 satellites showed the similar results with C20 and C30, and the final cuboid models of CAST and SECM were obtained by averaging each type of satellite. The final estimations, as well as the STD of the cuboid models, are shown in Table 5.

Validations
The validations of the cuboid model were based on the D0 variation, orbit accuracy, and clock offset accuracy of the BeiDou-3 satellites. The D0 variation was used to verify the systematic variation of C30. The orbit accuracy can reflect the performance of the SRP model, and the precision of 24 h orbit overlapping and SLR residuals were selected as the two indicators. Meanwhile, considering that the ECOM2 model can also improve the orbit accuracy of Galileo and QZSS, we also compared the Cuboid + ECOM1 with the ECOM2 model. For the clock offset accuracy, the fitting curve of the clock offset, RMS of the fitting residuals, and STD of the clock offset overlapping were selected as the indicators.

D0 Variation
The smaller the variation of the estimated D0 parameter with the β-angle, the better the performance of the SRP model, when it is close to the actual of the solar radiation pressure. Therefore, the variation of the D0 parameter with the β-angle can be one of the indicators for evaluating the SRP model. The results are shown in Figure 6.

D0 Variation
The smaller the variation of the estimated D0 parameter with the β-angle, the better the performance of the SRP model, when it is close to the actual of the solar radiation pressure. Therefore, the variation of the D0 parameter with the β-angle can be one of the indicators for evaluating the SRP model. The results are shown in Figure 6. For the ECOM2 model, the scatters of D0 were larger than those for ECOM1 and Cuboid + ECOM1. The main reason was that the function of the D-accelerations was not a constant of ECOM2 (Equation (2)) for an arc-length. Therefore, only the D0 variations of the ECOM1 and Cuboid + ECOM1 were compared. For C20, the STD of D0 was 1.28 and 0.85 nm/s 2 for the ECOM1 and Cuboid + ECOM1 solutions, respectively, and we found that the Cuboid + ECOM1 model can reduce the variation by about 34%. However, for 0 < |β| < 40°, a linear systematic variation, with a slope of 0.028 nm/s 2 /°, was found. This might be improved by accurate body-dimensions, rather than the "guess" values of the CAST's satellites. For 14 < β < 40°, a linear systematic variation, with a slope of 0.048 nm/s 2 /°, was found for C30, when the pure ECOM1 model was adopted for POD. After adding the cuboid model, the linear slope was reduced to 0.005 nm/s 2 /°. This also proved the previous inference on the D0 variation of C30 (Figure 3), i.e., that the systematic variation of D0 was greatly reduced by the S a αδ and S a ρ parameters of the a priori cuboid model. For the ECOM2 model, the scatters of D0 were larger than those for ECOM1 and Cuboid + ECOM1. The main reason was that the function of the D-accelerations was not a constant of ECOM2 (Equation (2)) for an arc-length. Therefore, only the D0 variations of the ECOM1 and Cuboid + ECOM1 were compared. For C20, the STD of D0 was 1.28 and 0.85 nm/s 2 for the ECOM1 and Cuboid + ECOM1 solutions, respectively, and we found that the Cuboid + ECOM1 model can reduce the variation by about 34%. However, for 0 < |β| < 40 • , a linear systematic variation, with a slope of 0.028 nm/s 2 / • , was found. This might be improved by accurate body-dimensions, rather than the "guess" values of the CAST's satellites. For 14 < β < 40 • , a linear systematic variation, with a slope of 0.048 nm/s 2 / • , was found for C30, when the pure ECOM1 model was adopted for POD. After adding the cuboid model, the linear slope was reduced to 0.005 nm/s 2 / • . This also proved the previous inference on the D0 variation of C30 (Figure 3), i.e., that the systematic variation of D0 was greatly reduced by the a αδ S and a ρ S parameters of the a priori cuboid model.

Satellite Orbits
The RMS of the 24 h orbit overlapping, which was the orbit of the middle of the day in the first 3d arc solution compared with those of the first day in the second 3d arc solution, and can also reflect the internal coincidence accuracy of the orbit [29]. A comparison of the RMS of the 24 h orbit overlapping using the three models is shown in Table 6, where pseudo-random noise (PRN) and average (AVG) were used. The RMS for each of satellite was given for the along-track (A), cross-track (C), and radial (R) components. From Table 6, for the ECOM2 model, the averages of RMS were 8.56, 5.31, and 2.71 cm for the A-, C-, and R-components, respectively. Compared with the ECOM1 model, the degradation was about 1 cm for the R-component. However, the improvement of the C-component was about 1 cm, and the 3-dimension (3D) RMS had a slight improvement of about 0.5 cm. For the Cuboid + ECOM1 model, the average improvements of the ACR-components were 0.99 (11.6%), 1.30 (20.8%), and 0.17 cm (9.5%), and the C-component had a higher improvement than the A-and R-components. In general, an average precision of 7.54, 4.96, and 1.62 cm of the 24 h orbit overlapping can be achieved for BeiDou-3 POD, with a 3D precision of 9.2 cm.
SLR is an independent technique operating in the optical frequency range, which is virtually free from the systematic effects related to ionosphere and troposphere delays, ambiguities, and clocks [7]. Therefore, SLR observations are contaminated by only a few error sources and are widely used to evaluate the accuracy of the GNSS orbit [30]. While the high accuracy of the orbits can be determined by the GNSS technique, some orbit modeling deficiencies remain in the SLR residuals, with the characteristic pattern first noted by Urschl et al. [31]. SLR residuals are the differences between the observed (using SLR) and computed distances (derived from the GNSS technique orbits) at the measurement epochs [8]. The station coordinates of the SLR stations are fixed to the a priori reference frame SLRF2008 [32]. The SLR observations are corrected for the relativistic effects, troposphere delays, and for the offset of the laser retroreflector arrays (LRA), related to the satellites' centers of mass [7]. For the BeiDou-3 MEO satellites, the LRA offsets are shown in Table 7. Residuals larger than ±0.4 m were removed, as gross error, when data editing. The one-way SLR residuals are shown in Figure 7. The blue line was the cubic polynomial fitting curve of the SLR residuals. From the fitting curve, it can be seen that: For the ECOM1 model result, a systematic variation between the lower and higher ε-angle was about 14 cm for the C20 and C30 satellites. When the ECOM2 and Cuboid + ECOM1 models were used, the variation was significantly reduced to about 5 cm. Thus, the ECOM2 and Cuboid + ECOM1 models can reduce the systematic variation of the SLR residuals. Compared with the RMS of the SLR residuals of the ECOM1 model, shown in Table 8, for the CAST's satellites, the STD was improved by about 0.78 cm (15.4%), and the RMS was slightly improved by about 0.54 cm (10.3%) for the ECOM2 model. For the Cuboid + ECOM1 model, the STD was improved by about 1.0 cm (19.3%), and the RMS was slightly improved by about 0.43 (8.7%). However, the offsets of C20 and C21 were increased for about 0.55 and 1.17 cm for the ECOM2 and Cuboid + ECOM1 models, respectively. The possible reasons were the unmodeled acceleration or unmodeled errors, which might also have been caused by the deviation of the satellite-dimensions of CAST, further leading a bias of the cuboid model. Thus, further work will need to be conducted in this direction. For the SECM's satellites, compared with the ECOM1 model, the RMS of the ECOM2 and Cuboid + ECOM1 models was improved by about 0.7 cm (10.9%) and 1.6 cm (25.6%), respectively. The Cuboid + ECOM1 model had a higher accuracy than the ECOM2 model. Compared with the RMS of the SLR residuals of the ECOM1 model, shown in Table 8, for the CAST's satellites, the STD was improved by about 0.78 cm (15.4%), and the RMS was slightly improved by about 0.54 cm (10.3%) for the ECOM2 model. For the Cuboid + ECOM1 model, the STD was improved by about 1.0 cm (19.3%), and the RMS was slightly improved by about 0.43 (8.7%). However, the offsets of C20 and C21 were increased for about 0.55 and 1.17 cm for the ECOM2 and Cuboid + ECOM1 models, respectively. The possible reasons were the unmodeled acceleration or unmodeled errors, which might also have been caused by the deviation of the satellite-dimensions of CAST, further leading a bias of the cuboid model. Thus, further work will need to be conducted in this direction. For the SECM's satellites, compared with the ECOM1 model, the RMS of the ECOM2 and Cuboid + ECOM1 models was improved by about 0.7 cm (10.9%) and 1.6 cm (25.6%), respectively. The Cuboid + ECOM1 model had a higher accuracy than the ECOM2 model.

Satellite Clock Offsets
The accuracy of the satellite clock offset is consistent with that of the satellite orbit, and the systematic error of the SRP model can also be reflected in the satellite clock offset. Therefore, the quadratic polynomial fitting was performed on each of the 3d-arc length of the satellite clock offset, and then the residuals were used to analyze the accuracy of the ECOM1, ECOM2, and Cuboid + ECOM1 models. The fitted residuals are shown in Figure 8, and the fitted RMS is shown in Table 9.    Comparing Figures 7 and 8, it can be seen that the clock residuals were consistent with those of the satellite orbit and also had a systematic variation with respect to the ε angle. The ECOM2 and Cuboid + ECOM1 model can reduce the systematic variation of the residuals of C20 and C30. For C20, the fitting curve of the ECOM2 had a smaller intercept and slope than Cuboid + ECOM1. However, the fitting RMS of Cuboid + ECOM1 was slightly smaller than that of ECOM2. For C30, the fitting curve of the Cuboid + ECOM1 had a smaller intercept, slope, and RMS than ECOM2. For the Cuboid + ECOM1 model, the improvements of RMS were about 0.06 (24%) and 0.04 ns (14%) for C30 and C20, respectively. In general, the average RMS of ECOM2 was larger than that of the Cuboid + ECOM1 model. The main reason was that the greater number of parameters of ECOM2 reduces their observability and result in less stable results than those of the ECOM1 model, especially for the small observations of the BeiDou-3 satellites.
Considering that the cuboid model can improve the precision of orbit overlapping, the clock offset overlapping was also analyzed. Considering that there is a smaller number of observations (Figure 2), the ratio of the effective clock offset epochs to the theoretical epochs (864 for three days, with an interval of 300 s) is low, for instance, C24 had a ratio of 71%. Therefore, to increase the number of valid epochs to calculate the STD and mean values of the clock offset overlapping, the overlapping arc was set to 48 h, rather than 24 h. The precision of the 48 h clock overlapping is shown in Table 10. The STD and offset were obtained by a quadratic difference between the satellite clock offset series and reference clock offset series. Compared with the ECOM1 model, the average improvements of the STD and offset were about 0.01 ns (4%) and 0.06 ns (11%), for the ECOM2 model, and were 0.02 ns (13%) and 0.06 ns (16%) for the Cuboid + ECOM1 model, respectively.

Conclusions
Considering that the metadata of the BeiDou-3 satellite is not fully released at present, we adopted the adjustable box-wing model (ABW) model, for which the optical properties are not very critical, to obtain the precise solar radiation accelerations. The cuboid parameters were directly estimated by the precise solar radiation accelerations. Compared with the method that uses estimated optical properties, together with mass and body-dimensions, to calculate the cuboid parameters, the directly estimated method can avoid the effect of the propagation of satellite mass and dimension errors.
However, this method for estimating the cuboid parameters is not applicable to the orbit-normal mode (ON), because Equations (6) and (7) of box-wing in the DYB-coordinate system are based on the yaw-steering (YS) mode. Therefore, when estimating the cuboid parameters, the results of |β| < 3 • were removed. However, when validating the cuboid model, it was found that the orbital accuracy during the period of |β| < 3 • is about 4.8 cm of the RMS of the satellite laser ranging (SLR) residuals, which is consistent with the accuracy of the YS mode. Thus, the a priori cuboid model is also applicable for the continuous yaw-steering (CYS) mode. For details concerning the ON mode, please refer to [10][11][12].
The cuboid model can improve the accuracy of precise orbit determination (POD), and the systematic variation of the SLR residuals between the lower and higher sun-satellite-earth angle (ε), was reduced from 14 cm to less than 5 cm. The cuboid model can also improve the clock offset accuracy, and the improvements of the RMS of the quadratic polynomial fitting residuals was about 0.06 (24%) and 0.04 ns (14%) for C30 and C20, respectively.
The phase center offsets (PCO) had a high correlation with the ECOM1 model, especially for the D0 parameter, although the BeiDou-2 IGSO and MEO had a smaller D0 variation than BeiDou-3 MEO, and the horizontal PCO scatters still had a systematic variation with respect to the β-angle [33]. A tight D0 constraint can reduce the high correlation between the horizontal PCO and D0 parameters and obtain more stable horizontal PCO estimations, where D0 can be provided by an a priori box-wing model when estimating Galileo PCOs [34]. Thus, the cuboid model might have a better performance in the estimations of PCO for the BeiDou-3 satellites. Currently, the accuracy of POD of the BeiDou-3 satellite is still limited by the number of tracking stations. Therefore, the accuracy of BeiDou-3 POD will be improved by more stations, and more observations, as well as accurate PCOs in the future.