1 Introduction

Terrestrial reference frames (TRFs) are represented by sets of coordinates at individual points of an Earth-fixed polyhedron. Such representations may be either static (epoch frames) or time-dependent (multi-year frames), and their principal aim is to realize a terrestrial reference system (TRS) for precise positioning on the deforming Earth. The theoretical specifications for the TRS origin, orientation and scale stem from a mixture of conventional stipulations, physical considerations and other application-oriented preferences, whereas their practical implementation relies on mathematical models, data processing strategies and datum-related constraints that are used in the frame estimation process by space geodetic techniques. Equally important to any TRF coordinate-based representation in a geodetic network is therefore its dual complement which refers to the coordinate system that is inherited by the underlying frame itself. The present study deals with TRF representations of the latter type and their accuracy assessment from the results of space geodetic solutions.

At first glance, the problem considered in this paper is rather straightforward and it can be stated in the following simple way: given the covariance (CV) matrix for a set of estimated coordinates, find the accuracy of the origin, orientation and scale of the coordinate system that is realized by these coordinates. In case of time-varying coordinates in multi-year frames, the accuracy assessment of the coordinate system extends to its initial state (at a reference epoch) and its stability relative to a set of stations with known kinematic behavior. The practical relevance of this problem is well known to geodesists, and it is dictated by the need to assess TRFs in relation to their fundamental role, namely to provide access to a well specified TRS in support of Earth observation, satellite navigation and other geodetic positioning applications.

Since any frame solution is supplied with a full CV matrix for the estimated coordinates, it is sensible to propagate their variances and co-variances to the elements of the realized coordinate system. Obviously, such propagation leads to an internal-only evaluation of the frame origin, orientation and scale, which is nevertheless essential to form a complete view of TRF quality. Its aim is to quantify, in a statistical sense, the common-mode errors (CMEs) of translational, rotational and scaling type that are present in the TRF coordinates, and thus to deduce the intrinsic accuracy of the coordinate system which is inherited by the frame solution. The CMEs are regarded here as a correlated part of the total estimation error in TRF coordinates that is concealed in their CV matrix coming out of the least-squares adjustment of space geodetic observations in global or regional networks. One of the goals in this paper is to formalize the assessment of internal CMEs in frame solutions by exploiting the error CV matrix of their estimated coordinates.

The quest to evaluate the TRF accuracy, especially for the origin and scale elements of the International Terrestrial Reference Frame (ITRF) and its successive realizations, has been a topic of primary interest to geodesists over the last decade or so. This is mainly driven by the need to maintain the ITRF as the fundamental basis for Earth system science and global change studies, and thus to meet the demands of scientific users for highly accurate and well predictable station coordinates. In the context of the Global Geodetic Observing System (GGOS), the specifications of multi-year frame accuracy for the most stringent applications (e.g., long-term monitoring of sea level rise) and their related products are known to be at the level of 1 mm and 0.1 mm/yr per each frame parameter (Plag and Pearlman 2009; Blewitt 2015) which is not fully achievable at the present time, neither in internal or external sense. Several studies have been actually performed in the past years to investigate the accuracy level of ITRF solutions using traditional approaches for assessing the origin and scale (in-)consistency between successive ITRFs and their technique-specific contributions, or external comparisons with ground/space geodetic observations and geophysical models that were not used in the ITRF construction (e.g., Altamimi et al. 2008; Collilieux and Altamimi 2013; Wu et al. 2011). Their findings demonstrated that the most recent solutions are more accurate than the previous in terms of origin and scale rate definition, yet additional analysis is still required to deliver a more precise error budget for the individual frame parameters and their long-term stability (Collilieux et al. 2014). What is more, each of these TRF evaluation methods carries its own advantages, but also limitations and drawbacks, and none of them is adequate by itself to provide a full picture of the true accuracy of the realized coordinate system. The present study contributes to the aforesaid efforts through a stochastic approach that evaluates the origin, orientation and scale accuracy of any TRF solution by an appropriate mapping of its coordinate CV matrix onto the Helmert parameter space. Such an approach complements the existing methods for the internal evaluation of TRFs, and it can be effectively applied either to secular or non-secular frame solutions from single- or multi-technique estimation schemes.

From a practical viewpoint, the inference of TRF accuracy at coordinate system level is not an effortless task to be resolved by simple error propagation on some definite analytic expression. It is essentially an inverse problem that requires a suitable projection from the coordinate parameter space to the Helmert parameter space based on the linearized similarity transformation model and other problem-specific criteria. A main challenge for geodesists is therefore the choice of a mapping rule that allows the conversion of coordinate errors in a terrestrial network to matching errors in the realized coordinate system. Several types of projectors may be considered for this problem, all of which aim at extracting the correlated part of the estimation error in TRF coordinates that is describable by small random perturbations in their coordinate system. Different projectors shall lead to different attributes in the filtered errors and likely dissimilar results in the inferred TRF accuracy at coordinate system level. The objective of this paper is to investigate not only the theoretical aspects of such projectors, but also the practical impact on the results obtained by their implementation in space geodetic solutions.

The rationale of our work is linked to other previous studies that have used the mathematical apparatus of inverse error mappings and covariance projectors to address key issues on reference frames and their quality analysis. The basic principles are rooted in Meissl’s inner error theory which was introduced in his 1962 paper as a novel framework for filtering the effect of an arbitrary set of parameters out of a given CV matrix, and on his subsequent work where the concepts of inner accuracy and free network adjustment were originally introduced into geodesy (Meissl 1965, 1969). Considering that the accuracy of estimated coordinates in geodetic networks (e.g., as formally obtained by their CV matrix after a least squares adjustment) not only reflects the accuracy of network shape characteristics but also the accuracy to which the underlying network has been attached to a particular coordinate system, the focus of Meissl’s work was to distinguish these two effects in order to evaluate the geometrical strength of geodetic networks regardless of their external datum definition. The theoretical background of this problem was further enriched by Baarda (1973) and Pope (1973, 1974) for the case of minimally constrained solutions, and also by Ebner (1974, 1975) for more general settings. Later on, the same problem was reconsidered by Blewitt et al. (1992) for removing datum-related dependencies from coordinate CV matrices of GPS global solutions in order to compare them with ITRF coordinates under a properly weighted Helmert transformation. Their work initiated the usage of covariance projectors for the filtering of frame parameters from space geodetic solutions and the efficient handling of datum constraints in combined frame estimation schemes; for more details see also Blewitt (1998). Another important work is Sillard and Boucher (2001) which focused on the sensitivity assessment of space geodetic techniques to estimable frame parameters (the so-called reference system effect according to their terminology) based on orthogonal decompositions of normal equations systems (NEQs). Their analysis offered an operational framework to identify the degrees-of-freedom of the TRS that are not practically reduced by space geodetic observations in terrestrial networks, and thus to detect the truly required minimal constraints for datum definition purposes. The present study is somewhat similar to the previous paper, yet it refers to a different inferential problem, namely the evaluation of the internal accuracy of the realized coordinate system in TRFs that are estimated by space geodetic techniques (and typically by an external set of datum constraints). In fact, the proposed covariance projector for solving this problem is different from the one suggested in Sillard and Boucher (2001), as it will be explained in more detail later in the paper.

The rest of the paper is organized as follows: Sect. 2 presents the theoretical background and explains the mapping problem of frame accuracy at coordinate system level; the influence of datum constraints and their weighting level on the CMEs of frame solutions is discussed in Sect. 3, whereas additional key aspects for scientific users of TRFs, namely the long-term degradation of multi-year frames and their forward-in-time accuracy assessment are given in Sect. 4; in Sect. 5, the theoretical framework of the previous sections is utilized to perform the internal accuracy evaluation of the realized origin, orientation and scale in the ITRF frame series based on the formal CV matrices for their estimated positions and velocities in the four technique subnetworks (DORIS, SLR, VLBI, GNSS); finally, in Sect. 6 some conclusions and prospects for future work are given.

2 Theoretical background

2.1 Coordinate frames and their common-mode errors

The notion of CMEs in geodetic frames is related to mapping small disturbances of a 3D Cartesian coordinate system to small disturbances in the coordinates of a terrestrial network with respect to that system (and vice versa). Hereafter, these disturbances shall be treated as zero-mean random variables with associated covariance structure which will be elucidated later in the paper. Their joint relation is expressible by the linear formula

$${{{\varvec{d}}{\varvec{X}}_{{\varvec{\theta}}} }} = {{{{\varvec{A}} \, {\varvec{d}}{\varvec{\theta}}}}}$$
(1)

which represents the classic 7-parameter or 14-parameter Helmert transformation (HT) in compact differential form. In case of static or epoch TRFs, the vector dθ contains the origin, orientation and scale offsets of the coordinate system (3 translations, 3 rotations, 1 scale factor), whereas in secular TRFs it is augmented by additional elements to account also for drifts in the time-variable part of the coordinate system (3 translation rates, 3 rotation rates, 1 scale rate). Accordingly, the vector \({\varvec{d}}{\varvec{X}}_{{\varvec{\theta}}}\) gives the variations of frame coordinates (either disturbed positions or disturbed positions and velocities) and A is the HT’s Jacobian matrix in compliance with the frame type to be analyzed in the considered network.

Every TRF solution from space geodetic data is affected by random estimation errors, a part of which is always describable through the linear disturbance model of Eq. (1). This parameterized error part constitutes the CMEs that exist in any set of frame coordinates obtained from inexact measurements via a least squares estimation process. Specifically it contains the sum of net-translational, net-rotational and net-dilatational random errors, and it reflects the (internal) accuracy of the coordinate system that is realized by the TRF solution itself. The remaining error part is independent of the coordinate system and it dictates the geometrical accuracy of the frame polyhedron that is defined by the station coordinates. In principle, the following error decomposition is therefore relevant to any TRF solution

$$ {\varvec{d}}{\varvec{X}} = \underbrace{{{\varvec{d}}{\varvec{X}}_{{\varvec{\theta}}} }}_ {{{{\varvec{A}} \, {\varvec{d}}{\varvec{\theta}}}}}\,\,+\,\,{\varvec{d}}\overline{\varvec{X}} $$
(2)

where dX denotes the random errors relating to the solution’s CV matrix \({\varvec{\varSigma}}_{{\varvec{X}}}\), while the meaning of the other terms is self-evident from the previous discussion. A simplified illustration of this error decomposition is shown in Fig. 1. Although the above equation is useless for direct computations, it provides a convenient framework to infer the statistical behavior of CMEs from the formal CV matrix of frame coordinate errors. The treatment of this “error mapping” problem emulates the least-squares inversion of Eq. (2), and it entails various open issues to be discussed next in this paper.

Fig. 1
figure 1

Graphical illustration of the error decomposition in geodetic TRFs according to Eq. (2). The additive components consist of a parametric part (CMEs) which signifies the internal accuracy of the realized coordinate system, and a nonparametric part (inner errors) that reflects the geometrical accuracy of the frame polyhedron as defined by the station coordinates

Full size image

2.2 Relevance of Eq. (2) in frame accuracy analysis

The error splitting in Eq. (2) is a modeling choice for analyzing the internal accuracy of TRFs at coordinate system level. Its rationale is based on the fact that frame coordinates define at the same time two distinct features, that is (1) the shape of an Earth-fixed polyhedron and (2) its relative placement with respect to a coordinate system, and thus, their estimation errors will contaminate both of these features. In fact, the only assumption that we make here is to delineate the resulting errors in those features with a simple additive model, without necessarily implying a null correlation between them.

In line with the above reasoning, Eq. (2) is applicable in all different kinds of TRF solutions, including minimally constrained (MC) solutions, over-constrained solutions, or even free solutions. However, the proper interpretation of the error term dθ should be carefully considered in each of the above cases. For example, in unconstrained solutions the CMEs reflect the frame accuracy at coordinate system level that is attainable only by the used observations in a given network. This is useful for assessing the sensitivity of space geodetic techniques to estimable frame parameters (origin, scale), as it was already proposed in the work of Sillard and Boucher (2001). Conversely, the CMEs in partially or fully constrained solutions are affected by the datum constraints and their interaction with the available observations during the frame estimation process. In this case, the term dθ reflects the random errors in the realized coordinate system (in the form of unknown offsets and drifts for the origin, orientation and scale), while its CV matrix is a useful metric to assess the frame’s compliance to formal accuracy standards for various stringent applications (e.g., determination of long-term sea level change, monitoring of vertical land motion, multi-technique satellite orbit determination, etc.).

It is essential to clarify that the merit of Eq. (2) is not to perform its least-squares inversion in an optimal statistical sense. Specifically, the error mapping problem that is considered here is not rectified by setting up an optimal estimator via the Gauss–Markov formulation, and it would be misleading to adopt such a view for the inference of CMEs in frame solutions. In fact, the sought solution for dθ should not be linked with the best (minimum variance) linear unbiased estimator, as this dictates that the error components \({\varvec{d}}{\varvec{X}}_{{\varvec{\theta}}}\) and \({\varvec{d}}\overline{\varvec{X}}\) become de facto uncorrelated. This is a false presumption since both components are supposed to be random quantities with joint (yet unknown) dependence on the frame coordinate errors. Their de-correlation is a mathematical artifice that can be used to resolve our error mapping problem, but it will also lead to an unorthodox view of the frame’s inner accuracy and, most likely, to under-estimated CMEs in the frame coordinates; more details and numerical examples are given in later sections.

2.3 Inference of CMEs in frame solutions

In general terms, the error mapping of Eq. (2) is resolved in a linear context as follows

$$ {\varvec{d}}{\varvec{X}} \, = \, \underbrace{{{\varvec{Q}}{\varvec{d}}{\varvec{X}}}}_{{{\varvec{d}}{\varvec{X}}_{{\varvec{\theta}}} }}\,\,+\,\,\underbrace{{({\varvec{I}} - {\varvec{Q}}) \, {\varvec{d}}{\varvec{X}}}}_{{{\varvec{d}}\overline{\varvec{X}}}} \, $$
(3)

where Q is a projection matrix (Q2 = Q) to be chosen on the basis of desired attributes for the individual error components. Since the aim here is to extract the hidden CMEs in TRFs, this matrix has to preserve the correlated part of frame coordinate errors which is describable by the disturbance model of Eq. (1) (and it has to filter out the remaining part). Consequently, a necessary property of the projection matrix is

$$ {\varvec{Q}} \, {\varvec{A}} \, = \, {\varvec{A}} $$
(4)

which ensures that the CMEs in dX are absorbed by the error component \({\varvec{d}}{\varvec{X}}_{{\varvec{\theta}}}\) (recall that A is the HT’s Jacobian matrix conforming to the type of frame, epoch or secular, to be analyzed). Another critical aspect of the projection matrix is that it restrains the inner part of frame coordinate errors (i.e., the complement of CMEs) through the formula

$$ {\varvec{Q}} \, {\varvec{d}}\overline{\varvec{X}} \, = \, {\mathbf{0}} $$
(5)

which is deduced from Eq. (3). This, in turn, confines the projector choices to those that can assign a meaningful geometrical interpretation to Eq. (5) and they ensure the invariance of \({\varvec{d}}\overline{\varvec{X}}\) from the underlying coordinate system.

A suitable class of projection matrices for the previous error decomposition is given by the general expression

$$ {\varvec{Q}} \, = \, {\varvec{A}} \, ({\varvec{A}}^{T} {\varvec{PA}})^{ - 1} {\varvec{A}}^{T} {\varvec{P}} $$
(6)

where P is a weight matrix for the frame coordinates. The use of this projector is equivalent to applying the weighted least-squares inversion of Eq. (2), and it implies the following inference for the CV matrix of the error term \({\varvec{d}}{\varvec{\theta}}\):

$${\varvec{\varSigma}}_{{\varvec{\theta}}} \, = \, ({\varvec{A}}^{T} {\varvec{PA}})^{ - 1} {\varvec{A}}^{T} {\varvec{P}} \,{\varvec{\varSigma}}_{{\varvec{X}}} {\varvec{PA}} \, ({\varvec{A}}^{T} {\varvec{PA}})^{ - 1} $$
(7)

The above matrix gives the internal accuracy of the realized origin, orientation and scale for any frame solution by exploiting the formal CV matrix of the estimated coordinates. The notion of inner errors in this case complies with the condition

$$ {\varvec{A}}^{T} {\varvec{P}} \, {\varvec{d}}\overline{\varvec{X}} \, = \, {\mathbf{0}} $$
(8)

which is compatible with Eq. (5), yet its geometrical interpretation is obscured by the unstipulated weight matrix P.

2.4 Choice of the weight matrix

To utilize the previous error mapping in practice, we should settle the choice of the weight matrix for the frame coordinates. Two possible options are considered here, namely

$$ {\varvec{P}} ={\varvec{\varSigma}}_{{\varvec{X}}}^{ - 1} \;\;{\text{or}}\;\;{\varvec{P}} = \varvec{\rm I} $$
(9)

The first option is the one suggested in Sillard and Boucher (2001) on the basis of its de-correlation property for the projected error components, that is

$$ {\varvec{d}}\overline{\varvec{X}} = \, ({\varvec{I}} - {\varvec{A}} \, ({\varvec{A}}^{T}{\varvec{\varSigma}}_{{\varvec{X}}}^{ - 1} {\varvec{A}})^{ - 1} {\varvec{A}}^{T}{\varvec{\varSigma}}_{{\varvec{X}}}^{ - 1} ) \, {\varvec{d}}{\varvec{X}} $$
(10a)
$$ {\varvec{d}}{\varvec{X}}_{{\varvec{\theta}}} = \, {\varvec{A}} \, ({\varvec{A}}^{T}{\varvec{\varSigma}}_{{\varvec{X}}}^{ - 1} {\varvec{A}})^{ - 1} {\varvec{A}}^{T}{\varvec{\varSigma}}_{{\varvec{X}}}^{ - 1} \, {\varvec{d}}{\varvec{X}} \, \doteq \, {\varvec{A}} \, {\varvec{d}}{\varvec{\theta}} $$
(10b)
$${\varvec{\varSigma}}_{{{\varvec{d}}\overline{\varvec{X}},{\varvec{d}}{\varvec{X}}_{{\varvec{\theta}}} }} = \, {\mathbf{0}} $$
(10c)

In this case, the frame accuracy at coordinate system level is evaluated by the CV matrix

$${\varvec{\varSigma}}_{{\varvec{\theta}}} \, = \, ({\varvec{A}}^{T}{\varvec{\varSigma}}_{{\varvec{X}}}^{ - 1} {\varvec{A}})^{ - 1} $$
(11)

which has minimum trace among all alternate inferences given by Eq. (7) for different weight matrix choices (the above matrix is identical with the result given in the third line of Eq. (18) of Sillard and Boucher (2001)). This suggests a tight evaluation of CMEs by assigning to them the smallest possible contribution to the TRF estimation error. Hence, Eq. (11) sets a useful theoretical bound in the sense of defining the “best case scenario” for the internal accuracy of the realized coordinate system. The disadvantage of this weighting option is that the resulting inner errors cannot admit a clear-cut geometrical interpretation and their invariance with respect to the frame’s coordinate system becomes dubious. In particular, the inner errors defined by Eq. (10a) satisfy the condition

$$ {\varvec{A}}^{T}{\varvec{\varSigma}}_{{\varvec{X}}}^{ - 1} {\varvec{d}}\overline{\varvec{X}} \, = \, {\mathbf{0}} $$
(12)

which does not ensure the absence of net-translational, net-rotational and net-scaling components in their values! As a result, the CV matrix in Eq. (11) is likely to give an optimistic view of frame accuracy that may not reflect the actual CMEs in the estimated coordinates.

The second option employs a unit weight matrix for the frame coordinates, and it is equivalent to applying the ordinary least-squares inversion of Eq. (2), that is

$$ {\varvec{d}}\overline{\varvec{X}} = \, ({\varvec{I}} - {\varvec{A}} \, ({\varvec{A}}^{T} {\varvec{A}})^{ - 1} {\varvec{A}}^{T} ) \, {\varvec{d}}{\varvec{X}} $$
(13a)
$$ {\varvec{d}}{\varvec{X}}_{{\varvec{\theta}}} = \, {\varvec{A}} \, ({\varvec{A}}^{T} {\varvec{A}})^{ - 1} {\varvec{A}}^{T} {\varvec{d}}{\varvec{X}} \, \doteq \, {\varvec{A}} \, {\varvec{d}}{\varvec{\theta}} $$
(13b)
$${\varvec{\varSigma}}_{{{\varvec{d}}\overline{\varvec{X}},{\varvec{d}}{\varvec{X}}_{{\varvec{\theta}}} }} = \, ({\varvec{I}} - {\varvec{A}} \, ({\varvec{A}}^{T} {\varvec{A}})^{ - 1} {\varvec{A}}^{T} ) \,{\varvec{\varSigma}}_{{\varvec{X}}} {\varvec{A}} \, ({\varvec{A}}^{T} {\varvec{A}})^{ - 1} {\varvec{A}}^{T} \ne \, {\mathbf{0}} $$
(13c)

It is essential to clarify that the choice \({\varvec{P}} = \varvec{\rm I}\) does not mean that the frame coordinates are assumed to have the same quality or being uncorrelated with each other, but it only dictates that the frame stations are equally contributing to the assessment of the realized coordinate system. The frame accuracy at coordinate system level is now given by the projected CV matrix

$${\varvec{\varSigma}}_{{\varvec{\theta}}} = \, ({\varvec{A}}^{T} {\varvec{A}})^{ - 1} {\varvec{A}}^{T}{\varvec{\varSigma}}_{{\varvec{X}}} {\varvec{A}} \, ({\varvec{A}}^{T} {\varvec{A}})^{ - 1} $$
(14)

which generally leads to different results compared to Eq. (11), whereas the respective inner errors fulfill the condition

$$ {\varvec{A}}^{T} \varvec{d\overline{X}} = {\mathbf{0}} $$
(15)

The above condition ensures the nullification of net-translational, net-rotational and net-scaling components in \({\varvec{d}}\overline{\varvec{X}}\), thus complying with the traditional perception of inner errors in geodetic networks. The independence between the CMEs and inner errors is assured by the orthogonality condition

$$ {\varvec{d}}\overline{\varvec{X}}^{T} {\varvec{d}}{\varvec{X}}_{{\varvec{\theta}}} = 0 $$
(16)

which is easily verified from Eqs. (13a) and (13b). A brief summary of the two weighting options of Eq. (9) and their associated properties is given in Table 1.

Table 1 Summary of projection schemes for the internal assessment of CMEs and the evaluation of frame accuracy at coordinate system level
Full size table

The projected error components \({\varvec{d}}\overline{\varvec{X}}\) and \({\varvec{d}}{\varvec{X}}_{{\varvec{\theta}}}\) by the second weighting option remain stochastically correlated with each other, yet this is not problematic in any way and it offers a more realistic treatment of the error mapping problem at hand. The stochastic correlations between CMEs and inner errors, as expressed by Eq. (13c), are a manifestation of the fact that the shape of the frame polyhedron and its relative placement in a coordinate system do not occur from independent procedures, but they always stem from the same estimation process and a single set of noisy coordinates (this fact is ignored by the first weighting option).

An additional aspect of the second option is that the uniform weighting of frame stations honors to full extent their geometrical configuration and its contribution to the realized coordinate system. On the other hand, the uniform weighting may be seen as a drawback in the sense that it equates the contribution of stations with different time spans of observation history and dissimilar quality in their estimated coordinates. This argument can be counteracted by the fact that an operational TRF is not just a selected core-network of a limited number of high-quality stations, but it consists of all stations that have participated in its estimation phase from space geodetic observations. Obviously, to avoid an “unfair” evaluation of frame accuracy due to the influence of few bad stations with low-accuracy coordinates, it is sensible to remove them in advance from the CV matrix \({\varvec{\varSigma}}_{{\varvec{X}}}\) before applying the covariance projector of Eq. (14) or Eq. (11).

An empirical “middle-ground” approach that reconciles the two weighting options of Eq. (9) could be also considered to support the TRF accuracy evaluation. It would seem reasonable, for example, to have diagonal weights in \({\varvec{P}}\) that reflect the positive impact a station might have on the reference frame (which can be considered less extreme than inclusion or rejection of a station in the frame evaluation process). Such an alternative approach of diagonal weighting has been addressed in Blewitt et al. (1992), yet it will not be pursued in the present study.

Based on the previous discussions, the weight matrix \({\varvec{P}} = \varvec{\rm I}\) is perceived as the most suitable choice for evaluating the (internal) accuracy of frame solutions at coordinate system level. In fact, any other weighting choice does not properly isolate the CMEs since the matching inner errors will not be free of net-translation, net-rotation and net-scaling components as per Eq. (15). The same preference was also endorsed by Rebischung (2014) with regard to a similar inferential problem, namely the estimation of implicit (datum related) parameters in linear adjustment models of space geodetic observations. Following a different approach from the one presented here, his study showed that the alternate weight matrix \({\varvec{P}} ={\varvec{\varSigma}}_{{\varvec{X}}}^{ - 1}\) discards the correlations between datum parameters, network geometry and other model-based factors, and it results in a tighter assessment of implicit parameters (and their accuracy level) in space geodetic solutions.

As a last remark, let us note that the frame accuracy assessment at coordinate system level relies on the correctness of the CV matrix \({\varvec{\varSigma}}_{{\varvec{X}}}\) which is obtained from the least squares adjustment of space geodetic observations in global or regional networks. It is therefore important that any scaling issues or other mis-modeling aspects of this matrix are clearly known before its proper exploitation through the covariance projection schemes that were described in this section.

2.5 Numerical example

An instructive example from the implementation of the previous error covariance projectors is given here. In particular, Eqs. (11) and (14) were used to infer the frame accuracy at coordinate system level of the weekly combined solutions in the EUREF Permanent GNSS Network (EPN) for the period 1999–2021 (GPS weeks 1000–2154). The plots in Fig. 2 depict the square roots of the diagonal elements of the weekly CV matrices {\({\varvec{\varSigma}}_{{\varvec{\theta}}}\)} which were computed by each weighting option based on the official SINEX files of the EPN weekly solutions. To allow easier visualization of the results, the accuracy of the orientation and scale parameters has been expressed in mm by multiplying their standard deviations with the mean Earth radius. Note that the EPN weekly solutions are tied to various ITRF realizations—mostly in accordance with the International GNSS Service (IGS) reference frame updates—and they have not all used the same alignment strategy nor the same fiducial set of IGS reference stations. For more details on the EPN post-processed products, see the information given at https://www.epncb.oma.be/.

Fig. 2
figure 2

Internal accuracy of the combined EPN weekly solutions at coordinate system level. The two plots refer to different weighting options for the weekly frame coordinates, namely (a) \({\varvec{P}} ={\varvec{\varSigma}}_{{\varvec{X}}}^{ - 1}\) and (b) \({\varvec{P}} = \varvec{\rm I}\). The alignment method of the EPN weekly solutions to the ITRF changed from heavy constraints (HCs) to minimal constraints (MCs) in the GPS week 1303 (December 2004)

Full size image

By comparing the results in Fig. 2, it is confirmed that the first option (\({\varvec{P}} ={\varvec{\varSigma}}_{{\varvec{X}}}^{ - 1}\), left plot) overestimates the frame accuracy compared to the second option (\({\varvec{P}} = \varvec{\rm I}\), right plot) throughout the entire time period. This is more emphatically seen in the results for the pre-2005 period, where the alignment of the EPN weekly solutions to ITRF is based on heavy constraints at IGS fiducial stations. For that period, the two weighting options show differences that reach up to two orders of magnitude for the internal accuracy of the frame parameters. These large differences should be attributed to the fact that CMEs and inner errors are strongly correlated with each other in over-constrained networks, therefore their de-correlation under the first weighting option (see Eq. 10c) yields a “false” CV matrix \({\varvec{\varSigma}}_{{\varvec{\theta}}}\) with minimum trace that does not reflect in a realistic way the frame accuracy at coordinate system level. After the switch of the EPN alignment strategy to minimal constraints, the two weighting options show much better consistency and their differences are gradually reduced from few mm down to sub-mm level.

A detailed analysis of the results obtained by the second weighting option in the EPN network, along with an explanation of the frame accuracy “offset” that is visible at the beginning of 2010 in both plots of Fig. 2, is given in a recent paper by Kotsakis and Chatzinikos (2022).

3 The role of datum constraints and their weighting

Let us consider the analytic form of the CV matrix of estimated frame coordinates that is usually stored in SINEX files, that is

$${\varvec{\varSigma}}_{{\varvec{X}}} = \, \hat{\sigma }^{2} \, ({\varvec{N}} + {\varvec{H}}^{T} {\varvec{WH}})^{ - 1} $$
(17)

where N is the normal matrix of the used observations (after reducing other auxiliary parameters from the original NEQ), H is the design matrix of datum constraints, W is the weight matrix of datum constraints, and \(\hat{\sigma }^{2}\) is the a posteriori variance factor of the estimated solution. The above expression applies to several (but not all) cases of frame estimation problems that appear in geodetic practice, and it relies on the following assumptions: (a) the datum constraints carry their own prior accuracy level that is specified by the matrix \({\varvec{W}}^{ - 1}\) (typically this is a diagonal matrix set by the user to control the impact of datum constraints into the final solution), and (b) the aforesaid accuracy level is propagated to the variances and co-variances of the estimated frame coordinates.

After considering Eq. (17), a crucial distinction needs to be made between the projected CV matrix \({\varvec{\varSigma}}_{{\varvec{\theta}}}\) and the inverse weight matrix \({\varvec{W}}^{ - 1}\) which is embedded in TRF solutions. The former matrix is computable by Eq. (14), and it gives the internal accuracy of the coordinate system that is realized through all frame stations and their estimated coordinates. The latter matrix refers to the prior accuracy of datum constraints which are typically enforced in a subset of frame stations. The former matrix is influenced by the latter without necessarily reproducing its a priori weights for the datum definition. To put it in simple words: the accuracy level at which the frame origin, orientation and scale are defined through a set of datum constraints in a subset of reference stations does not coincide with the accuracy at which those elements are being realized by the final solution over the entire frame network. Several other factors that are hidden in the CV matrix of Eq. (17), such as the network geometry, the type of constraints and the spatial distribution of fiducial stations, the noise level and datum content (if any) of the available observations, as well as their joint interaction during the frame estimation process, play a key role and they also affect the accuracy of the realized coordinate system.

There is a special case, however, where the impact of the aforesaid factors is rescinded and the internal accuracy of the realized coordinate system shall be controlled solely by the formal weighting of the datum constraints. This situation occurs if the so-called full inner constraints (a particular type of MCs with optimal properties for the entire network, see, e.g., Blaha 1971, 1982) are applied in frame estimation using the design matrix

$$ {\varvec{H}} \, = \, ({\varvec{A}}^{T} {\varvec{A}})^{ - 1} {\varvec{A}}^{T} $$
(18)

Any TRF solution by such datum constraints is essentially free from (internal) CMEs in the sense that the estimation errors of frame coordinates will not contain net-translation, net-rotation and net-scaling components. In this case, the observation noise and network geometry do not affect the accuracy of the realized coordinate system, a fact that is verified by Eq. (14) as follows

$$ \begin{aligned}{\varvec{\varSigma}}_{{\varvec{\theta}}} \, & = \, ({\varvec{A}}^{T} {\varvec{A}})^{ - 1} {\varvec{A}}^{T} \, \underbrace {{\hat{\sigma }^{2} ({\varvec{N}} + {\varvec{H}}^{T} {\varvec{WH}})^{ - 1} }}_{{{\varvec{\varSigma}}_{{\varvec{X}}} }} \, {\varvec{A}} \, ({\varvec{A}}^{T} {\varvec{A}})^{ - 1} \\ & = \, \hat{\sigma }^{2} ({\varvec{A}}^{T} {\varvec{A}})^{ - 1} {\varvec{A}}^{T} \, ({\varvec{N}} + {\varvec{A}}({\varvec{A}}^{T} {\varvec{A}})^{ - 1} {\varvec{W}}({\varvec{A}}^{T} {\varvec{A}})^{ - 1} {\varvec{A}}^{T} )^{ - 1} {\varvec{A}} \, ({\varvec{A}}^{T} {\varvec{A}})^{ - 1} \\ & = \, \hat{\sigma }^{2} \, {\varvec{W}}^{ - 1} \\ \end{aligned} $$
(19)

The above result stems from the matrix identity \({\varvec{H}}({\varvec{N}} + {\varvec{H}}^{T} {\varvec{WH}})^{ - 1} {\varvec{H}}^{T} = {\varvec{W}}^{ - 1}\) which holds for any minimally constrained normal matrix (see Kotsakis 2012), yet it has limited interest due to the practical irrelevance of full inner constraints for frame estimation problems. The datum definition is hardly applied on all frame stations and it usually involves a smaller subset of stations, or even separate subsets for each constrained element of the coordinate system. On the other hand, the covariance projector of Eq. (14) is applied over all frame stations and it will not reproduce the prior weights of datum constraints at the selected reference stations.

The special result of Eq. (19) is also verified if we use the alternative covariance projector of Eq. (11). The proof is straightforward after taking into account the fundamental property of the rank-deficient normal matrix, that is \({\varvec{NA}} = {\mathbf{0}}\) (Kotsakis 2012).

4 Forward prediction of internal accuracy in multi-year frames

A key aspect of a successful TRF is the ability to predict accurate station coordinates on the deformable Earth at any time needed (Blewitt 2015). For this reason, it is essential to know the accuracy level at which the coordinate system of multi-year frames is retained in their primary networks. It is actually expected that the internal accuracy of multi-year frames deteriorates in time, a fact that can degrade their long-term “predictability” for the dynamic state of the Earth system (e.g., sea level change, vertical land motion, ice melting rates, etc.).

To assess how the internal accuracy of multi-year frames varies in time, we need to combine the covariance mapping of Eq. (14) and the regularized kinematic model of frame coordinates. This shall be demonstrated here for the case of secular frames, but the same procedure can be also extended in non-secular frames with modeled seasonal displacements (e.g., ITRF2020). For the secular case, the regularized coordinates of the frame stations at any time instant are given by the linear model

$$ {\varvec{X}}(t) \, = \, {\varvec{X}}(t_{o} ) \, + \, (t - t_{o} ) \cdot \dot{\varvec{X}} $$
(20)

while their formal CV matrix is analytically expressed as

$$\begin{aligned} {\varvec{\varSigma}}_{{{\varvec{X}}(t)}} \, & = \,{\varvec{\varSigma}}_{{{\varvec{X}}(t_{o} )}} \, + \, (t - t_{o} )^{2}{\varvec{\varSigma}}_{{\dot{\varvec{X}}}} \, + \, (t - t_{o} ) \\ &\quad\times \left( {{\varvec{\varSigma}}_{{{\varvec{X}}(t_{o} ),\dot{\varvec{X}}}} \, +{\varvec{\varSigma}}_{{{\varvec{X}}(t_{o} ),\dot{\varvec{X}}}}^{T} } \right) \end{aligned} $$
(21)

where the individual terms at the right side stem from the full CV matrix of the multi-year frame solution. The latter can be projected at the coordinate system level by using the Jacobian matrix of the 14-parameter HT model according to the general scheme

$$\begin{aligned} & \varvec{\varSigma }_{\varvec{X}} = \left[ {\begin{array}{*{20}c} {\varvec{\varSigma }_{{\varvec{X}(t_{o} )}} } & {\varvec{\varSigma }_{{\varvec{X}(t_{o} ),\varvec{\dot{X}}}} } \\ {\varvec{\varSigma }_{{\varvec{\dot{X}},\varvec{X}(t_{o} )}} } & {\varvec{\varSigma }_{{\varvec{\dot{X}}}} } \\ \end{array} } \right]{\text{ }}\xrightarrow{{{\rm{Eq}}.\;(14)}}{\text{ }}\\ & \varvec{\varSigma }_{\varvec{\theta }} = \left[ {\begin{array}{*{20}c} {\varvec{\varSigma }_{{\varvec{\theta }(t_{o} )}} } & {\varvec{\varSigma }_{{\varvec{\theta }(t_{o} ),\varvec{\dot{\theta }}}} } \\ {\varvec{\varSigma }_{{\varvec{\dot{\theta }},\varvec{\theta }(t_{o} )}} } & {\varvec{\varSigma }_{{\varvec{\dot{\theta }}}} } \\ \end{array} } \right]\end{aligned} $$
(22)

The above result offers a detached representation of TRF accuracy that refers to the initial state (at a reference epoch to) and the stability of the realized coordinate system. Alternatively, we can apply the covariance projector of Eq. (14) directly to Eq. (21) by using the Jacobian matrix of the 7-parameter HT model, that is

$${\varvec{\varSigma}}_{{{\varvec{\theta}}(t)}} \, = \, ({\varvec{A}}^{T} {\varvec{A}})^{ - 1} {\varvec{A}}^{T}{\varvec{\varSigma}}_{{{\varvec{X}}(t)}} \, {\varvec{A}} \, ({\varvec{A}}^{T} {\varvec{A}})^{ - 1} $$
(23)

which is also expressible in the equivalent form (after simple algebraic manipulation)

$${\varvec{\varSigma}}_{{{\varvec{\theta}}(t_{o} + \delta t)}} \, = \,{\varvec{\varSigma}}_{{{\varvec{\theta}}(t_{o} )}} \, + \, \delta t^{2}{\varvec{\varSigma}}_{{\dot{\varvec{\theta }}}} \, + \, \delta t\left( {{\varvec{\varSigma}}_{{{\varvec{\theta}}(t_{o} ),\dot{\varvec{\theta }}}} \, +{\varvec{\varSigma}}_{{{\varvec{\theta}}(t_{o} ),\dot{\varvec{\theta }}}}^{T} } \right) $$
(24)

where the individual matrices at the right side of the last equation are exactly the same with the ones obtained under the composite mapping of Eq. (22).

As already stated, the TRF specifications for the most demanding users dictate that the origin, orientation and scale should be determined at an accuracy level of 1 mm, and they should remain stable over time at a rate of 0.1 mm/yr (Plag and Pearlman 2009). These specifications are aligned with the partitioned representation of TRF accuracy in Eq. (22), but they seem to overlook the added uncertainty in the realized coordinate system due to the cross-correlations between \({\varvec{\theta}}(t_{o} )\) and \(\dot{\varvec{\theta }}\) as shown in Eq. (24). Secular frames do not theoretically have a linear-in-time evolution of their internal accuracy and may even exhibit an interim improvement over time in case of strong negative correlations between \({\varvec{\theta}}(t_{o} )\) and \(\dot{\varvec{\theta }}\) (this is largely depended upon the cross-correlations between \({\varvec{X}}(t_{o} )\) and \(\dot{\varvec{X}}\)).

Considering the above remarks, the time-dependent CV matrix in Eq. (24) should be preferred for checking the compliance of the internal accuracy of multi-year frames with imposed quality specifications on their realized coordinate system. Its practical use shall be demonstrated with various examples in the following section.

5 Internal evaluation of the ITRF series

5.1 Rationale, data and methodology

The goal in this section is to evaluate the accuracy of the realized origin, orientation and scale in the ITRF series based on the formal CV matrices for their estimated positions and velocities in the four technique subnetworks (DORIS, SLR, VLBI, GNSS). These matrices are the ones obtained at the final step of the frame estimation process and they were extracted from the respective SINEX files which are available in the ITRF website (https://itrf.ign.fr/). The frame releases that are considered here include ITRF2000, ITRF2005, ITRF2008 and ITRF2014. The latest release ITRF2020 is quite distinct from the previous ones due to the non-secular modeling of station motions in the inter-technique combination phase (Altamimi et al. 2023), and it will be investigated in another paper. More details about the tested solutions and their estimation strategies can be found in the papers by Altamimi et al. (2002, 2007, 2011, 2016).

For our investigation, we exploit the covariance projector of Eq. (14) in conjunction with the Jacobian matrix of the 14-parameter HT model. The reasons for using this projector instead of the alternative form of Eq. (11) have been thoroughly discussed in Sect. 2. The computations were performed separately for each subnetwork, thus leading to multiple “views” of the ITRF accuracy at coordinate system level (i.e., a separate matrix \({\varvec{\varSigma}}_{{\varvec{\theta}}}\) is obtained for every technique subnetwork in each ITRF release). To prevent any misreadings of our findings, it is clarified in advance that what we evaluate here is the internal accuracy at which the ITRF origin, orientation and scale are realized in the technique subnetworks by their estimated coordinates (station positions and velocities) from successive ITRF solutions. The presented results convey the full statistical information of the CV matrices of these solutions, along with the impact of network geometry on the realized frame parameters. The inferred accuracy, for example, of the ITRF origin in the GNSS subnetwork does not reflect its likely offset (and drift) relative to the true mean CM, but rather it indicates the uncertainty level at which the ITRF origin (that is, the mean CM as sensed by SLR) is reproduced by the estimated coordinates of the GNSS stations. Forward in-time projections have been also computed by Eq. (24) to assess the degradation of the frame parameters in ITRF solutions, with emphasis on their decadal performance.

The present study is performed on the entire set of ITRF stations that appear in the SINEX files of the official solutions. No attempt is made to form “core” subnetworks based on selective criteria for the used stations, e.g., to exclude stations by considering their observation time span or other conditions related to geological site stability. Hence, the results shown here refer to the total ITRF network (and its technique-specific subnetworks) of the successive solutions, and they may be slightly pessimistic compared to the achievable frame accuracy in other custom-built clusters of ITRF stations. Even so, a number of low-accuracy stations were removed beforehand from our analysis to avoid undue distortions in the computed results (this is mostly crucial for the VLBI and SLR subnetworks). Their removal was performed by deleting their respective rows/columns from the full CV matrices in which those stations were detected. In particular, we removed from each subnetwork all points whose estimated positions and velocities had sigma values larger than 5 cm and 5 cm/yr, respectively. The total number of stations that were finally used for evaluating the internal accuracy in the ITRF series is given in Table 2.

Table 2 Number of stations in the technique subnetworks that are used for the internal accuracy assessment of the ITRF solutions
Full size table

The key aspects to look for in the results of this investigation are:

  • a gradual improvement of the internal accuracy from one ITRF solution to the next; and

  • a tendency to approach (or even exceed) the accuracy level of 1 mm and 0.1 mm/yr for the realized frame parameters in the technique-specific ITRF subnetworks. Of course, this is not sufficient to ensure the compliance with the GGOS-related specifications as it refers only to the internal accuracy at which the ITRF origin, orientation and scale are accessible through the individual subnetworks (and their estimated coordinates), yet it is essential to obtain an insightful view of the ITRF quality. Additional evaluations with external geodetic observations and/or geophysical models are obviously needed in order to assess the true accuracy level of the ITRF in conformity with its theoretical setting (see also the related discussion in Sect. 1).

Our findings seem to verify, to a large extent, both of the aforesaid aspects, and they are discussed in the next sections.

5.2 Frame realization accuracy in ITRF solutions

A series of CV matrices {\({\varvec{\varSigma}}_{{\varvec{\theta}}}\)} with dimensions 14 × 14 have been computed in accordance with the methodology that was outlined in the previous section. Herein, we focus on their diagonal elements which expose the accuracy of the initial state and rate-of-change of the realized frame parameters in different ITRF solutions (and their technique subnetworks). Note that the reference epoch of the initial state varies for each tested solution, namely 1997.0 (ITRF2000), 2000.0 (ITRF2005), 2005.0 (ITRF2008) and 2010.0 (ITRF2014).

As far as the frame origin, the growth of its internal accuracy is depicted in a number of plots given in Fig. 3. The improvement, from one ITRF solution to the next, is mostly visible in the stability of the origin components (translation rates along the three Cartesian axes) while similar behavior is also observed in the respective plots for the scale and orientation parameters shown in Figs. 4 and 5. A reasonable explanation is that the increasing time span of space geodetic observations (which are used in each ITRF solution) has greater influence on the stability, and less on the initial state, of the realized frame parameters. The ITRF2014 origin rates in the DORIS, SLR and GNSS subnetworks exhibit an internal accuracy around 0.1 mm/yr (or better), whereas in the VLBI subnetwork only the TZ rate conforms to this accuracy level. In previous solutions, especially ITRF2000 and ITRF2005, the internal accuracy of the origin rates is lower by a factor of 3 or even more, depending on the technique subnetwork (see Fig. 3). It is noted that the origin stability in the GNSS subnetwork seems to be better than 0.1 mm/yr even in the ITRF2005 and ITRF2008 solutions. This refers of course only to the internal accuracy of the realized origin by the ITRF coordinates (station positions and velocities) in the particular subnetwork, and it should be largely credited to its stronger spatial configuration and higher station density compared to the other technique subnetworks.

Fig. 3
figure 3

Internal accuracy of the realized frame origin in different ITRF solutions and their respective technique subnetworks

Full size image
Fig. 4
figure 4

Internal accuracy of the realized frame scale in different ITRF solutions and their respective technique subnetworks

Full size image
Fig. 5
figure 5

Internal accuracy of the realized frame orientation in different ITRF solutions and their respective technique subnetworks

Full size image

The initial state of the ITRF origin (as realized by the estimated coordinates in each subnetwork at the reference epoch) shows the same internal accuracy in the last two solutions, and also a sizeable improvement since the ITRF2000 release. The accuracy level of 1 mm (at to) is almost reached in ITRF2014 for the DORIS and SLR subnetworks, and partially for the VLBI subnetwork (only in the TZ component). In the GNSS subnetwork, all origin components at the reference epoch are realized with much higher internal accuracy (< 1 mm), even in pre-ITRF2014 solutions, mainly as a result of its denser global coverage over the other technique subnetworks.

The internal accuracy of scale and orientation in the ITRF solutions is plotted in Figs. 4 and 5, respectively. To allow easier reading of these plots, the sigmas for the frame parameters have been expressed in mm and mm/yr after re-scaling by the Earth’s mean radius. Overall, the internal accuracy of the ITRF scale and its rate appears to be weaker in the VLBI subnetwork. This is clearly seen, for example, in ITRF2014 where these parameters are realized at the level of 1.6 mm and 0.3 mm/yr, whereas in the other three subnetworks they exhibit an internal accuracy level that is better than 1 mm and 0.1 mm/yr. On the other hand, the orientation parameters and their rates are weaker in the SLR subnetwork and, to lesser extent, in the DORIS subnetwork (especially for RZ and its rate). The rotation angle RZ and its rate are systematically better in the VLBI subnetwork (compared to SLR and DORIS) for all tested ITRF solutions. Finally, the GNSS subnetwork offers the best realization accuracy for the scale and orientation parameters (and their rates) at a statistical level that is better than 1 mm and 0.1 mm/ yr. It is reminded, however, that this optimistic performance reflects only the internal accuracy as deduced by the full CV matrix of the estimated positions/velocities and it is mainly attributed to the large number of stations within the GNSS subnetwork (see Table 2).

5.3 Appraisal of the long-term accuracy in ITRF solutions

In this section, we examine the frame stability in ITRF solutions by considering the impact of correlations between the estimated positions and velocities in each technique subnetwork. For this purpose, we apply the error covariance expression of Eq. (24) at different intervals δt (up to 20 years), and then the results are shown as continuous function of time in a series of plots given in Fig. 6. These plots depict the square roots of the diagonal elements of the CV matrices \({\varvec{\varSigma}}_{{{\varvec{\theta}}(t_{o} + \delta t)}}\), and they reveal the temporal degradation of the realized coordinate system in each ITRF solution. The zero origin of the time axis in all plots identifies the reference epoch to of each tested solution. For the sake of economy, we purposely omit the assessment of ITRF2005 from the presented results. This does not cause any crucial loss of information since the omitted results are fully consistent with the ones already shown here, while their appearance would only clutter the presentation without any significant gain.

Fig. 6
figure 6

Forward in-time projection of the internal accuracy in ITRF solutions at coordinate system level. The results refer to the four technique subnetworks and the zero origin of the time axis in each plot identifies the reference epoch of the respective solution

Full size image

The long-term frame accuracy is generally improved from one ITRF solution to the next, especially when comparing the earliest and latest releases that are considered here. The only exception occurs in the VLBI subnetwork for the equatorial translations and scale parameters, which seem to retain the same stability level when switching from ITRF2008 to ITRF2014, see Fig. 6. Overall, the ITRF stability does not fully comply with the strict specification of “1 mm + 0.1 mm/yr” which aims to serve the needs of stringent scientific users. According to this specification, the internal accuracy of the ITRF, say after a decade from its reference epoch, should be at the 2 mm level or better (for each frame parameter), yet this is only partially verified by our results. As seen in Fig. 6, none of the pre-ITRF2014 solutions maintains such high stability level in the DORIS, VLBI and SLR subnetworks. The situation is improved in ITRF2014 where only a few frame parameters turn out to be “weakly realized” (at σ > 2 mm) after a period of ten years. These parameters are: the rotation angle RZ (for DORIS), the scale and the equatorial translations (for VLBI), and the three rotation angles RX, RY, RZ (for SLR).

In the VLBI subnetwork, the scale shows the weakest long-term performance among all other frame parameters in all ITRF solutions, whereas in the DORIS subnetwork a similar result occurs with the rotation angle RZ. In the case of the SLR subnetwork, the origin and scale parameters outperform the three rotation parameters in ITRF2005, ITRF2008 and ITRF2014 (but not in ITRF2000). Moreover, if we consider a longer time period, say δt = 15 yrs, then the internal accuracy of the realized scale in the VLBI subnetwork degrades to the level of 4–5 mm, in contrast to the SLR subnetwork where it remains a bit better than 2 mm (this comparison refers to ITRF2014), see Fig. 6.

The long-term accuracy of the frame parameters is significantly better in the GNSS subnetwork (compared to the other technique subnetworks) for all tested ITRF solutions. Specifically, the stability of ITRF2008 and ITRF2014 seems to fulfill the specification “1 mm + 0.1 mm/yr” for all frame parameters, whereas in ITRF2000 it is reached only for the orientation parameters; see Fig. 6. In ITRF2014, the accuracy of the realized origin, orientation and scale in the GNSS subnetwork remains better than 1 mm even after a period of 20 years from the initial reference epoch (2010.0). Again it is noted that such an optimistic assessment reflects only the internal accuracy as deduced by the full CV matrix of the estimated positions/velocities in each ITRF solution, and it is somewhat “biased” by the larger number of stations involved in the GNSS subnetwork.

6 Summary—conclusions

To investigate the accuracy of TRFs at coordinate system level, we introduced a stochastic approach based on a suitable projector for the CV matrix of frame coordinates that are obtained by space geodetic techniques. The theoretical aspects of this evaluation scheme have been addressed in detail in the present paper. The role of the proposed projector is to extract the CMEs (of translational, rotational and dilatational type) that are always hidden in estimated frame coordinates, and thus to deduce the uncertainty of the coordinate system which is realized by these coordinates. Our approach offers a versatile tool to assess the internal accuracy of TRFs and the compliance with application-specific quality standards of their long-term stability on the deforming Earth.

In the context of this study, we have carried out the accuracy evaluation of the realized origin, orientation and scale in the ITRF series (and their technique subnetworks), which revealed the advancements in the achievable frame stability from ITRF2000 up to ITRF2014. The accuracy improvement among the different ITRF solutions that is exposed in our results reflects the growth in quality and time span of space geodetic data, the enhancement of their modeling strategies, and the strengthening of the terrestrial network, that have occurred over the years in the ITRF estimation process. In terms of future work, we plan to repeat the same procedure for the internal accuracy assessment of the latest ITRF2020 solution and also of other available realizations of the International Terrestrial Reference System (i.e., DTRF, JTRF).