Quality Assessment of GNSS, SLR, VLBI, and Doris Inputs for ITRF2014 and ITRF2020 Using Trf Stacking Methods

1. Introduction

A unique standard terrestrial reference frame (TRF) is the fundamental to ensure inter-operability and consistency of geodetic products and to adequately exploit various measurements collected by ground-based sensors, or via artificial satellites for Earth science and geodesy applications [1,2,3], including plate tectonics [4,5], co-seismic and post-seismic deformations [6,7], global geophysical fluid dynamics [8,9], ice melting [10], accurately determine point positions, and the rate of sea level rise [11,12], etc. Generally, the most widely utilized reference framework is the International Terrestrial Reference Frame (ITRF) series [13], which has a history spanning four decades. Since its inception in 1988, the ITRF has undergone 14 versions, namely ITRF 88, 89, 90, 91, 92, 93, 94, 96, 97, 00, 05, 08, 14, and 20.

The current ITRF2020 encompasses the positions, velocities, and annual and semi-annual signals of globally distributed reference stations, integrating the strengths of the four Space Geodesy technologies: GNSS, SLR, VLBI, and DORIS [2]. Typically, analysis centers (AC) for GNSS, SLR, VLBI, and DORIS reprocess historical observations employing the latest optimal computation strategies and models. Then, the results from these centers are combined by the combination centers (CC) for the four technologies, IGS, IVS, ILRS, and IDS. These integrated outcomes are considered to be relatively high-quality and accurate solutions because of the latest processing strategies and models, and the correction of system errors of different analysis centers, and submitted as the input data for establishing and updating ITRF.

However, when utilizing these high-quality solutions, one must still consider outliers, inaccurate covariance, and datum precision evaluation. While SLR employs a loosely constrained solution, VLBI utilizes a normal equation, and DORIS derives a technical solution from internal constraints. Consequently, the solutions of these three technologies are not aligned to a fixed frame and traditional time series analysis methods prove insufficient in assessing the quality of these solutions. Additionally, traditional time series analysis methods cannot assess the datum information for all of the four techniques either. Therefore, we propose a TRF stacking approach following Altamimi, et al. [14] to achieve the previously mentioned objectives. TRF stacking accumulates epoch solutions to long-term TRF expressed by coordinates and velocities in a specific epoch, enabling us to evaluate observational qualities and the datum precision. Furthermore, to establish 2014 and 2020 ITRF, solutions of each Space Geodesy technique incorporated new strategies and models, we can investigate the impact of these strategies at the stacking TRF level.

In general, the solutions for each technique in 2020 exhibit significant improvements compared to the previous versions. For GNSS, the new IERS secular pole model was adopted in Repro3, replacing the previous mean pole model used in Repro2 from the original IERS 2010 Conventions. Consequently, Repro3 and IGS14 exhibit systematic biases, making alignment impractical. Therefore, a special reference frame, IGSR3, was utilized, aligning with the new IERS secular pole model, the igsR3.atx satellite PCOs, and the igsR3.atx ground antenna calibrations [15].

For SLR, a series of improvement measures were implemented to address its systematic biases [16], like enhancements to the satellite center of mass model corrections [17]. A comprehensive investigation into the site-specific systematic errors in the data collected since 1983 was conducted [18]. Regardless of reported issues or configuration changes at the sites, both the coordinates and systematic biases were simultaneously estimated for all of the 7 ACs and CCs (ASI, BKG, DGFI, ESA, GFZ, JCET, NSGF) [19].

For VLBI, some new and improved models for the IVS analysis have been applied like a constant drift of 5.8 µas/year on Galactic aberration. Additionally, mid-session-epoch solutions are calculated to present the best performance of VLBI, but they are transformed to 12-hour to ensure the consistency with the epochs of other techniques which are submitted for establishing the ITRF. [20,21].

When combining results from various analysis centers, IDS avoids introducing external reference frames and employs internal constraints during technical within-group combinations [22]. Consequently, the geocenter and scale results from DORIS observations directly reflect its precision in determining geocenter and scale. Compared to 2014, DORIS 2020 solution has undergone some new improvements: the application of new phase center variations for Alcatel DORIS ground antennas and the gradual replacement of Alcatel with Starec antennas over time; better modeling of the surface forces on the satellites [22,23].

The remaining sections of this paper are organized as follows: Section 2 presents the methods employed in this study, Section 3 introduces the input data and the results, finally, Section 4 provides the conclusions.

2. Methods

2.1. Stacking Model

A significant application of space geodetic observations is the continuous determination of ground station coordinates [24]. Using four technologies, GNSS, SLR, VLBI, and DORIS, coordinates sets and their internally consistent Earth Orientation Parameters (EOP) can be estimated based on observations collected at each measurement interval, daily and weekly [25]. Due to the complex surface deformation characteristics of the Earth, the coordinates of individual epochs not only include Earth’s crustal movements, terrain variations, and surface loads as geophysical signals but also involve various random errors. Consequently, in order to mitigate the influence of random factors and better apply physical signals, it is common to stacking observations from individual epochs over a certain time span to coordinates, velocities, and other modelable parameters (such as periodic signals) at a specific reference epoch. This process yields a reference framework through the stacking of observations, and we call it TRF Stacking method.

The transformation relationship between TRFs is typically represented using the Helmert similarity transformation model. For a station i of technique s, located in a specific epoch $t_s^i$, the coordinates ${\mathbf{X}}_s^i$ in the epoch reference frame are related to the stacking reference frame coordinates ${\mathbf{X}}_{ck}^i$ at the reference epoch $t_0$, as well as the Earth Orientation Parameters (EOP) in the two reference frames, according to the following transformation:

$${\mathbf{X}}_s^i={\mathbf{X}}_{ck}^i+(t_s^i-t_0){\dot{\mathbf{X}}}_{ck}^i+{\mathbf{T}}_s+{\mathbf{D}}_s{\mathbf{X}}_{ck}^i+{{\mathbf{R}}_s\mathbf{X}}_{ck}^i$$

(2.1)

$$\left\{\begin{array}{c}{\mathcal{x}}_s^p={\mathcal{x}}_c^p+{R2}_s\\ {\mathcal{y}}_s^p={\mathcal{y}}_c^p+{R1}_s\\ \begin{array}{c}{UT}_s={UT}_c-\frac{1}{f}{R3}_s\\ {\dot{\mathcal{x}}}_s^p={\dot{\mathcal{x}}}_c^p\\ \begin{array}{c}{\dot{\mathcal{y}}}_s^p={\dot{\mathcal{y}}}_c^p\\ {LOD}_s={LOD}_c\end{array}\end{array}\end{array}\right.$$

(2.2)

where $f=$1.002737909350795, ${\mathbf{T}}_s$ represents the three translation parameters, ${\mathbf{D}}_s$ denotes scale parameters, ${\mathbf{R}}_s$ signifies rotation parameters and ${\mathbf{R}}_s=\left[\begin{array}{ccc}0& -{R3}_s& {R2}_s\\ {R3}_s& 0& -{R1}_s\\ -{R2}_s& {R1}_s& 0\end{array}\right]$. If there is a need to model additional signals such as periodic terms for the stations, further model optimization can be performed in Equation 2.1. Corresponding adjustments should be made in the subsequent datum definition part [2]. However, the primary objective of using single- technique TRF stacking in this paper is not to analyze nonlinear features. Therefore, only linear terms are included.

2.2. Datum Definition

According to Eq. 2.1, the introduction of transformation parameters and the estimation of velocity terms result in a rank deficiency in the normal equations, with a rank deficiency of 14 [26]. It can also be considered as the datum for the stacking TRF, ${\mathbf{X}}_{ck}^i$ and ${\dot{\mathbf{X}}}_{ck}^i$, is undefined. The rank deficiency can to be addressed by datum defining.

According to the definition by the International Earth Rotation and Reference Systems Service (IERS), the datum for the International Terrestrial Reference System (ITRS) includes the origin, scale, orientation, and their variations over time [27].

Due to the distinct characteristics of various space geodetic techniques, such as SLR, GNSS, DORIS, which are sensitive to the Earth’s center of mass, the TRF origin (CM) can be defined [28,29,30].

The scale, by definition, should be the scale of a local TRF in the sense of general relativity. All of the four techniques are sensitive to the scale, and in practical applications, physical parameters (such as Earth’s gravitational constant GM and the speed of light c) and a certain relativistic correction model are employed to determine it. If different techniques use distinct physical parameter values and models, their scales may differ. However, adopting recommended values from the IERS specifications in data processing tends to minimize scale differences determined by various techniques.

The orientation of the current ITRS has evolved from the orientation defined by BIH1984.0, and new versions of ITRS maintain consistency with the previous version’s orientation by conforming to the no-net-rotation (NNR) condition, which ensures horizontal motion relative to the Earth’s surface [31].

In mathematics, the definition of a datum is typically realized through the imposition of constraints. Generally, two constraint methods are employed: when defining the datum using a set of external reference coordinates and their velocities, denoted as ${\mathbf{X}}_R$ and ${\dot{\mathbf{X}}}_R$, the minimum constraint approach is utilized; when preserving the intrinsic physical parameter information associated with ${\mathbf{X}}_c$, the minimum constraint is applied.

Minimum constraint means the minimization of Helmert transformation parameters between external reference frames, ${\mathbf{X}}_R$, ${\dot{\mathbf{X}}}_R$, and the Stacking TRF, ${\mathbf{X}}_c$, ${\dot{\mathbf{X}}}_c$ [32]. The relationship of ${\mathbf{X}}_R$, ${\dot{\mathbf{X}}}_R$ and ${\mathbf{X}}_c$, ${\dot{\mathbf{X}}}_c$ can be defined as

$$\left[\begin{array}{c}{\mathbf{X}}_R\\ {\dot{\mathbf{X}}}_R\end{array}\right]-\left[\begin{array}{c}{\mathbf{X}}_c\\ {\dot{\mathbf{X}}}_c\end{array}\right]=\left[\begin{array}{c}\begin{array}{cc}\mathbf{A}& 0\end{array}\\ \begin{array}{cc}0& \mathbf{A}\end{array}\end{array}\right]\left[\begin{array}{c}\mathsf{\theta }\\ \dot{\mathsf{\theta }}\end{array}\right]$$

(2.3)

$\left[\begin{array}{c}\mathsf{\theta }\\ \dot{\mathsf{\theta }}\end{array}\right]$ encompasses the seven transformation parameters and their rates, the matrix A includes the partial derivatives of the transformation parameters, and the design matrix A is defined as

$$\mathbf{A}=\left[\begin{array}{ccccccc}1& 0& 0& x_c^i& 0& z_c^i& -y_c^i\\ 0& 1& 0& y_c^i& -z_c^i& 0& x_c^i\\ 0& 0& 1& z_c^i& y_c^i& -x_c^i& 0\end{array}\right]$$

(2.4)

Expressing Equation 2.3 as ${\mathbf{X}}_1-{\mathbf{X}}_0=\mathbf{G}\mathsf{\vartheta }$, where $\mathsf{\vartheta }$ is the least squares solution and given by

$$\mathsf{\vartheta }={\left({\mathbf{G}}^{\mathrm{T}}\mathbf{G}\right)}^{-1}{\mathbf{G}}^{\mathrm{T}}\left({\mathbf{X}}_1-{\mathbf{X}}_0\right)$$

(2.5)

Letting ${\left({\mathbf{G}}^{\mathrm{T}}\mathbf{G}\right)}^{-1}{\mathbf{G}}^{\mathrm{T}}=\mathbf{B}$, considering $\mathsf{\vartheta }$ as the observed quantity with a weight ${\mathbf{P}}_{\mathsf{\vartheta }}={\left({\mathbf{Q}}_{\mathsf{\vartheta }}\right)}^{-1}$, and treating ${\mathbf{X}}_0$ as the parameter, the minimization condition $\mathsf{\vartheta }=0$ leads to the normal equation for Equation 2.5:

$${\mathbf{B}}^{\mathrm{T}}{\mathbf{P}}_{\mathsf{\vartheta }}\mathbf{B}\mathsf{\delta }{\widehat{\mathbf{X}}}_0={\mathbf{B}}^{\mathrm{T}}{\mathbf{P}}_{\mathsf{\vartheta }}\mathbf{B}({\mathbf{X}}_1-{\mathbf{X}}_0^{\left(0\right)})$$

(2.6)

$\mathsf{\delta }{\widehat{\mathbf{X}}}_0$ represents the correction to the parameters ${\mathbf{X}}_0$.

Internal constraint signifies that for all transformation parameters within the span, they are set to zero at the reference epoch $t_0$, including both the values and velocities.

$$\left\{\begin{array}{c}{\mathsf{\theta }}_k\left(t_0\right)=0\\ {\dot{\mathsf{\theta }}}_k=0\end{array}\right.$$

(2.7)

Following the derivation, we obtain:

$$\left\{\begin{array}{c}{\displaystyle \sum _{k\in K}}{\mathsf{\theta }}_k=0\\ {\displaystyle \sum _{k\in K}}\frac{{\mathsf{\theta }}_k}{{(t_k-t_0)}^{-1}}=0\end{array}\right.$$

(2.8)

2.3. Outliers and Variance Component Estimation

Based on the posterior residuals ${\mathbf{v}}_s$, the posterior variance factor can be computed by

$${\sigma }_s^2=\frac{{\mathbf{v}}_s^T{\mathbf{P}}_s{\mathbf{v}}_s}{f_s}$$

(2.9)

where ${\mathbf{P}}_s$ is the weight matrix of the observations (obtained by taking the inverse of the covariance matrix), and $f_s$ is

$$f_s=6n_s-\mathrm{t}\mathrm{r}\left({\mathbf{A}}_s{\mathbf{N}}^{-1}{\mathbf{A}}_s^T{\mathbf{P}}_s\right)$$

(2.10)

This approach can also be referred to Davies and Blewitt [33]. It can serve as a residual-based posterior weighting method until ${\sigma }_s^2$ approaches unity. Observations with normalized residuals (residual divide standard deviation) exceeding 4 are considered outliers and are subsequently removed during this rescale process.

Weighted Root Mean Square (WRMS) is computed by

$$\mathrm{W}\mathrm{R}\mathrm{M}\mathrm{S}=\frac{{\mathbf{v}}_s^T{\mathbf{P}}_s{\mathbf{v}}_s}{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}\left({\mathbf{P}}_s\right)}$$

(2.11)

3. Results and Analysis

3.1. Data Introduction and Preprocessing

The selected data for this study is presented in the Table 1. To ensure a stable estimation of velocities and prevent situations where only one observation exists in a segment, making velocity unsolvable, stations with an observation count greater than 50 are chosen. For VLBI, we opted to use data from a minimum of three stations per session, starting from the 600th GPS week onwards. In 2014, 2887 sessions are chosen, with 20 sessions excluded due to abnormal transformation parameters. In 2020, 4433 sessions are selected, and 11 sessions with abnormal transformation parameters were excluded. Mid-session epoch solutions are used rather than 12-hour epoch solutions. Regarding SLR, the 2014 solutions had fixed weekly solutions from the 678th GPS week onwards, so both the 2020 and 2014 SLR datasets are selected from weeks beyond the 678th GPS week. For GNSS, weekly solutions rather than daily solutions are utilized here. All weekly solutions of IDS-09 for 2014 and IDS-16 for 2020 are selected for DORIS. Core stations are selected referred to ITRF2014 and ITRF2020 core-networks, for DORIS 2014 we select more core stations, because there are only two DORIS stations in ITRF2014 core-network.

When conducting datum definition or orientation, the external reference frame coordinates of a set of core stations from each technique are employed, station and core station distributions are illustrated in Figure 1. Generally, the utilized station numbers of GNSS, SLR, VLBI, DORIS are 911, 56, 56, 146 (1656, 33, 30, 57) for 2014 and 1181, 63, 63, 175 (2773, 34, 25, 80) for 2020, where the numbers in parentheses indicate the count of segments due to discontinuities. We can observe that there is an abundance of GNSS stations, mainly concentrated in terrestrial regions, VLBI and SLR have a relatively limited number of stations, and the north-south distribution is uneven for all three techniques. DORIS stations, on the other hand, exhibit a more uniform north-south distribution, and often exists multiple stations within one site.

SLR inputs consist of solutions with loose constraints, and the substantial formal errors documented in the SINEX files result from the relaxed positional constraints implemented during laser range data reductions. These errors are indicative of rotational noise associated with the reference frame. To mitigate this effect, the covariance matrices of the solutions are adjusted by imposing the no-net-rotation (NNR) condition. No no-net-scale (NNS) condition or no-net-translation (NNT) is applied to preserve the intrinsic SLR origin and scale [34,35,36].

The NEQs for each VLBI session are inverted by incorporating NNT and NNR conditions into the a priori station position coordinates provided in the input SINEX files. To maintain the inherent scale of VLBI, no NNS condition is applied during the inversion of the NEQs.

The post-seismic deformations for all four techniques are initially subtracted based on the analysis results from the ITRF. The segmentation of the time series is also referred to the results obtained from the ITRF. Corresponding file download links can be accessed at the following address: https://itrf.ign.fr/.

3.2. Stacking Results

We perform outlier removal of the coordinate observations and posterior scaling of the covariance matrix through iterative Space Geodesy TRF stacking. Observations with normalized residuals exceeding 4 are considered outliers and are subsequently removed. This iterative process continues until no outliers can be detected and the variance scale factor approaches 1. The iteration results are tabulated in Table 2.

Two datum definition strategies are utilized for TRF stacking. Initially, to understand the characteristics of the origin and scale determination for GNSS, SLR, and DORIS technologies, we applied internal constraints to their translation and scale parameters to retain their internally determined origin and scale information. Since the results from GNSS are already fixed to ITRF2014 and IGSR3 for Repro2 and Repro3, the origin and scale information represents the observed origin and scale measured in the corresponding reference frame. Additionally, to study the frame characteristics of each technique relative to ITRF, we defined the translation and scale parameters through the ITRF core stations using minimum constraints. Both orientation schemes align with either ITRF2014 or ITRF2020. As VLBI is not sensitive to the origin, when solving with an internal constraint scheme, only scale parameters are constrained, and other reference datum information is still defined by ITRF.

Table 3 illustrates the 14 transformation parameters between the Stacking Space Geodesy TRFs of single technique and the ITRF2014 (ITRF2020). The transformation parameters are computed based on the core station networks. When the origin and scale are aligned with the ITRF, the 14 transformation parameters between the Stacking Space Geodesy TRF and ITRF generally approach zero. However, when the origin and scale of the Space Geodesy TRFs are defined using internal constraints, the transformation parameters then reflect the differences between the techniques and the ITRF.

3.3. Analysis of Translation and Scale Time Series

Figure 2 presents the series of translation and scale parameters between the observations and the Stacking Space Geodesy TRFs under two types of constraints: internal constraints and minimum constraints.

When internal constraints are applied, the translation and scale parameters of GNSS characterize the observational datum precision of GNSS in its aligned reference framework, specifically ITRF2008 and IGSR3. In the GNSS 2014 Stacking TRF, under the minimum constraint definition, the translation parameter Tx and Ty are approximately 2mm, with their velocities essentially close to zero. The 2014 GNSS input data reference frame is fixed to ITRF2008, after stacking and defining the reference frame to ITRF2014, the calculated Tx and Ty values reflect the differences between ITRF2008 and ITRF2014. The drift of the translation parameter in the Z direction and the scale parameter is primarily influenced by the technical characteristics of GNSS [37,38]. In 2020, GNSS adopted a secular polar motion model during the data processing, resulting in non-zero translation velocities for both Tx and Ty where ${\dot{\mathrm{T}}}_{\mathrm{X}}=0.07\mathrm{m}\mathrm{m}/\mathrm{y}$ and ${\dot{\mathrm{T}}}_{\mathrm{y}}=-0.05\mathrm{m}\mathrm{m}/\mathrm{y}$ (by linear fit). In general, although GNSS provides numerous stations and high precision in positioning, its datum observation lacks stability, especially in the Z component and scale [39]. The GNSS coordinate calculations exhibit a significant dependence on the reference framework. Consequently, notable differences exist between the outcomes obtained under the two types of datum definitions.

From Figure 2, it is evident that the translational parameter series of SLR exhibit excellent consistency between internal and minimum constraints. This is attributed to the stability in determining the origin by the SLR. The translation velocity of ${\dot{\mathrm{T}}}_{\mathrm{y}}$ in SLR 2014 is due to the drift at around 2010, which can be also found in ITRF2014 results. From the scale parameter time series, it is noticeable that there is a drift in the SLR scale around 1998. Therefore, for defining the scale in ITRF2020, only SLR data from after 1997.75 is chosen [2]. The consistency of scale between SLR and ITRF has improved, increasing from around -5mm in 2014 to approximately -1mm in 2020.

The determination of scale parameters in VLBI exhibits a significant level of noise, with a consistency of 3mm with ITRF2014 and 2mm with ITRF2020. A noticeable drift in VLBI scale occurred around 2013, prompting ITRF2020 to only incorporate data before 2013.75 when defining the scale using VLBI. Consequently, the scale velocity derived from fitting the VLBI scale parameter series with all epochs in 2020 differs by approximately 0.21mm/year from the velocity obtained by fitting the data up to 2013.75.

The origin determinations in the X and Y directions from DORIS demonstrate accuracy comparable to SLR. However, there is notable variability in the Z direction. The translation velocity of ${\dot{\mathrm{T}}}_{\mathrm{y}}$ in DORIS 2014 may be caused by the selected core station network. The consistency in scale between DORIS and ITRF is approximately within the range of 5-10 mm, and not stable enough compared to SLR and VLBI.

3.4. Analysis of Coordinate Residuals

Figure 3 illustrates the WRMS of the station coordinate residuals for each technique across their time span. Since the WRMS of the 2014 datasets are similar to those for the 2020 datasets during the same period, only 2020 WRMS are presented. It is evident that the accuracy of DORIS, SLR, and GNSS techniques gradually improve over time. The WRMS for the DORIS in the X and Y directions decreased from 30mm to 5mm, while the Z direction WRMS decreased from 25mm to 5mm. For SLR, the WRMS in the X, Y, and Z directions decreased from 20mm to better than 10mm. In the case of GNSS, the WRMS in the X and Y directions decreased from 4mm to 1.5mm, while the Z direction WRMS decreased from 5mm to 2mm. The station coordinates’ precision of DORIS surpasses SLR slightly after 2015, with GNSS exhibiting the highest coordinate precision. The variation in DORIS station coordinate precision is closely linked to changes in the DORIS satellite constellation [22,23]. VLBI shows significant variability in position precision within individual sessions and exhibit no significant change over time.

Furthermore, Table 4 displays the average WRMS for each technique that calculated separately for each stations. The precision of stations for GNSS, SLR, and DORIS in the 2020 version is improved compared to that in 2014.

3.5. XPO and YPO Residuals Compared to 14C04 and 20C04

As depicted in Figure 4, a comparison is made between the stacking polar motion parameters (XPO, YPO) and the corresponding C04 series, namely 14C04 and 20C04. Mean residuals and standard deviations are then calculated to characterize the observational precision of polar motion parameters for each technique. The results indicate that GNSS has the highest precision in polar motion measurement, with a difference standard deviationof approximately 0.05mas compared to C04. VLBI exhibits relatively large noise in polar motion observations. The polar motion observations of SLR show differences with C04 within 0.3mas. The precision of DORIS polar motion observations gradually improves, with differences from C04 within 0.2mas after the year 2015.

4. Conclusions

In this study, we employed the TRF stacking method to assess the precision of inputs from four techniques, GNSS, SLR, VLBI, and DORIS, in ITRF2014 and ITRF2020. The evaluation encompassed outlier removal, post-rescaling of covariance matrices, datum analysis, coordinate precision analysis, and polar motion precision analysis.

Roughly 0.5% or less of the SLR observations are identified as outliers, while the ratio of DORIS, GNSS, and VLBI observations are below 1%, around 2%, and ranging from 1% to 1.2%, respectively. In 2014, the post-rescaling covariance scale factors are 25.07, 27.25, 18.84, 6.98 for GNSS, SLR, VLBI, and DORIS, respectively. The results for JTRF2014 were 32.49, 27.04, 34.81, 7.29, respectively, while ITRF2014 yielded 5.29, 25.20, 23.04, 6.86, respectively. In 2020, the post-rescaling covariance scale factors are 8.95, 14.9, 16.8, 7.78, respectively, with ITRF2020 results being 6.76, 15.44, 133.63, 8.53, respectively. It is important to note that ITRF uses daily solutions for GNSS data instead of weekly solutions, and the epoch spans for VLBI and SLR selected by us differ from those in JTRF and ITRF. Furthermore, ITRF2020 uses 12-hour epoch solutions for VLBI data, whereas our analysis utilizes mid-session epoch solutions. These are the primary reasons for the differences in covariance matrix rescale factors.

Based on the datum analysis results, it is shown than the data processing strategy for SLR has improved the consistency of SLR scale with ITRF from -5mm in ITRF2014 datasets to -1mm in ITRF2020 datasets. VLBI exhibits a drift in scale around 2013, resulting in an approximate scale difference of 0.21mm/year. Therefore, when using VLBI for scale definition, this drift needs to be taken into consideration. The stability of DORIS origin observations has been enhanced, with the scale consistency with ITRF falling within the range of 5-10mm. The datum accuracy of DORIS significantly improves with the progression of time. The rates of X and Y translation components between GNSS 2020 input data and ITRF2020 indicate ${\dot{\mathrm{T}}}_{\mathrm{X}}=0.07\mathrm{m}\mathrm{m}/\mathrm{y}$ and ${\dot{\mathrm{T}}}_{\mathrm{y}}=-0.05\mathrm{m}\mathrm{m}/\mathrm{y}$.

The decreasing standard deviations of the Polar motion parameter (XPO, YPO) discrepancies between Stacking TRFs and 14C04 (20C04) indicating an improvement in the precision of polar motion observations for all of the four techniques. The standard deviations of (XPO, YPO) offsets from 2014 to 2020 are: (0.059, 0.049) to (0.044, 0.041) for GNSS, (0.279, 0.273) to (0.238, 0.247) for SLR, (0.739, 0.731) to (0.522, 0.524) for VLBI and (0.666, 0.633) to (0.422, 0.395) for DORIS in miliarcseconds.

From the perspective of the WRMS in station coordinates, Since the inception of the technique, the station coordinate WRMS of DORIS decreased from 30mm to 5mm for X and Y components, and 25mm to 5mm for Z component; SLR WRMS decreased from 20mm to better than 10mm (X, Y and Z); GNSS WRMS decreased from 4mm to 1.5mm (X and Y) and 5mm to 2mm (Z); while VLBI shows no significant change.