Estimation of the Antenna Phase Center Correction Model for the BeiDou-3 MEO Satellites

Satellite antenna phase center offsets (PCOs) and phase variations (PVs) for BeiDou-3 satellites are estimated based on the tracking data of the Multi-GNSS Experiment (MGEX) and the international GNSS Monitoring and Assessment System (iGMAS) network. However, when estimating the (PCOs) of BeiDou-3 medium Earth orbit (MEO) satellites by pure Extending the CODE Orbit Model (ECOM1), the x-offset estimations of the PCOs have a systematic variation of about 0.4 m with the elevation of the Sun above the orbital plane (β-angle). Thus, a priori box-wing solar radiation pressure (SRP) model of BeiDou-3 MEO was assisted with ECOM1. Then, the satellite type-specific PCOs and common PVs were obtained. The estimations of PCOs and PVs were compared with the MGEX PCOs from the precise orbit and clock offset. When the MGEX PCOs were used, the root mean square (RMS) of 24 h overlap was 6.76, 4.36, 1.46 cm, in along-track, cross-track, and radial directions, respectively; the RMS and standard deviations (STD) of the 24 h clock offset overlap were 0.28 and 0.15 ns; the fitting RMS of the 72 h clock offset of the quadratic polynomial was 0.243 ns. After comparing this with the estimated PCOs and PVs, the RMS of the 24 h orbit overlap was decreased by 6.5 mm (10.54%), 1.8 mm (4.4%), and 1.1 mm (8.03%) in the along-track, cross-track, and radial directions, respectively; the RMS and STD of the 24 h clock offset overlap were decreased by 0.024 ns (8.6%) and 0.020 ns (13.1%), respectively; the fitting RMS of the 72 h clock offset of the quadratic polynomial was reduced by about 0.016 ns (6.5%).


Introduction
As of September 2019, 23 satellites of the BeiDou-3 satellites were launched into orbit, which are 20 medium Earth orbit (MEO) satellites, one Geostationary Earth Orbit (GEO) satellite, and two Inclined Geosynchronous Orbit (IGSO) satellites. By the end of 2020, the global constellation of BeiDou-3 will be complete, and will provide a global navigation satellite system (GNSS) positioning, navigation, and timing services, together with GPS, GLONASS, and Galileo [1,2]. Precision orbit and clock offsets in the same reference frame are a prerequisite for multi-GNSS applications. In multi-GNSS data processing, in addition to the same International Earth Rotation and Reference Systems Service (IERS) protocol, the unification of the framework is mainly reflected in the consistency of the antenna phase center correction (PCC) and the coordinates of the core station of the framework in Solution INdependent EXchange Format (SINEX). The PCC consists of the antenna phase center offsets (PCOs) and phase variations (PVs).
The International GNSS Service (IGS) [3] has released the phase center offsets (PCOs) and phase variations (PVs) for GPS, GLONASS, and Galileo satellites in igs14.atx, which is aligned to the IGS2014.

PCO Parameters Model
The satellite antenna PCO is a vector from the mass center to phase center, which consists of x-offset, y-offset, and z-offset in a satellite body-fixed coordinate system [22]. The correction of PCO parameters on the observed distance is shown in Equation (1): ∆ρ(α, η) = dx· sin α sin η + dy· cos α sin η + dz· cos η (1) where α is the azimuth and η is the nadir angle seen from the satellite to the station, and the azimuth α is started from the y-axis toward the x-axis of the satellite body-fixed system when viewing from the station. dx, dy, and dz are x-offset, y-offset, and z-offset, respectively.

PV Parameters Model
Due to the high correlation between the satellite antenna PVs and z-offset, raw PVs were estimated that correspond to the following sum [19,21]: where ∆z was the correction of the a priori value of the z-offset, to prevent the normal equation from being singular and considering the correlation with the satellite clock offsets, an a priori constraint is added: The PVraw(η i ) are estimated satellite-by-satellite as piece-wise constants. After PVraw(η i ) parameters are obtained, a separate adjustment model, shown as Equation (2), is established to derive PV(η i ) and ∆z parameters. The criterion of least square adjustment is as follows: where a is a constant parameter, and n is the maximum integer nadir-angle of PV(η i ). The residuals of the least square adjustment are PVs, and the PVs and z-offset datum is n i=0 PV(η i ) = 0. This model mainly used for GPS and GLONASS to estimate satellites phase center correction [4,19].

Solar Radiation Pressure Model
The empirical Center for Orbit Determination in Europe (CODE) orbit model is expressed in a Sun-oriented reference frame, where D points to the Sun, Y goes along with the solar panel axis, and B completes a right-handed system [10]. The ECOM1 model can be expressed as: where D 0 , Y 0 , B 0 are three constant parameters, B c and B s are cosine and sine terms in B direction, and u is the argument of latitude. ECOM1 was widely adopted for the SRP of GPS and GLONASS for most of the Multi-GNSS Experiment (MGEX) analysis centers [7,23].
Although the ECOM2 model [24] was better used than the ECOM1 for precise orbit determination of the cuboid satellites [15,25], it was not applied for estimating PCOs, because the increased parameters of second-and fourth-order harmonic terms for the D-component caused a higher correlation between PCOs and SRP parameters [5].

A Priori Box-Wing Model for BeiDou-3 MEO
During the nominal yaw-steering period [22], the box-wing model is reformulated as Equation (6) [11]: where a D and a B are the accelerations in D and B directions, respectively; ε is the sun-satellite-earth angle; a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , a 7 are the model coefficients.
The a priori box-wing model can provide a priori value to the pure ECOM1 model. It can impose reasonable constraints on the ECOM1 parameters when estimating PCOs and obtain more stable x-offset parameters. There are two satellite manufacturers for the BeiDou-3 satellites: the China Academy of Space Technology (CAST) and the Shanghai Engineering Center for Microsatellites (SECM) [26,27]. The optical properties and shape sizes of the two types of satellites are different, so the coefficients of the a priori box-wing model are also different. The a priori box-wing model of the two types of BeiDou-3 MEO satellites was estimated by Yan et al. [15], and the coefficients are listed in Table 1.

Strategies
The distribution map of the stations of precise orbit determination for the BeiDou-3 satellites is shown in Figure 1. Among them, 56 MGEX/IGS stations and 17 international GNSS Monitoring and Assessment System (iGMAS) stations provide BeiDou-3 and GPS observations. The general data processing strategies for estimating PCOs, PVs, and precise orbit determination of the BeiDou-3 satellites are listed in Table 2, mainly including observation models, error models, and parameter estimation models.
For special data processing, the differences will be described before the analysis of results. The conventional MGEX values of the block-specific PCOs of the BeiDou-3 satellites are listed in Table 3.
There are four frequencies, named C01, C02, C06, and, C07, in the conventional MGEX PCOs model. Duplicate entries for C01 and C02 support the use of RINEX 3.01 and 3.02 observation codes of B1I, and the code of the subsequent RINEX versions, respectively. The conventional MGEX PCOs are both the initial value when they were estimated, and the reference model when the estimated PCOs were validated.
The distribution map of the stations of POD is shown in Figure 1. Among them, 56 MGEX/IGS 142 stations and 17 international GNSS Monitoring and Assessment System (iGMAS) stations provide

PCOs and PVs Estimations
Considering that the factors affecting the three components are different, we divided the PCOs into horizontal and vertical PCOs. The horizontal PCOs are composed of the x-offset and y-offset, which are affected by the SRP model. Here, the impact of different SRP models on the x-offset are analyzed.

Horizontal PCO Estimations
When PCOs were estimated by the pure ECOM1 (Figure 2a), the x-offset time series of the C21 satellite produced a systematic variation of about 40 cm with the β-angle, which is similar to the results of the Galileo satellites [5]. Therefore, the a priori box-wing model ( Table 1) was used to assist the ECOM1 model ( Figure 2b). Although more stable x-offset estimations were obtained, there were still large scatters when the absolute value of the β-angle was large. Further, when a constraint of 1 nm/s 2 was added to the D 0 parameter (Figure 2c), the x-offset series appeared to be a stationary sequence that was almost irrelevant to the β-angle.   Comparing the average x-offset of the three SRP models (Figure 2a-c), it was found that there was a significant difference of about 11 cm between the pure ECOM1 and the a priori model-assisted ECOM1 solutions. For the results using the a priori model, the difference between the constraint and non-constraint was 1.3 cm. For the other CAST satellites, similar results were found. However, no significant systematic variations in the x-offset were found for SECM satellites, such as C25, C26, C27, C28, C29, and C30, which had smaller ranges of β-angle of ±60 • and ±40 • . For the three SRP models, there was no significant effect on the y-offset time series.  Figure 2b, corresponding to the larger scatters of x-offset, |β| > 60 • , the results were the gray-shaded area in Figure 3, and the RMS of the orbit overlap in along-track (A), cross-track (C), and radial (R) directions were computed and shown in Table 4. Compared with the results of estimating PCOs and non-constraint of D 0 parameter in scheme (b), the RMS of scheme (c) was reduced by 0.65 (6.8%), 6.36 (59.2%), and 0.4 (16.3%) cm, in A-, C-, and R-components, respectively. Therefore, it can be inferred that the correlation between the x-offset and the C-component of the orbit was higher than the A-and R-components for |β| > 60 • . After adding the constraint to the D 0 parameter, the RMS of orbit overlap was close to the reference result (scheme a), and the differences between scheme (c) and scheme (a) were 0.79, −1.72, and 0.25 cm, for A-, C-, and R-components, respectively.    After obtaining the daily horizontal PCOs time series, the average of each satellite was calculated as the estimation of this satellite. The standard deviation (STD) was also used as an indicator to measure the repeatability of the daily PCOs series. The specific results are shown in Table 5. From the STD values, a centimeter precision level of horizontal PCOs was obtained. Considering that the PCOs of the same type of satellite are close, such as C19-C24 for CAST satellites, and C25-C30 for SECM, to improve the reliability of the PCOs, the average of the same type of satellite was taken as the final type-specific satellite PCOs. The STD of inner-type PCOs was also calculated. Compared with other satellites, C32, C33, and C35 were launched into orbit relatively late and had a short period of valid observations, resulting in poor accuracy of PCOs, the three satellites were excluded when calculating the final type-specific PCOs. The final PCOs and STDs are shown in Table 6.

Vertical PCO and PV Estimations
The main factors affecting the z-offset are the PVs and satellite clock offset. Thus, the PVs and z-offset were estimated in two steps. Firstly, the z-offset was estimated while the PVs were fixed as zero for all nadir-angles, then, the z-offset was fixed as the estimation, and the PVraws were estimated as a linear piece-wise constant model. Subsequently, PVs were determined by removing the cosine variation and the constant according to Equation (4), which yielded final PVs and correction to the z-offset. The z-offset estimations and the STD values of BeiDou-3 MEO satellites in the first step are listed in Table 7 when the PVs were fixed to the value of zero. The z-offset estimations of C21 and C30 were selected as the representatives of CAST and SECM satellites and are shown in Figure 4. For the CAST satellites, the largest difference in the z-offset was 156 mm and the average value of STD was 124 mm, which indicated a good consistency between z-offsets of satellites. For the SECM satellites, except for the C25, C29, and C35 satellites, the maximum difference of z-offset was 139 mm and the average value of the STD was 121mm, which indicated a good consistency of z-offset between the same type of satellite. Therefore, the average values of 2368.38 and 1596.38 mm were fixed as the z-offset estimations to solve PVs for CAST and SECM satellites, respectively.    For each satellite, PVraws were obtained from the average of the daily PVraws. The day 232 repeatability was 1-2 mm at the 0° nadir-angle, 1 mm at 1° and 13°, and less than 1 mm for the other 233 nadir-angles. The PVraws of each satellite are shown in Figure 5. For the average PVraws of all 234 satellites, the STD was 1.6 mm at the 0° nadir-angle, 1 mm at 1° and 13°, and less than 1 mm for the 235 others. The reasons were mainly that the observations at 0 and 1° nadir angles were few, and the 236 weight of the observations at the 13° was reduced. For each satellite, PVraws were obtained from the average of the daily PVraws. The day repeatability was 1-2 mm at the 0 • nadir-angle, 1 mm at 1 • and 13 • , and less than 1 mm for the other nadir-angles. The PVraws of each satellite are shown in Figure 5. For the average PVraws of all satellites, the STD was 1.6 mm at the 0 • nadir-angle, 1 mm at 1 • and 13 • , and less than 1 mm for the others. The reasons were mainly that the observations at 0 and 1 • nadir angles were few, and the weight of the observations at the 13 • was reduced.
For each satellite, PVraws were obtained from the average of the daily PVraws. The day repeatability was 1-2 mm at the 0° nadir-angle, 1 mm at 1° and 13°, and less than 1 mm for the other nadir-angles. The PVraws of each satellite are shown in Figure 5. For the average PVraws of all 234 satellites, the STD was 1.6 mm at the 0° nadir-angle, 1 mm at 1° and 13°, and less than 1 mm for the 235 others. The reasons were mainly that the observations at 0 and 1° nadir angles were few, and the 236 weight of the observations at the 13° was reduced.   )   C19  C20  C21  C22  C23  C24  C25  C26  C27  C28  C29  C30  C32  C33  C34 C35 average Based on Equation (4), a separate least-squares adjustment was carried out to derive the unmodeled z-offset and PVs. After the adjustment, the PVs are shown in Figure 6a. Then, the average corrections of z-offsets are shown in Figure 6b. The final PVs and the z-offsets of the BeiDou-3 MEO satellites are listed in Tables 8 and 9. Considering that the z-offsets were different from other satellites, the estimations of C25 and C29 were given separately, and the type-specific z-offsets and STDs were given. For the CAST and SECM, the z-offset estimations had a good internal consistency with STDs of 35 and 51 mm, respectively.  Table 8 and Table 9. Considering that the z-offsets were different from other 243 satellites, the estimations of C25 and C29 were given separately, and the type-specific z-offsets and

Validations
To verify the impact of different satellite antenna phase center correction models on the precise orbit determination and clock offset estimation, three inspectors of overlap orbit, overlap clock offset, and clock offset fitting were selected. Three schemes were designed for comparison: Scheme 1: Precise orbit determination by the conventional MGEX PCOs model (Table 3), named S1; Scheme 2: Precise orbit determination by the estimated PCOs-only model (Table 9), named S2; Scheme 3: Precise orbit determination by the estimated PCOs and PVs models (Tables 8 and 9), named S3; The selected time was from the day of the year (doy) of 001 to 180 in 2019. Precise orbit determination was implemented according to schemes 1, 2, and 3, respectively. Compared to Table 2, the difference in the strategy was the PCC of the BeiDou-3 MEO satellites.

Orbit Precision
Since the 3-day arcs were adopted for the precision orbit determination, the difference in the overlap between the adjacent two solutions can be used to measure the internal coincidence of the orbit. Specifically, the orbit of the middle day of the 3-day solution is taken as a reference to measure the precision of the third day of the previous 3-day solution [33,34], and the overlap RMS values are given in along-track (A), cross-track (C), and radial (R) directions, are shown in Table 10. It can be seen from Table 10 that the precision of the orbit overlap of the C35 satellite was worse than that of other satellites, in which the RMS of A-direction reached 10 cm, and the radial direction also reached 2.7 cm. The observations of the C35 were less than the others, resulting in poor accuracy of precise orbit determination. The averages of RMS in the three schemes were computed. Compared with the MGEX PCOs (S1), the precision of the overlap with the estimated PCOs-only (S2) was increased by 6.1 mm (9.08%), 1.8 mm (4.10%), and 1.3 mm (8.64%) in the A-, C-, and R-directions, respectively; the precision of the overlap using the estimated PCOs and PVs (S3) increased by 6.5 mm (10.54%), 1.8 mm (4.4%), and 1.1 mm (8.03%) in the A-, C-, and R-directions, respectively. The improvement in the A-direction was more significant than the C-and R-directions when the estimated PCOs were used.
To quantitatively analyze the number of effective observations, and the contribution of the MGEX and iGMAS networks for the BeiDou-3 satellites, 36 and 14 stations were selected from MGEX and iGMAS, respectively. The time was randomly selected as day 150 of the year. The number of effective observations after POD is analyzed in Figure 7. It can be seen that C19, C20, C21, C22, and C28 had the most observations, with about 8300 from MGEX, and about 2900 from iGMAS. C23-C27 and C29-C34 had about 4500 observations from MGEX and 2700 from iGMAS. C35 had a smaller number of observations, with about 690 from MGEX. An expected ratio of the number of observations of MGEX and iGMAS was about 2.57 (36/14). A measured ratio of the number of observations of MGEX and iGMAS was computed by the actual number after POD. C23-C27 and C29-C34 had a lower magnitude (measured ratio < expected ratio) for about 1.57 of the number of observations. In particular, C35 had the lowest measured ratio of 0.23. The loss of the pseudorange and phase measurements of the B3I frequency mainly caused a lower measured ratio than expected. C19, C20, C21, C22, and C28 had a normal level of number of observations from the MGEX, with a measured ratio of 2.75. The main reason for this is that the iGMAS receiver has undergone software and hardware updates and has a better ability to receive BeiDou-3 satellite observations. However, the receivers of most of the MGEX stations have not been updated, and some of the receivers that have been replaced or updated have maintained good reception capabilities during the data processing period in this study. As of November 2019, after most of the MGEX receivers have been upgraded, data reception has returned to normal.

303
The precision of the 24 h clock offset overlap was taken to evaluate the accuracy of different PCC 304 models (S1, S2, and S3). The STD and RMS were obtained by a quadratic difference between the 305 satellite clock offset series and reference clock offset series [35], where the reference clock of BREW 306 was chosen. The specific clock offset precision is shown in Table 11.

Clock Offset Precision
The precision of the 24 h clock offset overlap was taken to evaluate the accuracy of different PCC models (S1, S2, and S3). The STD and RMS were obtained by a quadratic difference between the satellite clock offset series and reference clock offset series [35], where the reference clock of BREW was chosen. The specific clock offset precision is shown in Table 11. The RMS and STD of the C35 satellite were 0.34 ns and 0.25 ns, which showed a poorer accuracy than other satellites. One possible reason was the lack of observations (seen from Figure 7), another is a poorer accuracy of PCOs of the C35 compared to other satellites of SECM. After getting sufficient observations, further analysis of the specific reasons for this is required. Compared to the MGEX PCOs result (S1), the average RMS and STD of the estimated PCOs-only (S2) were decreased by 0.016 ns (5.8%) and 0.017 ns (11.3%), respectively. For the estimated PCOs and PVs result (S3), the RMS and STD of the clock offset were decreased by 0.024 ns (8.6%) and 0.020 ns (13.1%), respectively. The RMS of quadratic polynomial fitting of 72 h clock offset overlap is shown in Table 12. observations were more than other satellites of the same type, and the minimum fitting RMS of 0.153 ns was achieved; the remaining SECM satellites had a fitted RMS of 0.171-0.269 ns. Overall, the clock offset fitting RMS of the SECM satellite was smaller than the result of the CAST satellite. Compared with the results of MGEX PCOs (S1), the estimated PCOs-only (S2) were used to reduce the RMS of the clock offset by 0.01 ns (3.9%), and the estimated PCOs and PVs (S3) were used to reduce fitting RMS by about 0.016 ns (6.5%).

Conclusions
When estimating the phase center offsets (PCOs) in-orbit, the x-offset parameter may have systematic variation within the β-angle. In this case, the average of the longer-term data was generally used as the final estimation of the x-offset. Considering that the period of the β-angle is half a year for the MEO satellites, at least one year of data is usually required. The systematic variation of x-offset with the β-angle is mainly due to the correlation between the solar radiation pressure (SRP) model. The absolute value of the β-angle is larger, the correlation is stronger. To solve this problem, a priori SRP model is added to assist the ECOM1, and apply a suitable D 0 constraint to reduce the correlation. Therefore, a stationary sequence of the x-offset estimations that were weakly related to the β-angle was obtained, and the final estimation of the x-offset can be obtained for a relatively short period, rather than a half year.
In this work, the average of the daily x-offset estimation of the China Academy of Space Technology (CAST) satellite was −0.321 ± 0.092 m, with a systematic variation of about 0.4 m by the ECOM1. After adding the a priori SRP model [15], a constraint of 1 nm/s 2 was also applied to the D 0 parameter, and a stationary sequence with the average of −0.227 ± 0.044 m was obtained. Comparing the two averages and STDs, the average of x-offset estimations had a difference of about 0.1 m, and the stability was improved by 50%. For the y-offset parameter, whether or not the a priori SRP model was added had no significant effect on the final estimation.
The estimated PCOs and PVs were compared with the Multi-GNSS Experiment (MGEX) PCOs from the precise orbit and clock offset. The Root Mean Square (RMS) of the 24 h orbit overlap with estimated PCOs and PVs was decreased by 6.5 mm (10.54%), 1.8 mm (4.4%), and 1.1 mm (8.03%) in the along-track, cross-track, and radial directions, respectively; the RMS and standard deviation (STD) of the 24 h clock offset overlap were decreased by 0.024 ns (8.6%) and 0.020 ns (13.1%), respectively; the fitting RMS of the 72 h clock offset of the quadratic polynomial was reduced for about 0.016 ns (6.5%).
From the comparison results, the precision of the orbit and clock offset with the estimated PCOs and PVs was improved. However, the improvements were mainly determined by the difference between the two models. Comparing the estimated PCOs and the MGEX PCOs in the x-offset component, the differences were 24 and 16 mm for CAST and the Shanghai Engineering Center for Microsatellites (SECM) satellites, respectively; for the y-offset, the difference was less than 9 mm; for the z-offset, the differences were 270 and 402 mm for CAST and SECM satellites, respectively. Considering that the BeiDou-3 Geostationary Earth Orbit (GEO) satellite ground station has too few observation data, the later launch of the BeiDou-3 Inclined Geosynchronous Orbit (IGSO) satellite results in a shorter data period, which is not enough to provide high-precision PCC estimations. These two types of satellites are not covered in this study.
Comparing the estimated PCOs of the Galileo Full Operational Capability (FOC) in-orbit [5] with the ground-based accurate phase center correction (PCC) model released by the European GNSS Agency (GSA), the differences between the two models were less than 5 mm in the x-offset and y-offset parameters, and the difference in the z-offset is about 152-304 mm. The similar difference in the z-offset of the BeiDou-3 MEO and Galileo might be caused by the unmodeled errors, and result in a scale error with the International Terrestrial Reference Frame 2014 (ITRF2014) [36]. This possible error might be the PCC model error at the receiver-end. Currently, IGS provides receiver antenna PCC correction for GPS and GLONASS frequencies only [4]. For BeiDou-2, BeiDou-3, and Galileo satellites, the GPS L1 and L2 model approximations are generally selected. Therefore, after the receiver antenna PCC calibrations of the Galileo and BeiDou frequencies are published, more accurate Galileo and BeiDou satellites antenna PCC parameters can be estimated in orbit.
Author Contributions: X.Y. and G.H. provided the initial idea for this study. X.Y. and L.W. designed the experiments. Z.Q. and S.X. performed the experiments. X.Y. analyzed the data and wrote the paper. G.H. and Q.Z. gave valuable advice on writing the paper.