Celestial frame tie from simulations using phase referencing to GNSS satellites

Aims. For decades now, researchers have been looking for a way to tie the kinematic and dynamic reference frames. Certain worldwide organizations have looked to using co-location in space, combining various techniques. Given the long list of possible applications of the Global Navigation Satellite System (GNSS), it is worthwhile investigating the connection between the most accurate and stable International Celestial Reference Frame (ICRF) and the Earth-centered Celestial Inertial reference frame (ECI) used in GNSS data processing. Methods. We simulated phase-referencing observations of GNSS satellites and nearby radio source calibrators to realize the connection between the two celestial reference frames. We designed two schemes for observation plans. One scheme is to select the satellite target when it can be observed by the greatest number of stations in order to obtain high-precision positioning. During each scan, we employ four regional networks to simultaneously track four chosen satellites. The alternative scheme is to observe satellite orbits of as many satellites as possible on different daily observations. In addition, to test the two schemes, we used Monte Carlo methods to generate 1000 groups of random errors in the simulation. Results. Finally, we estimate the right ascension and declination offsets ( ∆ α , ∆ δ ) of GNSS satellites in the ICRF, and then derive frame tie parameters based on those results: three global rotation angles ( A 1 , A 2 , A 3 ). The celestial angular offset results assessed from the former scheme show that this scheme leads to high precision of namely 1 mas, while the parameters of the frame tie determined from the second scheme can achieve an improved precision of better than 1.3 µ as.


Introduction
As humans are now contemplating the physics of deep space, observational targets correspondingly extend from the Earth and lunar satellites to the planets and their satellites, and from stars to extragalactic radio sources, especially with the development of the Global Navigation Satellite System (GNSS), satellite laser ranging (SLR), very long baseline interferometry (VLBI), and the Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS) system. It now therefore necessary to define a stable inertial celestial reference frame to describe the positions or motions of those observation targets. The celestial reference frames mainly contain the kinematic celestial reference frames and the dynamic celestial reference frames. The earliest kinematic celestial reference frames were the optical frames (e.g., FK4/FK5), which are determined by the positions of galactic stars. When VLBI became the most precise tool in astrometry in the fall of 1995, the International Celestial Reference Frame (ICRF) was produced by the International Astronomical Union (IAU) working group on reference frames, which is a no-net-rotation (NNR) kinematic quasi-inertial reference system based on the angular positions of extragalactic radio sources observed using VLBI. The dynamic planetary frames are defined by the equations of motion of the major celestial bodies in the Solar System. Many attempts were made in the 1980s and 1990s to understand the two types of celestial reference frames and their relationship. One of the aims of the HIPPARCOS mission was to link the optical catalog to the extragalactic radio frame with 1 mas positions and 1 mas yr −1 proper motions of 120 000 stars (Froeschle & Kovalevsky 1982;Feissel & Mignard 1998;Jacobs et al. 1998;Lestrade et al. 1995). Later on, the planetary and radio frames were tied by employing differential VLBI observations of various spacecraft and angularly nearby quasars (Newhall et al. 1986). In addition, several other techniques were used to determine the relative orientation of the planetary ephemerides and the radio frames. For instance, pulsar timing occultation of radio sources by planetary objects or spacecraft radio transmissions to obtain the positions relative to extragalactic radio sources were conducted and then compared to the corresponding coordinates determined from the spacecraft orbit about the target body (Jacobs et al. 1993). An extragalactic planetary frame tie was established using lunar laser ranging (LLR) and VLBI data analysis (Niell 1988;Finger & Folkner 1992;Folkner et al. 1994). The tie between the radio frame and other ephemerides was also obtained by conducting pulsar timing observations (Bartel et al. 1996).
The IAU encouraged the International Earth Rotation Service (IERS) to take appropriate measures in conjunction with the IAU working group on reference frames in order to maintain A6, page 1 of 9 the ICRF and its ties to reference frames at other wavelengths (Feissel & Mignard 1998). The second ICRF was established based on VLBI observations at S -and X-bands of 3414 extragalactic radio sources, which include 295 "defining" sources that determine the orientation of the frame axis (Fey et al. 2015). The noise floor of the ICRF2 is about 40 µas in position and 10 µas in the stability of the axis (Sotuela et al. 2012). The latest celestial reference frame ICRF3 was adopted by the IAU in 2018, taking the acceleration of the center of mass of the Solar System into account, including a total of 4588 sources in different bands (Charlot et al. 2020). Here, we aim to link the latest ICRF to the L-band celestial reference frame for GNSS processing (ECI).
Phase referencing is a method traditionally used in radio astronomy to measure or improve the position of weak radio sources or sources with uncertain positions by employing strong radio sources as calibrators. Most of the errors in the signal propagation and from the ground segments are eliminated by switching observations between the target and the calibrator with a small angular separation over a short time interval. Phase referencing plays a vital role in connecting reference frames (Lestrade et al. 1995) and has become an important observing method in deep space exploration, that is, Moon exploration and Mars exploration (Lanyi et al. 2005;Zhou et al. 2014;Haas et al. 2017). When the angular distance between the target and calibrator is smaller than the beam width of the VLBI telescope, two kinds of signals can be received at VLBI antennas simultaneously, which is a technique referred to as same beam interferometry (Gao et al. 2019). This method was applied to position Chang'E 3 and the lunar satellite SELENE. Furthermore, phase referencing synthetic imaging in astronomy has been successfully applied in the precise positioning of the Yutu patrol instrument of the Chang'E 3 project (Tong et al. 2014;Klopotek et al. 2019;Zhou et al. 2014Zhou et al. , 2015Li et al. 2013). With a wide range of applications, namely in geodesy, satellite altimetry, natural disaster monitoring and early warning, GNSS has grown into an indispensable precise positioning and navigation tool used daily in different fields. GNSS orbits are modeled based on daily arcs and are affected by different perturbations. Therefore, it is essential and promising to tie the celestial reference frame for GNSS data processing (ECI) to the most accurate and stable ICRF. The concept of tracking or observing GNSS satellites by VLBI can be traced back to the 1990s for calibrating the phase center of GPS satellites (Hase 1999).
From then on, VLBI telescopes were used to track GNSS satellites directly. The dynamic celestial reference frame for satellite data processing can be connected to the ICRF in this way. In addition, satellites such as APOD and E-GRASP were designed to carry various space geodetic tools on board. Simulations and measurements have been obtained for frame connections (Anderson et al. 2018;Sun et al. 2018). So far, several test experiments have been carried out, tracking GNSS signals directly (particularly GLONASS) without phase referencing, using a single VLBI baseline or a small network (Liu et al. 2014;Haas et al. 2014Haas et al. , 2017Kodet et al. 2014;Tornatore et al. 2011Tornatore et al. , 2014Plank et al. 2014Plank et al. , 2017. Nevertheless, these experiments mainly focused on tests on a small scale to check the feasibility of tracking Earth satellites with VLBI telescopes. Until now, phase-referencing experiments with VLBI stations to track GNSS satellites have only rarely been carried out on a large scale. Therefore, we propose the use of phase referencing on GNSS satellites to link the celestial reference frame for GNSS data processing to the ICRF. Furthermore, the satellite coordinates can be expressed in a quasi-inertial frame ICRF, which will improve the orbit and the Earth orientation parameters (EOPs).
However, we still face many difficulties when using the above technique, such as those caused by limitations in regards to the observation bands of VLBI receivers and the fast motion of satellites compared to the slow slew speed of antennas. It is therefore necessary to upgrade the software and hardware tools to realize VLBI observations of GNSS satellites. First, the traditional VLBI antenna equipment is used to observe distant objects in deep space, that is, extragalactic radio sources. Due to the great distance of the objects, the signals arriving at the stations are regarded as parallel. Therefore, the vector direction of the radio sources on each station is identical. GNSS satellites are close to the Earth, the signals arrive at each station in the form of spherical waves. Accordingly, the vector directions of a satellite seen from different antennas are different. Second, a high slew speed is required for tracking fast-moving satellites. The VLBI stations currently used worldwide are still relatively old, large, and slow, except for the newly built small VGOS antennas. In addition, extragalactic radio sources are usually observed in the S and X bands with wide bandwidth and weak signals. In contrast, GNSS satellites transmit L-band signals with narrow bandwidth and high strength. Hence, a long integration time is required to extract the observations of the extragalactic radio source in the L-band, and a signal amplifier is essential. It is necessary to change the receiver for the satellite observations, and signal attenuation should to be implemented in order to prevent the strong satellite signals from damaging the receiver ). Moreover, a theoretical VLBI delay model for GNSS observations needs to be introduced on account of the different observation modes of VLBI and GNSS (Sekido & Fukushima 2006). In addition, the calculation of relativistic effects for extragalactic radio sources and GNSS satellites is different. Finally, systematic errors also vary from one observational target to another. In the data processing for the extragalactic radio sources, one should adequately analyze the radio source structure (Xu et al. 2019), clock offsets, clock jumps, and phase ambiguity. But in the case of GNSS satellites, phase-center offsets, solar radiation pressure, and other models need to be considered.
In this work, we schedule and simulate phase-referencing observations of GNSS satellites in order to determine the celestial angular position offsets of GNSS satellites in the ICRF. Furthermore, we establish a link between the ECI and the ICRF. As one of the GNSS systems, GPS is used as an example. For comparison, we designed two schemes to perform the related simulations and estimations. In the first scheme, we select the satellite with the highest common visibility as the observation target. Four regional VLBI networks from Asia, Europe, Australia, and America are employed to simultaneously observe four different satellites and their corresponding angularly nearby calibrators. In the second scheme, observations of each GPS satellite are simulated for a whole day and are accumulated for days to get a better geometric distribution. Based on the simulated observations in the two schemes, the celestial angular position offsets are estimated. The rotation angles are derived as frame tie parameters through Monte Carlo simulations. The results from all regional networks or from all GPS satellites are then combined to derive a more feasible frame tie with both schemes, respectively.
A6, page 2 of 9 Fig. 1. Geometrical relationship diagram. R is the radius of the Earth, H is the satellite altitude, b is the baseline, ρ is the distance between the satellite and the station, γ is the angle between the signal path and the radius, ε is the elevation, θ is half of the largest opening angle, and β is the angle between R and the geocentric vector of the satellite.

Scheme one, based on the simultaneous maximal common visibility of GNSS satellites
In order to derive a frame tie from phase referencing of GNSS satellites using VLBI telescopes, scheduling and simulations are essential to verify the feasibility, to attempt to solve problems of systematic errors, to find the optimum approach for the experimental design, and to explore the best algorithm for the parameter estimation. The GPS, for example, has an altitude of approximately 20 000 km. The relations are expressed in Fig. 1. First, we can calculate β with β = arcsin b/2 R . Then we can derive ρ with the law of cosines The opening angle is the angle of a satellite to points st1 and st2, which are theoretical points indicating the maximal baseline on which the satellite is still visible to two stations on the ground; it is at its maximum when the altitude is fixed and the satellite passes the perpendicular bisector of the baseline. The value of half of the largest opening angle can be obtained with the rule of sines R sin θ = ρ sin β . Simultaneously, the angle γ is calculated using the equation R+H sin γ = R sin θ . The maximal elevation angle is then obtained with the relationship ε + γ = π/2 (Plank 2013).
As the frame tie at each epoch is derived based on at least three points to complete the transformation, we designed four regional station networks to track four satellites at the same time. When the baselines are short, the opening angle of a satellite to two stations is small as the satellite altitude is unchanged. At the same time, the elevation angles of the satellite or nearby radio source calibrators are high, as shown as Fig. 2. For example, when the baselines are 2000 km, the elevation angles are approximately 80 degrees and the opening angle is about 5 degrees. In this case, the geometrical distribution of the individual satellite is poor for high-precision positioning. To compensate for this insufficiency, we consider using observations from baselines that are longer than 2000 km. Also, we select satellites based mainly on the simultaneous common visibility (choosing the satellite that can be observed by the largest number of VLBI stations simultaneously) under the condition that the angular separation (calculated with the right ascension and declination) between the satellite target and the radio source calibrator is less than 5 degrees.
In contrast, high-quality calibrators are chosen according to the following criteria: quasars with flux density above 0.2 Janskys; signal-to-noise ratio (S/N) above 20; spherical angular distance between the target and calibrator of less than 5 • ; and position uncertainty (right ascension and declination separately) in the ICRF3 of the radio source of better than 0.5 mas.
Extragalactic radio sources at different wavelengths have different flux densities. Information about the flux density in the S and X bands is provided on the NRAO website 1 . Assuming that the VLBI stations can receive signals from satellites and radio source calibrators in the L band, the flux density of quasars in the L band can be calculated using the two-point spectral index with the information in the S and X bands (Shabala et al. 2014). Currently, the slew speed of available VLBI antennas is about 0.4 • −3 • per second. The slew speed of VLBI2010 antennas lies at 6 • −12 • per second (Schuh & Behrend 2012). Within 12.5 s, the slowest antenna can realize nodding observations between targets and calibrators if the angular separation is less than 5 • . In addition to the duration of slewing of the antennas to the different targets, the minimal scan length sc needs to be considered, which can be computed using Eq. (1): where S /N is the signal-to-noise ratio as introduced above; SEFD 1 and SEFD 2 denote the system equivalent flux densities for two stations of a baseline; Samprate is the unified sampling rate for the whole array; N chan represents the channel number; and Bits is the bandwidth. For instance, with one-bit sampling and the standard S /N, flux density, and SEFD, the scan length can reach 10 s with the European network. The VGOS antennas are small with high slew speed, meaning that the switching time is short.

Scheme two: based on a long-term optimal geometric distribution of GNSS satellites in ICRF
In order to obtain the frame tie with the highest precision by improving the geometric distribution of satellites in ICRF, we designed an alternative scheme with the goal of observing orbital arcs with maximal lengths for each satellite on each day. We utilize 1000 Monte Carlo simulations to get high feasibility. A long-term (at least one month) simulation is carried out in order to observe as many GPS satellites as possible. By tracking different satellites on different days, we can obtain the geometric distribution of the entire GPS constellation over 32 days when all 32 satellites of the GPS constellation are available. In addition, we use three months of simulations for further validation. Consequently, the celestial frame tie parameters are obtained not only for individual satellites but for the entire GPS constellation by accumulating observations over a longer period.

Scheduling method
Using the two-line elements (TLEs) 2 from the North American Aerospace Defense Command (NORAD) CelesTrak, the satellite orbits are integrated and provide the input for the VLBI antenna pointing. The accuracy of the satellite positions derived from the TLE is better than 10 km over a five-day period, which corresponds to an angular position of approximately 100 arcsec. As the beam width is a few arcminutes, it is sufficient for antenna pointing. However, observing sources at the edge of the beam might effect the observable. We therefore use a more precise orbit from the IGS products in our study to schedule the phase-referencing experiment. According to our study, the "stopand-go" mode will lower the quality of the observation. The ideal observation mode is to move the antenna during the observation of the satellite. As the satellite has a high S/N, the scan length is short, and so we can use VGOS stations (assuming they are capable of observing GNSS satellites) to observe the satellite for the short scan time (calculated using Eq. (1)). The IGS website 3 provides short-term predicted orbits at the level of a few centimeters Noll (2010). We use the predicted orbits to schedule the phase referencing of the GPS satellites for a full day. We selected calibrators from the ICRF3-S/X catalog. As the duration of the "calibrator-target-calibrator" in the phase-referencing procedure depends on the performance of the VLBI telescope (e.g., slew speed), angular separation, and scan length -which is less than 2 min -, we determine the interval of the phase-referencing observations to be 5 min ("quasar-satellite-quasar" mode; 2 min for quasar, assuming that VGOS stations can track the satellite for the scan length time for the integration and then switch to the quasar within 1 min in total).

Simulation models
Geometric models and main correction models should be considered in the simulation. However, errors from signal propagation and ground instruments in phase referencing are mostly canceled out by differential delays between satellites and calibrators. We can calculate the ionospheric delay errors (depending on frequency) of phase referencing observations based on the ionospheric products derived from GNSS observations and VLBI observations with an accuracy of 2.2 TECU (Maennel & Rothacher 2016) or based on a TEC map with a precision of 3.7-3.9 TECU (Sekido et al. 2003). In our work, only the main errors of phase-referencing observations are discussed, including differences in clock offsets, atmospheric delays, relativistic effects between the object and calibrator, and errors caused by the satellite orbit and position uncertainties of radio sources and white noise, as presented in Eq. (2): where ϕ geo is the theoretical phase delay; ϕ clk are clock offsets; ϕ atm represents the effects of different neutral atmospheric delays; ϕ rel is the relativistic effects; ϕ orb is the delay caused by errors in satellite orbit; ϕ sou is the delay from quasar position uncertainties; and ϕ wn is white noise. Among those errors, we use Allan standard deviation (ASD) estimation to simulate the characteristics of clocks. Clock offsets are generated using random walk and integrated random walk models with the parameter AS D = 1 × 10 −14 @50 min (Herring et al. 1990). We simulate zenith wet delays with turbulence models. The turbulence theory (Treuhaft & Lanyi 1987) is applied with a dedicated strategy proposed by Nilsson (Nilsson et al. 2014), which includes the following parameters: site-dependent structure constants Cn, effective heights h, wind velocity components (Ve and Vn), initial zenith wet delay (zwd0), correlation interval (dhseg), and a height increment for the numerical integration (δh). The value of Cn is found empirically for each antenna, except for VGOS stations, for which we use 2.5 × 10 −7 . Meanwhile, the other parameters are set to a standard value. We introduce the Vienna mapping function (Boehm et al. 2006) for the tropospheric delays of satellites and radio source calibrators in the simulation. We set the white noise of phase referencing for satellites and quasars to 4 ps, assuming that the VGOS requirements (measurement error of 4 ps) are met. The orbit errors for satellites are set to 2 cm, and the position offsets of radio sources are 0.5 mas according to the data-selection criteria.

Estimation methodology
The approach for estimating the reference frame tie is the determination of rotation angles between the celestial reference frame for GNSS data processing (ECI) and the ICRF. The transformation between the two celestial reference frames is expressed as the following matrix: (3) Equation (4) (Lestrade et al. 1995) describes the relationship between rotation angles and celestial angular offsets (offsets of right ascension α and declination δ) of satellites: The angular positions (α and δ) of satellites are transferred to the celestial reference frame (ECI) using EOPs from IGS website. α 1 and δ 1 are the original right ascension and declination of GPS satellites in the ECI and α 2 and δ 2 are corrected celestial angular position offsets of GPS satellites in ICRF3. A 1 , A 2 , and A6, page 4 of 9 L. Liu et al.: Celestial frame tie from simulations using phase referencing to GNSS satellites Fig. 3. Regional VLBI station networks in Asia, Europe, Australia, and America indicated with red stars, blue squares, green dots, and purple triangles, respectively, with station names denoted in the same colors.
A 3 are rotation angles between the celestial reference frame of the ECI and the ICRF3.
The first step is to determine the celestial angular offsets of satellites between two sets of celestial reference frames. Satellite positions from the IGS SP3 product, station coordinates from ITRF2008, and EOPs from IGS products and IERS are fixed in our work. The parameters that we estimate are the angular positions and offsets of the satellite in the ICRF with the observation model in Eq. (2). In this study, we introduce the advanced VLBI delay model for GPS satellites (Sekido & Fukushima 2006). For Earth-based VLBI observations, this model claims to achieve the precision of one picosecond for sources with an altitude of 100 km above the Earth's surface (Sekido & Fukushima 2006). Additionally, the VLBI delay model from IERS conventions 2010 4 is adopted for AGN calibrators. The partial derivatives of phase delay with respect to satellite position dτ dr are then derived by applying the related model. Based on these latter, we estimate instantaneous state vectors at each epoch using the singular value decomposition (SVD) least squares method. Celestial position offsets defined as offsets of angular position describe the shifts of right ascension (α) and declination (δ) from their prior positions, which are calculated according to the relation between three-dimensional Cartesian coordinates and spherical coordinates. Furthermore, we deduce partial derivatives of the angular positions with respect to three rotation angles (A 1 , A 2 , A 3 ) according to Eq. (4). After these steps, rotation angles are obtained as frame tie parameters.

Simulation results based on scheme one
As designed in scheme one, we use the four aforementioned regional station networks (shown in Fig. 3) to track four selected satellites simultaneously. Figure 4 shows the distribution of satellites and calibrators in the ICRF over 1 day, where satellites and sources are chosen according to the set of criteria used for the first scheme outlined above. In this figure, the Australian network shows an even coverage in the southern celestial hemisphere, and other networks provide coverage mainly in the northern hemisphere. Figure 5 shows the baseline lengths. American and European networks have more than 50 baselines including many short ones and other networks have less baselines. Figure 6 presents the number of observations from the Asian, European, Australian, and American networks for each scan during May 1, 2021, which are calculated in the simulation process. The strong correlation between the number of observations and the number of baselines is revealed and shown in Figs. 5 and 6. Networks with more baselines provide more common visibility and result in more observations.
Phase delay residuals are computed as simulated observations minus the theoretical value (O-C). The majority of the total residuals come from the simulated main stochastic errors stated above. Overall, most residuals are better than ±1 ns. Figure 7 shows that the estimated celestial angular offsets fluctuate by several milliarcseconds. The results are related to the baseline lengths, the number of observations, and the error sources added in the simulated observations. The four regional networks provide the precision better than 5 mas. The significant amounts of bias in the results are from the simulated errors. When the baseline lengths are smaller than 2000 km, the precision on positioning is low. The American network has the highest number of observations, and therefore provides best precision after short baselines are excluded. Asian networks have a smaller number of observations, and therefore the angular offset results are affected by significant amounts of bias.
We carried out 1000 Monte-Carlo simulations to improve the feasibility. Our results reveal that the derived rotation angles vary with simulated errors, and the formal errors from the same network are at the same level. Table 1 presents the average rotation angles from 1000 Monte Carlo simulations accompanied with their formal errors. The celestial rotation angle results from each regional network are derived with a precision of several to A6, page 5 of 9 A&A 671, A6 (2023) Fig. 4. Selected GPS targets and angularly nearby calibrators with four regional VLBI networks. Satellite targets in dots and angular nearby radio source calibrators in stars observed by networks from Asia, Europe, Australia, and America are described in red, blue, green, and purple.   7. Celestial angular offsets from four regional networks. Right ascension and declination offsets are shown in the upper and lower panels, respectively. Regional VLBI networks from Asia, Europe, Australia, and America are marked separately in red stars, blue squares, green triangles, and purple dots in both panels. hundreds of microarcseconds. The highest precision of frame tie parameters is estimated based on data from the American network. A 1 , A 2 , and A 3 angles determined from all regional networks are (7.8 ± 2.1 µas,-37.5 ± 3.5 µas, −66.7 ± 8.3 µas), and the precision of rotation angles is at a level of 1 µas. right ascension (mas) time (hours on 2021-5-1) network, including observations from stations of all four regional networks in scheme one. We can obtain this geometric distribution for the entire GPS constellation when a different satellite is successively monitored for 32 days, because GPS satellites repeat nearly the same ground track twice per day (32 GPS satellites are available). For the sake of clarity, the 32 days of observations are merged into one graph. In the second scheme, we simulate phase referencing observations for all available GPS satellites and source calibrators for a long period. The right ascension and declination offsets in the ICRF are estimated, and the results are shown in Fig. 9. In this plot, the precision of the declination offsets is generally better than 10 mas, which is slightly worse than the accuracy of right ascension offsets. The second scheme produces poorer angular position results for each satellite than the first scheme, which is probably attributable to the fact that, in the second scheme, the satellite targets are not selected based on their common visibility. The precision of the angular positions is lowered by the lower number of observations per scan. Frame-tie parameters for each satellite are estimated using Eq. (4) and averaged from 1000 Monte Carlo simulations. Figure 10 shows rotation angles determined from each satellite. In the figure, rotation angles are between ±0.5 mas, and the uncertainties are less than 180 µas. These results are still poorer than those in the first scheme and the reason is because the distribution of the single satellite is not as good as the four optimally selected satellites in the first scheme. This is clearly evident when comparing the track of the individual satellites in Fig. 8 with all the observed satellites in Fig. 4. Additionally, three-month simulations are carried out to verify the stability. Celestial frame tie parameters for each satellite deduced from long-term continuous simulations are demonstrated in Fig. 11. The rotation angles are between ±0.6 mas, and the precision is better than 140 µas. Angular offset results from each satellite over a three-month period are marginally superior to those from 1000 Monte Carlo simulations for each satellite, but are at approximately the same level.  Rotation angles with results from all available GPS satellites for one month (May, June, or July) and three months together (All). A 1 , A 2 , and A 3 are in red, blue, and green bars. Figure 12 shows derived frame-tie parameters from individual month (May, June, and July) and from all three months ("All"). This figure demonstrates that the precision of rotation angles from a one-month simulation is better than approximately 3.5 µas. Frame-tie parameters from three months are (A 1 = −13.1 ± 0.8 µas, A 2 = 6.8 ± 1.2 µas, A 3 = 20.4 ± 1.3 µas), indicating that rotation angles can be obtained from scheme two at a level of 1 microarcsec. These frame-tie results are slightly improved thanks to the better geometric distribution of the GPS satellites in the ICRF in the long term.

Conclusion
We simulate phase-referencing observations of GPS satellites by VLBI telescopes to determine frame-tie parameters for linking the dynamic satellite reference frame to the ICRF. In this study, we design related experiments with two schemes to schedule and simulate phase referencing observations of GPS satellites and angularly nearby radio source calibrators. In the first scheme, we estimate celestial angular offsets of GPS satellites in the ICRF at 1 mas level. After 1000 Monte Carlo simulations, we derive three rotation angles as frame-tie parameters individually from each regional network with submilliarcsecond precision. When all observations from four regional networks are combined, we acquire frame-tie parameters with a precision of better than 10 µas. In the second scheme, we compute the angular offsets of individual GPS satellites in the ICRF with a precision of a few milliarcseconds. In addition, we obtain celestial frame-tie parameters at a level of 0.1 mas for each satellite by performing 1000 Monte Carlo simulations. Frame-tie parameters from three-month cumulative observations of individual satellites are marginally superior to those from 1000 simulations. Combining simulations for one month and three months as designed in scheme two, we get more precise frame-tie parameters than those in the first scheme thanks to a better geometric distribution over the long term. Finally, the three-month simulation based on the second scheme can obtain a precision of better than 1.3 µas, the highest precision that we achieve. In summary, both of the designed schemes can achieve a precision of 1 microarcsecond. Furthermore, we recommend scheme one to the geodetic or astronomical community for short-term (i.e., one day) or even epoch-based frame-tie derivation. For now, the VLBI 24 h sessions for GPS observations over long timescales are not realistic because the standard VLBI observations have been scheduled during this period. However, this scheme can be used to obtain a high-precision and high-stability frame tie from long-term (i.e., one month) observations when more VGOS stations are established. In addition, we can get a good geometric distribution with six GPS orbit planes through a six-day observation by choosing one satellite from each orbital plane on each day. Furthermore, problems due to systematic uncertainties in individual satellites or misalignment of the whole GPS constellation might be solved in the future.