# Terrestrial reference frames and their internal accuracy at coordinate system level

Christopher Kotsakis, Miltiadis Chatzinikos

<figure><img src="image_1.png" /></figure>

## To cite this version:

Christopher Kotsakis, Miltiadis Chatzinikos. Terrestrial reference frames and their internal accuracy at coordinate system level. Journal of Geodesy, 2023, 97, 10.1007/s00190-023-01801-6. insu-04475256

HAL Id: insu-04475256

https://insu.hal.science/insu-04475256

Submitted on 23 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Terrestrial reference frames and their internal accuracy at coordinate system level

Christopher Kotsakis, Miltiadis Chatzinikos

<figure><img src="image_1.png" /></figure>

## To cite this version:

Christopher Kotsakis, Miltiadis Chatzinikos. Terrestrial reference frames and their internal accuracy at coordinate system level. Journal of Geodesy, 2023, 97, 10.1007/s00190-023-01801-6. insu-04475256

HAL Id: insu-04475256

https://insu.hal.science/insu-04475256v1

Submitted on 23 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Terrestrial reference frames and their internal accuracy at coordinate system level

Christopher Kotsakis<fnref n="1" /> Ⓒ · Miltiadis Chatzinikos<fnref n="2" /> Ⓒ

Received: 15 February 2023 / Accepted: 25 October 2023 / Published online: 20 November 2023
© The Author(s) 2023

## Abstract

The accuracy assessment of terrestrial reference frames (TRFs) at coordinate system level is a key task to ensure their successful use in Earth studies, satellite navigation and other geodetic positioning applications. Currently, the TRF quality 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. To evaluate the conformity of the internal accuracy of modern TRFs to such requirements, an appropriate mapping is needed to convert frame coordinate errors (and their CV matrix) in a terrestrial network to matching errors (and their CV matrix) in the realized coordinate system. Several projection schemes may be considered for this mapping 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. The goal of the present paper is to investigate the inference problem of frame accuracy at coordinate system level, and to discuss not only the theoretical aspects of the required covariance projectors, but also the practical impact on the results obtained by their implementation in space geodetic solutions. For this purpose, a relevant case study is performed to evaluate the accuracy 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).

**Keywords** Reference frames · Coordinate system · Common-mode errors · Frame accuracy · Covariance projection

## 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

<footnote>Christopher Kotsakis
kotsaki@topo.auth.gr</footnote>

<footnote>Miltiadis Chatzinikos
miltiadis.xatzinikos@obspm.fr</footnote>

<footnote n="1">¹ Department of Geodesy and Surveying, Faculty of Engineering, Aristotle University of Thessaloniki, University Box 440, 54124 Thessaloniki, Greece</footnote>

<footnote n="2">² SYRTE, Observatoire de Paris-PSL, Sorbonne Université, CNRS UMR8630, LNE, 61 avenue de l'Observatoire, 75014 Paris, France</footnote>
# Terrestrial reference frames and their internal accuracy at coordinate system level

Christopher Kotsakis<fnref n="1" /> Ⓒ · Miltiadis Chatzinikos<fnref n="2" /> Ⓒ

Received: 15 February 2023 / Accepted: 25 October 2023 / Published online: 20 November 2023
© The Author(s) 2023

## Abstract

The accuracy assessment of terrestrial reference frames (TRFs) at coordinate system level is a key task to ensure their successful use in Earth studies, satellite navigation and other geodetic positioning applications. Currently, the TRF quality 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. To evaluate the conformity of the internal accuracy of modern TRFs to such requirements, an appropriate mapping is needed to convert frame coordinate errors (and their CV matrix) in a terrestrial network to matching errors (and their CV matrix) in the realized coordinate system. Several projection schemes may be considered for this mapping 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. The goal of the present paper is to investigate the inference problem of frame accuracy at coordinate system level, and to discuss not only the theoretical aspects of the required covariance projectors, but also the practical impact on the results obtained by their implementation in space geodetic solutions. For this purpose, a relevant case study is performed to evaluate the accuracy 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).

**Keywords** Reference frames · Coordinate system · Common-mode errors · Frame accuracy · Covariance projection

## 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

<footnote>Christopher Kotsakis
kotsaki@topo.auth.gr</footnote>

<footnote>Miltiadis Chatzinikos
miltiadis.xatzinikos@obspm.fr</footnote>

<footnote n="1">¹ Department of Geodesy and Surveying, Faculty of Engineering, Aristotle University of Thessaloniki, University Box 440, 54124 Thessaloniki, Greece</footnote>

<footnote n="2">² SYRTE, Observatoire de Paris-PSL, Sorbonne Université, CNRS UMR8630, LNE, 61 avenue de l'Observatoire, 75014 Paris, France</footnote>
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 position-
ing 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 assess-
ing 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 construc-
tion (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 anal-
ysis is still required to deliver a more precise error budget
for the individual frame parameters and their long-term sta-
bility (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 real-
ized 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 solu-
tion by an appropriate mapping of its coordinate CV matrix
onto the Helmert parameter space. Such an approach com-
plements 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 accu-
racy 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 coor-
dinate 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 differ-
ent 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 theoreti-
cal aspects of such projectors, but also the practical impact on
the results obtained by their implementation in space geode-
tic solutions.

The rationale of our work is linked to other previous stud-
ies 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 prin-
ciples are rooted in Meissl's inner error theory which was
introduced in his 1962 paper as a novel framework for filter-
ing 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 net-
works (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 par-
ticular coordinate system, the focus of Meissl's work was to
distinguish these two effects in order to evaluate the geometri-
cal 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
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 position-
ing 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 assess-
ing 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 construc-
tion (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 anal-
ysis is still required to deliver a more precise error budget
for the individual frame parameters and their long-term sta-
bility (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 real-
ized 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 solu-
tion by an appropriate mapping of its coordinate CV matrix
onto the Helmert parameter space. Such an approach com-
plements 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 accu-
racy 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 coor-
dinate 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 differ-
ent 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 theoreti-
cal aspects of such projectors, but also the practical impact on
the results obtained by their implementation in space geode-
tic solutions.

The rationale of our work is linked to other previous stud-
ies 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 prin-
ciples are rooted in Meissl's inner error theory which was
introduced in his 1962 paper as a novel framework for filter-
ing 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 net-
works (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 par-
ticular coordinate system, the focus of Meissl's work was to
distinguish these two effects in order to evaluate the geometri-
cal 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

$$dX_{\theta} = A d\theta \qquad (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\theta$ 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 $dX_{\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

$$dX = \underset{A d\theta}{dX_{\theta}} + d\overline{X} \qquad (2)$$

where $dX$ denotes the random errors relating to the solution's CV matrix $\Sigma_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.

## 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
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

$$dX_{\theta} = A d\theta \qquad (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\theta$ 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 $dX_{\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

$$dX = \underset{A d\theta}{dX_{\theta}} + d\overline{X} \qquad (2)$$

where $dX$ denotes the random errors relating to the solution's CV matrix $\Sigma_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.

## 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
<figure><img src="image_3.png" /><figcaption>**Fig. 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</figcaption></figure>

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\theta$ 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\theta$ 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\theta$ should not be linked with the best (minimum variance) linear unbiased estimator, as this dictates that the error components $dX_\theta$ and $d\bar{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

$$dX = \underbrace{QdX}_{dX_\theta} + \underbrace{(I - Q)dX}_{d\bar{X}} \quad (3)$$
<figure><img src="image_3.png" /><figcaption>**Fig. 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</figcaption></figure>

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\theta$ 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\theta$ 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\theta$ should not be linked with the best (minimum variance) linear unbiased estimator, as this dictates that the error components $dX_\theta$ and $d\bar{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

$$dX = \underbrace{QdX}_{dX_\theta} + \underbrace{(I - Q)dX}_{d\bar{X}} \quad (3)$$
where $Q$ is a projection matrix ($Q^2 = 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

$$Q A = A \qquad (4)$$

which ensures that the CMEs in $dX$ are absorbed by the error component $dX_\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

$$Q d\bar{X} = 0 \qquad (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 $d\bar{X}$ from the underlying coordinate system.

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

$$Q = A (A^T P A)^{-1} A^T P \qquad (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 $d\theta$:

$$\Sigma_{\theta} = (A^{T} P A)^{-1} A^{T} P \Sigma_{X} P A (A^{T} P A)^{-1} \qquad (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

$$A^T P d\bar{X} = 0 \qquad (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

$$P = \Sigma_X^{-1} \quad \text{or} \quad P = I \qquad (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

$$d\bar{X} = (I - A (A^T \Sigma_X^{-1} A)^{-1} A^T \Sigma_X^{-1}) dX \qquad (10a)$$

$$dX_{\theta} = A (A^T \Sigma_X^{-1} A)^{-1} A^T \Sigma_X^{-1} dX \doteq A d\theta \qquad (10b)$$

$$\Sigma_{d\bar{X}, dX_{\theta}} = 0 \qquad (10c)$$

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

$$\Sigma_{\theta} = (A^T \Sigma_X^{-1} A)^{-1} \qquad (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

$$A^T \Sigma_X^{-1} d\bar{X} = 0 \qquad (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

$$d\bar{X} = (I - A (A^T A)^{-1} A^T) dX \qquad (13a)$$

$$dX_{\theta} = A (A^T A)^{-1} A^T dX \doteq A d\theta \qquad (13b)$$

$$\Sigma_{d\bar{X}, dX_{\theta}} = (I - A (A^T A)^{-1} A^T) \Sigma_X A (A^T A)^{-1} A^T \neq 0 \qquad (13c)$$

It is essential to clarify that the choice $P = I$ does not mean that the frame coordinates are assumed to have the same quality or being uncorrelated with each other, but it
where $Q$ is a projection matrix ($Q^2 = 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

$$Q A = A \qquad (4)$$

which ensures that the CMEs in $dX$ are absorbed by the error component $dX_\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

$$Q d\bar{X} = 0 \qquad (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 $d\bar{X}$ from the underlying coordinate system.

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

$$Q = A (A^T P A)^{-1} A^T P \qquad (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 $d\theta$:

$$\Sigma_{\theta} = (A^{T} P A)^{-1} A^{T} P \Sigma_{X} P A (A^{T} P A)^{-1} \qquad (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

$$A^T P d\bar{X} = 0 \qquad (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

$$P = \Sigma_X^{-1} \quad \text{or} \quad P = I \qquad (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

$$d\bar{X} = (I - A (A^T \Sigma_X^{-1} A)^{-1} A^T \Sigma_X^{-1}) dX \qquad (10a)$$

$$dX_{\theta} = A (A^T \Sigma_X^{-1} A)^{-1} A^T \Sigma_X^{-1} dX \doteq A d\theta \qquad (10b)$$

$$\Sigma_{d\bar{X}, dX_{\theta}} = 0 \qquad (10c)$$

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

$$\Sigma_{\theta} = (A^T \Sigma_X^{-1} A)^{-1} \qquad (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

$$A^T \Sigma_X^{-1} d\bar{X} = 0 \qquad (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

$$d\bar{X} = (I - A (A^T A)^{-1} A^T) dX \qquad (13a)$$

$$dX_{\theta} = A (A^T A)^{-1} A^T dX \doteq A d\theta \qquad (13b)$$

$$\Sigma_{d\bar{X}, dX_{\theta}} = (I - A (A^T A)^{-1} A^T) \Sigma_X A (A^T A)^{-1} A^T \neq 0 \qquad (13c)$$

It is essential to clarify that the choice $P = I$ does not mean that the frame coordinates are assumed to have the same quality or being uncorrelated with each other, but it
<table><caption>**Table 1** Summary of projection schemes for the internal assessment of CMEs and the evaluation of frame accuracy at coordinate system level</caption><thead><tr><th></th><th>Projection scheme I</th><th>Projection scheme II</th></tr></thead><tbody><tr><td>Weight matrix for frame coordinates</td><td>P = Σ<sub>X</sub><sup>-1</sup></td><td>P = I</td></tr><tr><td>Projection matrix for frame coordinate errors</td><td>Q = A (A<sup>T</sup> Σ<sub>X</sub><sup>-1</sup> A)<sup>-1</sup> A<sup>T</sup> Σ<sub>X</sub><sup>-1</sup></td><td>Q = A (A<sup>T</sup> A)<sup>-1</sup> A<sup>T</sup></td></tr><tr><td><strong>Frame error projection</strong></td><td></td><td></td></tr><tr><td>Estimation errors of frame coordinates</td><td>dX, Σ<sub>X</sub></td><td>dX, Σ<sub>X</sub></td></tr><tr><td>CMEs</td><td>dX<sub>θ</sub> = QdX = A dθ</td><td>dX<sub>θ</sub> = Q dX = A dθ</td></tr><tr><td>Inner errors</td><td>dX&#773; = (I - Q)dX</td><td>dX&#773; = (I - Q)dX</td></tr><tr><td>Frame errors at coordinate system level</td><td>dθ = (A<sup>T</sup> Σ<sub>X</sub><sup>-1</sup> A)<sup>-1</sup> A<sup>T</sup> Σ<sub>X</sub><sup>-1</sup> dX</td><td>dθ = (A<sup>T</sup> A)<sup>-1</sup> A<sup>T</sup> dX</td></tr><tr><td>Frame accuracy at coordinate system level</td><td>Σ<sub>θ</sub> = (A<sup>T</sup> Σ<sub>X</sub><sup>-1</sup> A)<sup>-1</sup></td><td>Σ<sub>θ</sub> = (A<sup>T</sup> A)<sup>-1</sup> A<sup>T</sup> Σ<sub>X</sub> A (A<sup>T</sup> A)<sup>-1</sup></td></tr><tr><td><strong>Mathematical properties</strong></td><td></td><td></td></tr><tr><td>Nullification of net-translation, net-rotation and net-scaling components (in inner errors)</td><td>No</td><td>Yes</td></tr><tr><td>Stochastic correlation btw CMEs and inner errors</td><td>No</td><td>Yes</td></tr><tr><td>Algebraic independence btw CMEs and inner errors</td><td>dX&#773;<sup>T</sup> Σ<sub>X</sub><sup>-1</sup> dX<sub>θ</sub> = 0</td><td>dX&#773;<sup>T</sup> dX<sub>θ</sub> = 0</td></tr></tbody></table>

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

$$ \Sigma_{\theta} = (A^{T} A)^{-1} A^{T} \Sigma_{X} A (A^{T} A)^{-1} \quad (14) $$

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

$$ A^T d\bar{X} = 0 \quad (15) $$

The above condition ensures the nullification of net-translational, net-rotational and net-scaling components in $d\bar{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

$$ d\bar{X}^T dX_{\theta} = 0 \quad (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.

The projected error components $d\bar{X}$ and $dX_{\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 $\Sigma_X$ before applying the covariance projector of Eq. (14) or Eq. (11).
<table><caption>**Table 1** Summary of projection schemes for the internal assessment of CMEs and the evaluation of frame accuracy at coordinate system level</caption><thead><tr><th></th><th>Projection scheme I</th><th>Projection scheme II</th></tr></thead><tbody><tr><td>Weight matrix for frame coordinates</td><td>P = Σ<sub>X</sub><sup>-1</sup></td><td>P = I</td></tr><tr><td>Projection matrix for frame coordinate errors</td><td>Q = A (A<sup>T</sup> Σ<sub>X</sub><sup>-1</sup> A)<sup>-1</sup> A<sup>T</sup> Σ<sub>X</sub><sup>-1</sup></td><td>Q = A (A<sup>T</sup> A)<sup>-1</sup> A<sup>T</sup></td></tr><tr><td><strong>Frame error projection</strong></td><td></td><td></td></tr><tr><td>Estimation errors of frame coordinates</td><td>dX, Σ<sub>X</sub></td><td>dX, Σ<sub>X</sub></td></tr><tr><td>CMEs</td><td>dX<sub>θ</sub> = QdX = A dθ</td><td>dX<sub>θ</sub> = Q dX = A dθ</td></tr><tr><td>Inner errors</td><td>dX&#773; = (I - Q)dX</td><td>dX&#773; = (I - Q)dX</td></tr><tr><td>Frame errors at coordinate system level</td><td>dθ = (A<sup>T</sup> Σ<sub>X</sub><sup>-1</sup> A)<sup>-1</sup> A<sup>T</sup> Σ<sub>X</sub><sup>-1</sup> dX</td><td>dθ = (A<sup>T</sup> A)<sup>-1</sup> A<sup>T</sup> dX</td></tr><tr><td>Frame accuracy at coordinate system level</td><td>Σ<sub>θ</sub> = (A<sup>T</sup> Σ<sub>X</sub><sup>-1</sup> A)<sup>-1</sup></td><td>Σ<sub>θ</sub> = (A<sup>T</sup> A)<sup>-1</sup> A<sup>T</sup> Σ<sub>X</sub> A (A<sup>T</sup> A)<sup>-1</sup></td></tr><tr><td><strong>Mathematical properties</strong></td><td></td><td></td></tr><tr><td>Nullification of net-translation, net-rotation and net-scaling components (in inner errors)</td><td>No</td><td>Yes</td></tr><tr><td>Stochastic correlation btw CMEs and inner errors</td><td>No</td><td>Yes</td></tr><tr><td>Algebraic independence btw CMEs and inner errors</td><td>dX&#773;<sup>T</sup> Σ<sub>X</sub><sup>-1</sup> dX<sub>θ</sub> = 0</td><td>dX&#773;<sup>T</sup> dX<sub>θ</sub> = 0</td></tr></tbody></table>

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

$$ \Sigma_{\theta} = (A^{T} A)^{-1} A^{T} \Sigma_{X} A (A^{T} A)^{-1} \quad (14) $$

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

$$ A^T d\bar{X} = 0 \quad (15) $$

The above condition ensures the nullification of net-translational, net-rotational and net-scaling components in $d\bar{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

$$ d\bar{X}^T dX_{\theta} = 0 \quad (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.

The projected error components $d\bar{X}$ and $dX_{\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 $\Sigma_X$ before applying the covariance projector of Eq. (14) or Eq. (11).
<figure><img src="image_4.png" /><figcaption>**Fig. 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) $\mathbf{P} = \mathbf{\Sigma}_X^{-1}$ and (b) $\mathbf{P} = \mathbf{I}$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>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)</figcaption></figure>

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 $\mathbf{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 $\mathbf{P} = \mathbf{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 $\mathbf{P} = \mathbf{\Sigma}_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 $\mathbf{\Sigma}_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 $\{\mathbf{\Sigma}_{\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/.

By comparing the results in Fig. 2, it is confirmed that the first option ($\mathbf{P} = \mathbf{\Sigma}_X^{-1}$, left plot) overestimates the frame accuracy compared to the second option ($\mathbf{P} = \mathbf{I}$, right plot)
<figure><img src="image_4.png" /><figcaption>**Fig. 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) $\mathbf{P} = \mathbf{\Sigma}_X^{-1}$ and (b) $\mathbf{P} = \mathbf{I}$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>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)</figcaption></figure>

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 $\mathbf{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 $\mathbf{P} = \mathbf{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 $\mathbf{P} = \mathbf{\Sigma}_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 $\mathbf{\Sigma}_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 {$\mathbf{\Sigma}_{\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/.

By comparing the results in Fig. 2, it is confirmed that the first option ($\mathbf{P} = \mathbf{\Sigma}_X^{-1}$, left plot) overestimates the frame accuracy compared to the second option ($\mathbf{P} = \mathbf{I}$, right plot)
throughout the entire time period. This is more emphati-
cally 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 Σθ 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 strat-
egy 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 expla-
nation 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 esti-
mated frame coordinates that is usually stored in SINEX files,
that is

$$
\Sigma_X = \hat{\sigma}^2 (N + H^T W H)^{-1} \quad (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 **σ̂**² is the a posteriori vari-
ance factor of the estimated solution. The above expression
applies to several (but not all) cases of frame estimation prob-
lems 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 **W**⁻¹
(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 Σθ and the inverse
weight matrix W⁻¹ 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 con-
straints which are typically enforced in a subset of frame

stations. The former matrix is influenced by the latter with-
out 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 fidu-
cial 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

$$
\mathbf{H} = (\mathbf{A}^T \mathbf{A})^{-1} \mathbf{A}^T \qquad (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 ver-
ified by Eq. (14) as follows

$$
\begin{align}
\Sigma_{\theta} &= (A^T A)^{-1} A^T \underbrace{\hat{\sigma}^2 (N + H^T W H)^{-1}}_{\Sigma_X} A (A^T A)^{-1} \nonumber \\
&= \hat{\sigma}^2 (A^T A)^{-1} A^T (N + A(A^T A)^{-1} W (A^T A)^{-1} A^T)^{-1} A (A^T A)^{-1} \nonumber \\
&= \hat{\sigma}^2 W^{-1} \tag{19}
\end{align}
$$

The above result stems from the matrix identity **H**(N + **H**^T**W** **H**)^{-1} **H**^T = **W**^{-1} which holds for any minimally con-
strained normal matrix (see Kotsakis 2012), yet it has limited
interest due to the practical irrelevance of full inner con-
straints 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
throughout the entire time period. This is more emphati-
cally 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 Σ<sub>θ</sub> 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 strat-
egy 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 expla-
nation 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 esti-
mated frame coordinates that is usually stored in SINEX files,
that is

$$
\Sigma_X = \hat{\sigma}^2 (N + H^T W H)^{-1} \quad (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 *σ̂*² is the a posteriori vari-
ance factor of the estimated solution. The above expression
applies to several (but not all) cases of frame estimation prob-
lems 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 *W*⁻¹
(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 Σ<sub>θ</sub> and the inverse
weight matrix W<sup>-1</sup> 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 con-
straints which are typically enforced in a subset of frame

stations. The former matrix is influenced by the latter with-
out 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 fidu-
cial 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

$$
H = (A^T A)^{-1} A^T \qquad (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 ver-
ified by Eq. (14) as follows

$$
\begin{align}
\Sigma_{\theta} &= (A^T A)^{-1} A^T \underbrace{\hat{\sigma}^2 (N + H^T W H)^{-1}}_{\Sigma_X} A (A^T A)^{-1} \nonumber \\
&= \hat{\sigma}^2 (A^T A)^{-1} A^T (N + A(A^T A)^{-1} W (A^T A)^{-1} A^T)^{-1} A (A^T A)^{-1} \nonumber \\
&= \hat{\sigma}^2 W^{-1} \tag{19}
\end{align}
$$

The above result stems from the matrix identity $H(N + H^T W H)^{-1} H^T = 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 $NA = 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

$$X(t) = X(t_o) + (t - t_o) \cdot \dot{X} \quad (20)$$

while their formal CV matrix is analytically expressed as

$$\Sigma_{X(t)} = \Sigma_{X(t_o)} + (t - t_o)^2 \Sigma_{\dot{X}} + (t - t_o) \times \left( \Sigma_{X(t_o), \dot{X}} + \Sigma_{X(t_o), \dot{X}}^T \right) \quad (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

$$\Sigma_X = \begin{bmatrix} \Sigma_{X(t_o)} & \Sigma_{X(t_o), \dot{X}} \\ \Sigma_{\dot{X}, X(t_o)} & \Sigma_{\dot{X}} \end{bmatrix} \xrightarrow{\text{Eq. (14)}} \\ \Sigma_{\theta} = \begin{bmatrix} \Sigma_{\theta(t_o)} & \Sigma_{\theta(t_o), \dot{\theta}} \\ \Sigma_{\dot{\theta}, \theta(t_o)} & \Sigma_{\dot{\theta}} \end{bmatrix} \quad (22)$$

The above result offers a detached representation of TRF accuracy that refers to the initial state (at a reference epoch $t_o$) 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

$$\Sigma_{\theta(t)} = (A^T A)^{-1} A^T \Sigma_{X(t)} A (A^T A)^{-1} \quad (23)$$

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

$$\Sigma_{\theta(t_o+\delta t)} = \Sigma_{\theta(t_o)} + \delta t^2 \Sigma_{\dot{\theta}} + \delta t \left( \Sigma_{\theta(t_o), \dot{\theta}} + \Sigma_{\theta(t_o), \dot{\theta}}^T \right) \quad (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 $\theta(t_o)$ and $\dot{\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 $\theta(t_o)$ and $\dot{\theta}$ (this is largely depended upon the cross-correlations between $X(t_o)$ and $\dot{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
property of the rank-deficient normal matrix, that is $NA = 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

$$X(t) = X(t_o) + (t - t_o) \cdot \dot{X} \quad (20)$$

while their formal CV matrix is analytically expressed as

$$\Sigma_{X(t)} = \Sigma_{X(t_o)} + (t - t_o)^2 \Sigma_{\dot{X}} + (t - t_o) \times \left( \Sigma_{X(t_o), \dot{X}} + \Sigma_{X(t_o), \dot{X}}^T \right) \quad (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

$$\Sigma_X = \begin{bmatrix} \Sigma_{X(t_o)} & \Sigma_{X(t_o), \dot{X}} \\ \Sigma_{\dot{X}, X(t_o)} & \Sigma_{\dot{X}} \end{bmatrix} \xrightarrow{\text{Eq. (14)}} \\ \Sigma_{\theta} = \begin{bmatrix} \Sigma_{\theta(t_o)} & \Sigma_{\theta(t_o), \dot{\theta}} \\ \Sigma_{\dot{\theta}, \theta(t_o)} & \Sigma_{\dot{\theta}} \end{bmatrix} \quad (22)$$

The above result offers a detached representation of TRF accuracy that refers to the initial state (at a reference epoch $t_o$) 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

$$\Sigma_{\theta(t)} = (A^T A)^{-1} A^T \Sigma_{X(t)} A (A^T A)^{-1} \quad (23)$$

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

$$\Sigma_{\theta(t_o+\delta t)} = \Sigma_{\theta(t_o)} + \delta t^2 \Sigma_{\dot{\theta}} + \delta t \left( \Sigma_{\theta(t_o), \dot{\theta}} + \Sigma_{\theta(t_o), \dot{\theta}}^T \right) \quad (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 $\theta(t_o)$ and $\dot{\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 $\theta(t_o)$ and $\dot{\theta}$ (this is largely depended upon the cross-correlations between $X(t_o)$ and $\dot{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
<table><caption>**Table 2** Number of stations in the technique subnetworks that are used for the internal accuracy assessment of the ITRF solutions</caption><thead><tr><th></th><th>ITRF2000</th><th>ITRF2005</th><th>ITRF2008</th><th>ITRF2014</th></tr></thead><tbody><tr><td>DORIS</td><td>78</td><td>113</td><td>130</td><td>155</td></tr><tr><td>SLR</td><td>122</td><td>68</td><td>104</td><td>107</td></tr><tr><td>VLBI</td><td>93</td><td>87</td><td>82</td><td>98</td></tr><tr><td>GNSS</td><td>338</td><td>302</td><td>560</td><td>1054</td></tr></tbody></table>

The numbers shown here do not reflect the station discontinuities that appear in the respective SINEX file of each technique

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 projec- tor of Eq. (14) in conjunction with the Jacobian matrix of the 14-parameter HT model. The reasons for using this pro- jector 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 sys- tem level (i.e., a separate matrix Σθ 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 tech- nique 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 uncer- tainty 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 solu- tions. 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 achiev- able 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 dis- tortions 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 esti- mated 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.

The key aspects to look for in the results of this investiga- tion 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 *inter- nal 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 {Σθ} 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).
<table><caption>**Table 2** Number of stations in the technique subnetworks that are used for the internal accuracy assessment of the ITRF solutions</caption><thead><tr><th></th><th>ITRF2000</th><th>ITRF2005</th><th>ITRF2008</th><th>ITRF2014</th></tr></thead><tbody><tr><td>DORIS</td><td>78</td><td>113</td><td>130</td><td>155</td></tr><tr><td>SLR</td><td>122</td><td>68</td><td>104</td><td>107</td></tr><tr><td>VLBI</td><td>93</td><td>87</td><td>82</td><td>98</td></tr><tr><td>GNSS</td><td>338</td><td>302</td><td>560</td><td>1054</td></tr></tbody></table>

The numbers shown here do not reflect the station discontinuities that appear in the respective SINEX file of each technique

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 $\Sigma_{\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.

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 *inter- nal 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 {$\Sigma_{\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).
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Internal accuracy of the realized frame origin in different ITRF solutions and their respective technique subnetworks</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 4** Internal accuracy of the realized frame scale in different ITRF solutions and their respective technique subnetworks</figcaption></figure>

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.

The initial state of the ITRF origin (as realized by the estimated coordinates in each subnetwork at the reference
<figure><img src="image_6.png" /><figcaption>Fig. 3 Internal accuracy of the realized frame origin in different ITRF solutions and their respective technique subnetworks</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4 Internal accuracy of the realized frame scale in different ITRF solutions and their respective technique subnetworks</figcaption></figure>

As far as the frame origin, the growth of its internal accu-
racy 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 ori-
entation parameters shown in Figs. 4 and 5. A reasonable
explanation is that the increasing time span of space geode-
tic 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 inter-
nal accuracy around 0.1 mm/yr (or better), whereas in the

VLBI subnetwork only the *TZ* rate conforms to this accu-
racy 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.

The initial state of the ITRF origin (as realized by the
estimated coordinates in each subnetwork at the reference
<figure><img src="image_8.png" /><figcaption>Fig. 5 Internal accuracy of the realized frame orientation in different ITRF solutions and their respective technique subnetworks</figcaption></figure>

epoch) shows the same internal accuracy in the last two solu-
tions, and also a sizeable improvement since the ITRF2000
release. The accuracy level of 1 mm (at t_o) is almost reached
in ITRF2014 for the DORIS and SLR subnetworks, and par-
tially for the VLBI subnetwork (only in the TZ component).
In the GNSS subnetwork, all origin components at the refer-
ence 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 sub-
networks.

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 param-
eters 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 parame-
ters 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 system-
atically 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 solu-
tions by considering the impact of correlations between the
estimated positions and velocities in each technique sub-
network. 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
Σθ(t_o+δ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 t_o 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
<figure><img src="image_8.png" /><figcaption>Fig. 5 Internal accuracy of the realized frame orientation in different ITRF solutions and their respective technique subnetworks</figcaption></figure>

epoch) shows the same internal accuracy in the last two solu-
tions, and also a sizeable improvement since the ITRF2000
release. The accuracy level of 1 mm (at t_o) is almost reached
in ITRF2014 for the DORIS and SLR subnetworks, and par-
tially for the VLBI subnetwork (only in the TZ component).
In the GNSS subnetwork, all origin components at the refer-
ence 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 sub-
networks.

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 param-
eters 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 parame-
ters 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 system-
atically 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 solu-
tions by considering the impact of correlations between the
estimated positions and velocities in each technique sub-
network. 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
Σθ(t_o+δ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 t_o 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
<figure><img src="image_9.png" /><figcaption>**Fig. 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</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 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</figcaption></figure>
the ones already shown here, while their appearance would
only clutter the presentation without any significant gain.

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 equato-
rial 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 sci-
entific 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 sub-
networks. 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 $\delta 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 sig-
nificantly 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 accu-
racy 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 suit-
able 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 projec-
tor 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 coor-
dinate system which is realized by these coordinates. Our
approach offers a versatile tool to assess the internal accu-
racy of TRFs and the compliance with application-specific
quality standards of their long-term stability on the deform-
ing Earth.

In the context of this study, we have carried out the accu-
racy 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 Inter-
national Terrestrial Reference System (i.e., DTRF, JTRF).

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s00190-023-01801-6.

**Acknowledgements** The Associate Editor and the three anonymous reviewers are gratefully acknowledged for their valuable comments and suggestions which have improved the final version of the paper.

**Author contributions** CK proposed the main idea of this contribution and designed the general structure of the paper. MC wrote the appropriate code and performed all numerical computations for the frame accuracy assessment in the EPN weekly solutions and the ITRF series. Both authors analyzed the results and worked on their classification and graphical representation. CK wrote the main bulk of the paper and MC reviewed and commented regularly on its contents. Both authors edited the final version of the manuscript.

**Funding** Open access funding provided by HEAL-Link Greece.

**Data availability** The authors declare that the data supporting the find-
ings of this study are available within the article (and its electronic
supplementary material). These datasets were derived from the follow-
ing public domain resources: https://igs.bkg.bund.de/root_ftp/EUREF/
products/, https://itrf.ign.fr/en/solutions/itrf2014#itrf2014-files, https://
itrf.ign.fr/en/solutions/itrf2008#itrf2008-files, https://itrf.ign.fr/en/solu-
tions/itrf2005#itrf2005-files, https://itrf.ign.fr/en/solutions/itrf2000#it
rf2000-files.
the ones already shown here, while their appearance would
only clutter the presentation without any significant gain.

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 equato-
rial 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 sci-
entific 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 sub-
networks. 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 $\delta 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 sig-
nificantly 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 accu-
racy 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 suit-
able 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 projec-
tor 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 coor-
dinate system which is realized by these coordinates. Our
approach offers a versatile tool to assess the internal accu-
racy of TRFs and the compliance with application-specific
quality standards of their long-term stability on the deform-
ing Earth.

In the context of this study, we have carried out the accu-
racy 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 Inter-
national Terrestrial Reference System (i.e., DTRF, JTRF).

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s00190-023-01801-6.

**Acknowledgements** The Associate Editor and the three anonymous reviewers are gratefully acknowledged for their valuable comments and suggestions which have improved the final version of the paper.

**Author contributions** CK proposed the main idea of this contribution and designed the general structure of the paper. MC wrote the appropriate code and performed all numerical computations for the frame accuracy assessment in the EPN weekly solutions and the ITRF series. Both authors analyzed the results and worked on their classification and graphical representation. CK wrote the main bulk of the paper and MC reviewed and commented regularly on its contents. Both authors edited the final version of the manuscript.

**Funding** Open access funding provided by HEAL-Link Greece.

**Data availability** The authors declare that the data supporting the find-
ings of this study are available within the article (and its electronic
supplementary material). These datasets were derived from the follow-
ing public domain resources: https://igs.bkg.bund.de/root_ftp/EUREF/
products/, https://itrf.ign.fr/en/solutions/itrf2014#itrf2014-files, https://
itrf.ign.fr/en/solutions/itrf2008#itrf2008-files, https://itrf.ign.fr/en/solu-
tions/itrf2005#itrf2005-files, https://itrf.ign.fr/en/solutions/itrf2000#it
rf2000-files.
## Declarations

**Conflict of interest** The authors have no conflicts of interest that are relevant to the content of this article.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Altamimi Z, Sillard P, Boucher C (2002) ITRF2000: a new release of the international terrestrial reference frame for earth science applications. J Geophys Res 107(B10):2214

Altamimi Z, Collilieux X, Legrand J, Garayt B, Boucher C (2007) ITRF2005: a new release of the international terrestrial reference frame based on time series of station positions and earth orientation parameters. J Geophys Res 112:B09401

Altamimi Z, Collilieux X, Métivier L (2011) ITRF2008: an improved solution of the international terrestrial reference frame. J Geod 85:457–473

Altamimi Z, Rebischung P, Metivier L, Collilieux X (2016) ITRF2014: a new release of the international terrestrial reference frame modeling nonlinear station motions. J Geophys Res Solid Earth. https://doi.org/10.1002/2016JB013098

Altamimi Z, Rebischung P, Collilieux X, Metivier L, Chanard K (2023) ITRF2020: an augmented reference frame refining the modeling of nonlinear station motions. J Geod 97:47. https://doi.org/10.1007/s00190-023-01738-w

Altamimi Z, Collilieux X, Boucher C (2008) Accuracy assessment of the ITRF datum definition. In: Xu P, Liu J, Dermanis A (eds) VI Hotine-Marussi symposium on theoretical and computational geodesy. International Association of Geodesy Symposia, vol 132. Springer, Berlin. https://doi.org/10.1007/978-3-540-74584-6_16

Baarda W (1973) S-Transformations and criterion matrices. Netherlands Geodetic Commission, Publications on Geodesy (New Series), vol 5. Delft, p 1

Blaha G (1982) Free networks: minimum norm solution as obtained by the inner adjustment constraint method. Bull Geod 56:209–219

Blaha G (1971) Inner adjustment constraints with emphasis on range observations. Department of Geodetic Science, The Ohio State University, OSU Report No. 148, Columbus, Ohio

Blewitt G, Heflin MB, Webb FH, Lindqwister UJ, Malla RP (1992) Global coordinates with centimeter accuracy in the international terrestrial reference frame using the global positioning system. Geophy Res Lett 19:853–856

Blewitt G (1998) GPS data processing methodology: from theory to applications. In: GPS for geodesy, Springer, pp 231–270

Blewitt G (2015) Terrestrial reference frame requirements for studies of geodynamics and climate change. In: van Dam T (eds) REFAG 2014. International Association of Geodesy Symposia, vol 146. Springer, Cham, pp 209–216. https://doi.org/10.1007/1345_2015_142

Collilieux X, Altamimi Z (2013) External evaluation of the origin and scale of the international terrestrial reference frame. In: Altamimi Z, Collilieux X (eds) Reference frames for applications in geosciences. International Association of Geodesy Symposia, vol 138. Springer, Berlin. https://doi.org/10.1007/978-3-642-32998-2_5

Collilieux X et al (2014) External evaluation of the terrestrial reference frame: report of the task force of the IAG sub-commission 1.2. In: Rizos C, Willis P (eds) Earth on the edge: science for a sustainable planet. International Association of Geodesy Symposia, vol 139. Springer, Berlin. https://doi.org/10.1007/978-3-642-37222-3_25

Ebner H (1974) Eine Theorie zur Analyse von Kovarianzmatrizen. Zeitschrift fur Vermessungswesen, pp 453–461

Ebner H (1975) Analysis of covariance matrices. Deutsche Geodätische Kommission, Reihe B, Heft Nr 214, München, pp 111–121

Kotsakis C (2012) Reference frame stability and nonlinear distortion in minimum-constrained network adjustment. J Geod 86:755–774. https://doi.org/10.1007/s00190-012-0555-6

Kotsakis C, Chatzinikos M (2022) Frame accuracy of combined EPN weekly coordinate solutions. In: International Association of Geodesy Symposia. Springer, Berlin. https://doi.org/10.1007/1345_2022_143

Meissl P (1965) Über die innere Genauigkeit dreidimensionaler Punkthaufens. Zeitschrift Fur Vermessungswesen 90:109–118

Meissl P (1962) Die innere Genauigkeit eines Punkthaufens. Österreichische Zeitschrift für Vermessungswesen 50, Nr. 5, S. 159–165, Nr. 6, S, pp 186–194

Meissl P (1969) Zusammengfassung und Ausbau der inneren Fehlertheorische eines Punkthaufens. Deutsche Geodätische Kommission, Reihe A, Heft Nt 61, München, pp 8–21

Plag HP, Pearlman M (eds) (2009) Global geodetic observing system—meeting the requirements of a global society on a changing planet in 2020. Springer Berlin, p 332. https://doi.org/10.1007/978-3-642-02687-4

Pope AJ (1973) The use of the “solution space” in the analysis of geodetic network adjustments. Presented at the IAG symposium on computational methods in geometric geodesy, Oxford, September 2–8, 1973

Pope AJ (1974) Transformation of covariance matrices due to changes in minimal control. AGU Fall Meeting, San Francisco, December 9, 1971

Rebischung P (2014) Can GNSS contribute to improving the ITRF definition? PhD Thesis, Institut National de l' Information Geographique et Forestiere, IGN/LAREG, Paris, France

Sillard P, Boucher C (2001) A review of algebraic constraints in terrestrial reference frame datum definition. J Geod 75:63–73

Wu X, Collilieux X, Altamimi Z, Vermeersen BLA, Gross RS, Fukumori I (2011) Accuracy of the international terrestrial reference frame origin and earth expansion. Geophys Res Lett 38:L13304. https://doi.org/10.1029/2011GL047450
**Declarations**

**Conflict of interest** The authors have no conflicts of interest that are relevant to the content of this article.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

**References**

Altamimi Z, Sillard P, Boucher C (2002) ITRF2000: a new release of the international terrestrial reference frame for earth science applications. J Geophys Res 107(B10):2214

Altamimi Z, Collilieux X, Legrand J, Garayt B, Boucher C (2007) ITRF2005: a new release of the international terrestrial reference frame based on time series of station positions and earth orientation parameters. J Geophys Res 112:B09401

Altamimi Z, Collilieux X, Métivier L (2011) ITRF2008: an improved solution of the international terrestrial reference frame. J Geod 85:457–473

Altamimi Z, Rebischung P, Metivier L, Collilieux X (2016) ITRF2014: a new release of the international terrestrial reference frame modeling nonlinear station motions. J Geophys Res Solid Earth. https://doi.org/10.1002/2016JB013098

Altamimi Z, Rebischung P, Collilieux X, Metivier L, Chanard K (2023) ITRF2020: an augmented reference frame refining the modeling of nonlinear station motions. J Geod 97:47. https://doi.org/10.1007/s00190-023-01738-w

Altamimi Z, Collilieux X, Boucher C (2008) Accuracy assessment of the ITRF datum definition. In: Xu P, Liu J, Dermanis A (eds) VI Hotine-Marussi symposium on theoretical and computational geodesy. International Association of Geodesy Symposia, vol 132. Springer, Berlin. https://doi.org/10.1007/978-3-540-74584-6_16

Baarda W (1973) S-Transformations and criterion matrices. Netherlands Geodetic Commission, Publications on Geodesy (New Series), vol 5. Delft, p 1

Blaha G (1982) Free networks: minimum norm solution as obtained by the inner adjustment constraint method. Bull Geod 56:209–219

Blaha G (1971) Inner adjustment constraints with emphasis on range observations. Department of Geodetic Science, The Ohio State University, OSU Report No. 148, Columbus, Ohio

Blewitt G, Heflin MB, Webb FH, Lindqwister UJ, Malla RP (1992) Global coordinates with centimeter accuracy in the international terrestrial reference frame using the global positioning system. Geophy Res Lett 19:853–856

Blewitt G (1998) GPS data processing methodology: from theory to applications. In: GPS for geodesy, Springer, pp 231–270

Blewitt G (2015) Terrestrial reference frame requirements for studies of geodynamics and climate change. In: van Dam T (eds) REFAG 2014. International Association of Geodesy Symposia, vol 146. Springer, Cham, pp 209–216. https://doi.org/10.1007/1345_2015_142

Collilieux X, Altamimi Z (2013) External evaluation of the origin and scale of the international terrestrial reference frame. In: Altamimi Z, Collilieux X (eds) Reference frames for applications in geosciences. International Association of Geodesy Symposia, vol 138. Springer, Berlin. https://doi.org/10.1007/978-3-642-32998-2_5

Collilieux X et al (2014) External evaluation of the terrestrial reference frame: report of the task force of the IAG sub-commission 1.2. In: Rizos C, Willis P (eds) Earth on the edge: science for a sustainable planet. International Association of Geodesy Symposia, vol 139. Springer, Berlin. https://doi.org/10.1007/978-3-642-37222-3_25

Ebner H (1974) Eine Theorie zur Analyse von Kovarianzmatrizen. Zeitschrift fur Vermessungswesen, pp 453–461

Ebner H (1975) Analysis of covariance matrices. Deutsche Geodätische Kommission, Reihe B, Heft Nr 214, München, pp 111–121

Kotsakis C (2012) Reference frame stability and nonlinear distortion in minimum-constrained network adjustment. J Geod 86:755–774. https://doi.org/10.1007/s00190-012-0555-6

Kotsakis C, Chatzinikos M (2022) Frame accuracy of combined EPN weekly coordinate solutions. In: International Association of Geodesy Symposia. Springer, Berlin. https://doi.org/10.1007/1345_2022_143

Meissl P (1965) Über die innere Genauigkeit dreidimensionaler Punkthaufens. Zeitschrift Fur Vermessungswesen 90:109–118

Meissl P (1962) Die innere Genauigkeit eines Punkthaufens. Österreichische Zeitschrift für Vermessungswesen 50, Nr. 5, S. 159–165, Nr. 6, S, pp 186–194

Meissl P (1969) Zusammengfassung und Ausbau der inneren Fehlertheorische eines Punkthaufens. Deutsche Geodätische Kommission, Reihe A, Heft Nt 61, München, pp 8–21

Plag HP, Pearlman M (eds) (2009) Global geodetic observing system—meeting the requirements of a global society on a changing planet in 2020. Springer Berlin, p 332. https://doi.org/10.1007/978-3-642-02687-4

Pope AJ (1973) The use of the “solution space” in the analysis of geodetic network adjustments. Presented at the IAG symposium on computational methods in geometric geodesy, Oxford, September 2–8, 1973

Pope AJ (1974) Transformation of covariance matrices due to changes in minimal control. AGU Fall Meeting, San Francisco, December 9, 1971

Rebischung P (2014) Can GNSS contribute to improving the ITRF definition? PhD Thesis, Institut National de l' Information Geographique et Forestiere, IGN/LAREG, Paris, France

Sillard P, Boucher C (2001) A review of algebraic constraints in terrestrial reference frame datum definition. J Geod 75:63–73

Wu X, Collilieux X, Altamimi Z, Vermeersen BLA, Gross RS, Fukumori I (2011) Accuracy of the international terrestrial reference frame origin and earth expansion. Geophys Res Lett 38:L13304. https://doi.org/10.1029/2011GL047450