A sequential estimation approach to terrestrial reference frame determination

Available online at www.sciencedirect.com

ScienceDirect Advances in Space Research xxx (2019) xxx–xxx www.elsevier.com/locate/asr

A sequential estimation approach to terrestrial reference frame determination Claudio Abbondanza ⇑, Toshio M. Chin, Richard S. Gross, Michael B. Heflin, Jay W. Parker, Benedikt S. Soja, Xiaoping Wu Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA Received 11 November 2019; accepted 13 November 2019

Abstract We review the main concepts underlying the determination of terrestrial reference frames (TRFs) through a recursive algorithm based on Kalman Filtering and Rauch-Tung-Striebel (RTS) smoothing which is currently adopted at Jet Propulsion Laboratory (JPL) to compute sub-secular frame products (JTRFs). We contextualize the TRF determination in the state-space framework and we emphasize connections between frame state, its observability through space-geodetic frame inputs and the similarity transformation which is central to frame definition. We elaborate on the notion of sub-secular frame, enabled by our approach, in constrast to standard TRF products which, secular by construction, are designed to represent the long-term mean physical properties of the frame. Comparisons of JTRF solutions to standard products such as the International Terrestrial Reference Frame (ITRF) suggest high-level consistency in a longterm sense with time derivatives of the Helmert transformation parameters connecting the two TRFs below 0:18 mm/yr. We discuss advantages and limitations of JPL approach to TRF determination and outline lines of inquiries that are currently being researched as part of JTRF development plan. Ó 2019 Published by Elsevier Ltd on behalf of COSPAR.

Keywords: Terrestrial reference frames; GNSS; VLBI; SLR; DORIS; Kalman filtering

1. Introduction The sequential estimation of global Terrestrial Reference Frames (TRFs) from space-geodetic (SG) products is rooted in state-space theory. Optimal state estimation is an expansively researched discipline applicable to all areas of science and technology concerned with the mathematical modelling of a dynamical system and its evolution over time. In contrast to frequency domain analyses based on Fourier and Laplace transforms, state-space theory adopts a time-domain framework which is particularly fit to characterize the behaviour of time-varying systems.

⇑ Corresponding author.

E-mail address: [email protected] (C. Abbondanza).

Among the first forms of discrete-time mean-square recursive filters developed, Kalman filtering (Ka´lma´n, 1960) is one of the elective approaches in recursive optimal state estimation. Widely applied to a variety of disciplines and applications ranging from oceanography (Fukumori and Malanotte-Rizzoli, 1995), global sea level rise (Hay et al., 2013), precise satellite orbit determination (Lichten and Border, 1987), GPS positioning (Zumberge et al., 1997), VLBI data analyses (Herring et al., 1990), to Earth rotation (Gross et al., 1998), combination of loosely constrained positions inferred from SG and terrestrial geodesy (Dong et al., 1998), and studies of correlations of the radial component of position time series at SG co-located sites (Collilieux et al., 2007), Kalman filtering allows for optimal state estimation of a dynamical system by assimilating noisy observations when an adequate stochastic descrip-

https://doi.org/10.1016/j.asr.2019.11.016 0273-1177/Ó 2019 Published by Elsevier Ltd on behalf of COSPAR.

Please cite this article as: C. Abbondanza, T. M. Chin, R. S. Gross et al., A sequential estimation approach to terrestrial reference frame determination, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.11.016

2

C. Abbondanza et al. / Advances in Space Research xxx (2019) xxx–xxx

tion of the system is supplied. During the data assimilation, Kalman filtering sequentially modifies the set of variables describing the state of a dynamical system by minimizing misfits between what is observed and to what the model physics predicts. From the conceptual standpoint, the state of a system embeds the set of variables providing a complete representation of its status at a given epoch and allowing for sufficiently accurate prediction of its future behaviour. Typical examples in satellite geodesy and navigation are the variables describing the state of an Earth orbiter/spacecraft, such as its position, velocity, acceleration, and angular orientation (see e.g. Tapley et al., 2004). The state in a TRF determination problem has to enable an appropriate description of the complex –and essentially non-linear– spatiotemporal patterns of the Earth deformation, while ensuring that a consistent frame definition is realized and maintained over time. This manuscript contextualizes the TRF estimation problem in Space Geodesy within the conceptual framework of state-space theory and reviews the key notions underlying the sequential frame determination through Kalman filtering adopted by Jet Propulsion Laboratory (JPL) to generate JPL frame products (JTRF, see e.g. Wu et al., 2015; Abbondanza et al., 2017). Section 2 articulates the problem of TRF determination in a sequential fashion by emphasizing the formulation of an appropriate state-space model as well as the central role of the similarity transformation within the combination scheme and briefly recalling the main concepts of Kalman filter and RTS smoother. In Section 3, some of the distinctive properties of JPL frames are illustrated with examples drawn from JTRF2014 – the official JPL product based on the analysis of the ITRF2014 inputs– and discussed in relation to ITRF (Altamimi et al., 2016) and DTRF (Seitz et al., 2016). Section 4 outlines the JTRF development and improvement plan as efforts to construct accurate JPL frame products continue. 2. Algorithm description Building on the methodology outlined in Wu et al. (2015), JPL constructs frame products by assimilating SG inputs and local tie observations with a discrete-time Kalman filter and Rauch-Tung-Striebel (RTS) smoother. 2.1. Basic notions If xi denotes the set of variables describing the frame state at time index i and zi the observations of the dynamical system (deforming solid Earth and its variable rotation, in this instance), an optimal solution for the sequential problem can be formalized as follows: min

x1 ;x2 ;...;xk

k n X

2

2

kxi  Uxi1 kQ1 þ kzi  H i xi kR1

i¼1

i

i

o

ð1Þ

where the propagator U encodes the dynamic structure of the system and allows the state xi to be propagated forward in time, H i is the linear operator mapping the frame state into the observation space, Qi is the covariance of the state process noise, Ri is the measurement covariance, and k  kA denotes the ‘2 -vector norm generalized according to the symmetric and positive-definite matrix A. The optimal solution of Eq. (1) can be recursively determined by applying a Kalman filter algorithm to:  xk ¼ U xk1 þ wk ð2Þ zk ¼ H xk þ vk where wk and vk are zero-mean, mutually uncorrelated Gaussian processes with covariance matrices Qk and Rk . Kalman filtering computes the optimal estimate of the pair ðxk ; P k Þ based upon the recursive sequence of a time update step, during which the frame state and its covariance (P k ) are extrapolated forward in time (i.e. predicted): ( xkjk1 ¼ U xk1jk1 ð3Þ P kjk1 ¼ P k1jk1 U P tk1jk1 þ Qk followed by a measurement update, during which the predicted state (xkjk1 ) and its covariance (P kjk1 ) get optimally updated based on new sets of measurements available: 8  1 > K k ¼ P kjk1 H tk H k P kjk1 H tk þ Rk > > > < yk ¼ zk  H k xkjk1 ð4Þ > xkjk ¼ xkjk1 þ K k y k > > > : P kjk ¼ ðI  K k H k ÞP kjk1 xkjk1 (P kjk1 ) in Eqs. (3) and (4) indicates the prior estimate of x (P) at time k based on the measurements z available up to time index k  1, whereas xkjk (P kjk ) in Eq. (4) denotes the optimal posterior estimate of x (P) obtained by assimilating observations available up to index k; the symbol ðÞt denotes the transpose operator. y k is the innovation conveyed by the new measurements zk , which directly affect the updated state xkjk through the gain matrix K k . A comprehensive treatment of the analytical developments for discrete dynamical systems leading to Equations 3 and 4 can be found in e.g. Gelb (1974), Stengel (1994), Simon (2006), Bierman (2006). Optimal state estimates based on the entire set of available observations are achieved through smoothing which, in its most elemental form, entails running the filter backward in time and combining forward and backward estimates. Rather than adopt forward-backward filtering, JPL makes use of a fixed time-step RTS smoother (see e.g. Chapter 9 of Simon, 2006). 2.2. State-space model and time update To contextualize the description of Section 2.1 in the framework of TRF analysis, an appropriate state-space model is required. The TRF state definition must allow for the non-linear nature of the solid Earth deformation

Please cite this article as: C. Abbondanza, T. M. Chin, R. S. Gross et al., A sequential estimation approach to terrestrial reference frame determination, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.11.016

C. Abbondanza et al. / Advances in Space Research xxx (2019) xxx–xxx

while ensuring that a consistent frame definition is realized and maintained over time. This is achieved by including in the state vector entities descriptive of the Earth deformational state along with the parameters of a similarity transformation used to relate each single-technique frame input to the combined output. Drawing from the theoretical framework outlined in Chapter 4 of the International Earth Rotations and Reference Systems Service (IERS) Conventions (Petit and Luzum, 2010), JPL constructs a sequential approach to TRF determination grounded on the notion of similarity, a shape-preserving transformation that can be used to relate single-technique frame inputs to combined frame outputs without inducing, in principle, any deformation of the underlying geometric structures of the observables (see e.g. Sillard and Boucher, 2001). Two sets xi and x0i of three-dimensional Cartesian coordinates in two distinct Euclidean affine spaces (i.e. frames), ideally describing an equivalent Earth surface deformational state, are functionally related through a similarity transformation: x0i ¼ xi þ bi þ kxi þ eijk rj xk

ð5Þ

where bi denotes the i-th component of the translation vector, k defines the distance ratio (i.e. scale), ri represents the infinitesimal clockwise rotation angle about the coordinate axis i and eijk is the Levi-Civita permutation symbol (repeated indices in the products imply summation). Eq. (5) is derived by linearising the three-dimensional rotation and scale transformation and is therefore valid under the assumption of smallness of k and ri (Altamimi et al., 2002). In TRF literature, Eq. (5) is often referred to as Helmert transformation and the set ðbk ; k; rk Þ as Helmert transformation parameters (see e.g. Collilieux et al., 2010). Frame inputs commonly utilized in TRF reductions provide time-dependent observations of Earth rotational and deformative state at discrete points on its surface as determined through global networks of Very-Long Baseline Interferometry (VLBI), Global Navigations and Satellite Systems (GNSS), Satellite Laser Ranging (SLR), and Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS). Solution INdependent EXchange Format (SINEX)1 files are conventionally adopted to exchange frame data and contain integrated station positions (reported as Cartesian coordinates) at different temporal resolutions ranging from weekly to daily intervals, as well as daily Earth Orientation Parameters (EOPs) (see e.g. Gross, 2015; Dehant et al., 2015) along with their full covariance. The state-space description of the dynamic system (cf Eq. (2)) adopted in JPL frame analyses can be formalized as: 1

A detailed description of SINEX format file, officially endorsed by IERS, is accessible at the following wepbage: https://www.iers.org/IERS/ EN/Organization/AnalysisCoordinator/SinexFormat/sinex.html.

2

Xk

3

2

7 6 6Vk 7 7 6 6Fk 7 ¼ 7 6 7 6 4 Hk 5 Ek |fflfflffl{zfflfflffl} xk

I

Dt

0

I 0

0 F

3

0 0

3 2

X k1

3

2

wxk

3

7 6 7 6 7 0 0 7 6 V k1 7 6 wvk 7 7 6 7 6 7 7 6 7 6 0 07 7  6 F k1 7 þ 6 wf k 7 7 6 7 6 7 0 0 I 0 5 4 Hk1 5 4 whk 5 wek 0 0 0 0 I Ek1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflffl{zfflfflfflfflffl} |fflfflffl{zfflfflffl} U xk1 Wk 6 60 6 60 6 6 40

ð6Þ where the elements of the state vector xk are ideally partitioned into (i) station-related parameters such as positions (X k ), linear trends (V k ), Fourier modes (F k ) describing periodic oscillations of the station positions, and (ii) global network parameters, such as Helmert transformation parameters (Hk ) and EOPs (Ek ). The propagator U (cf Eq. (1)) is also referred to, in specialized literature, as state transition matrix. With non-zero elements mostly confined to the diagonal band, the state transition matrix U is sparsely banded. Focussing on the sub-structure of U associated with the station-related parameters, Eq. (6) can be specialized as: 2

xi;k

3

21

Dt

0

6 7 6 0 6 vi;k 7 6 0 1   6 7 6 6 7 6 0 0 2eDts cos 2pDt s 6 i;k 7 6 Ts 6 7¼6 6s 7 6 1 6 i;k1 7 6 0 0 6 7 6 6q 7 6 0 4 i;k 5 4 0 0 qi;k1 0 0 0 2 3 3 2 xi;k1 wxi;k 6 7 7 6 6 vi;k1 7 6 wvi;k 7 6 7 7 6 6 7 7 6 6 si;k1 7 6 wsi;k 7 6 6 7 7 6 þ6 7 7 6 si;k2 7 6 0 7 6 7 7 6 6q 7 7 6 4 i;k1 5 4 wqi;k 5 qi;k2 0

0

0

0

0

e

2Dt s

0

0

0

0   Dt 2e s cos 2pDt Tq

0

1

0

3

7 7 7 0 7 7 7 0 7 7 7 2Dt e s 7 5 0

0

ð7Þ

where  the index i ¼ 1; 2; 3 denotes the Cartesian component in the Euclidean affine space,  k is the time index,  Dt is the fixed time step set to 1 week,  xi;k represents the non-periodic position of the station i at the epoch k,  vi;k is the instantaneous velocity (trend variable),    the triplets ðsi;k ; si;k1 ; si;k2 Þ and qi;k ; qi;k1 ; qi;k2 are the elements of a second-order regression model discretizing the differential equation of a stochastic damped oscillator (Chin et al., 2004): the former relates to the annual periodic mode, the latter to the semi-annual periodic mode;  T s and T q designate the annual and semi-annual periods, respectively.  s, the e-folding time describing the attenuation in the harmonic oscillator, is set to a large value so that no damping is allowed.

Please cite this article as: C. Abbondanza, T. M. Chin, R. S. Gross et al., A sequential estimation approach to terrestrial reference frame determination, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.11.016

4

C. Abbondanza et al. / Advances in Space Research xxx (2019) xxx–xxx

The station-related state parameters ðX k ; V k ; F k Þ allow the Cartesian coordinates to be decomposed into the sum of linear trends and periodic oscillations, the former describing geophysical phenomena such as e.g. plate motions and glacial isostatic adjustments, the latter relating to seasonal oscillations due to surficial mass loading. While several Fourier modes of oscillations can be included in this formulation, only annual and semi-annual modes are being currently adopted in JPL products, since these represent the dominant periodic variations found in the frame inputs (see e.g. Collilieux et al., 2007; Ray et al., 2008). It has to be noted that the Fourier state parameters of Eq. (7) may not only reflect solid Earth loading signals, but also systematic effects characterized by seasonal modulations of disparate nature affecting the frame inputs. This is the case of GNSS data, whose station positions at seasonal frequencies might be affected by, for instance, thermoelastic deformation signatures (see e.g. Xu et al., 2017) and thermally induced monument tilting (see e.g. Yan et al., 2009). The state vector xk is completed by 42 daily EOP parameters describing the Earth rotation (namely xp ; y p ; UT and their time derivatives x_ p ; y_ p ; LOD) for each day within the weekly interval and by Helmert transformation parameters grouped in sets each of which to be individually applied to the technique-specific global SG networks. As will be described in Section 2.3, the role of the similarity transformation parameters is crucial in that they enable each technique-specific frame input to be mapped into the frame output (i.e. combined frame) during the measurement update (cf Eqs. (9)). SG observing systems are only able to sense the total station motion ~xi;k , which can be interpreted as a linear combination of the station-related state parameters: ~xi;k ¼ ‘i;k þ si;k þ qi;k þ wxi;k

ð8Þ

where ‘i;k ‘i;k ¼ xi;k1 þ vi;k1  Dt denotes the non-periodic displacement of the station at epoch k. Eq. (8) accounts for all of the elements into which the station motion has been decomposed and will be utilized during the measurement update described in Section 2.3. The elements of the vector W k of Eq. (6) relate to the stochastic nature of the state vector parameters, whose associated process noise is modeled as white-noise driven random walk. By exploiting the recursive nature of Eq. (6), at the generic epoch n, each of the terms of the process noise vector W k can be represented as sequence of partial sums nn ¼

n X

wj

j¼1

where wj are indipendent, identically distributed Gaussian variables. We note that such an equation provides one of

the commonly adopted definitions for discrete random walk driven by a Gaussian process (see e.g. Section 10.1 in Papoulis and Pillai, 2001). JPL computational code enables process noise to be activated on each of the state vector parameters, although, in current TRF reductions, process noise is solely turned on for station positions, EOPs, and Helmert transformation parameters. By setting the variances of the Gaussian processes wvi ;k , wsi ;k and wqi ;k of Eq. (7) to non-zero values, velocities and seasonal terms are allowed to stochastically vary in time, thus accommodating time-variable trends and amplitudes/phases of the seasonal oscillators (Davis et al., 2012). The filter is initialized by utilizing the a priori values specified in Table 1, along with the uncertainties adopted to fill the diagonal terms of the covariance matrix at the first step of the assimilation. No dynamic is currently utilized to describe similarity and EOP parameters, which are both modelled as pure random walk whose associated r is set to a value sufficiently large to accommodate their observed week-to-week variations (cf Table 1). The sequential estimator allows for flexibility as to the structure of the state vector X k and its process noise W k . For instance, Fourier modes can be selectively removed and the stochastic perturbations to state parameters turned off. When the state vector does not include the oscillators and the stochastic perturbations to non-periodic displacements (wxi;k ) and linear trends (wvi;k ) are not adopted (i.e. the variance associated to each of the elements constituting the vector W k is set to 0), the sequential algorithm degenerates into a recursive least square estimator and piece-wise linear frames can be recovered. The current set-up of JPL frame products make use of linear trend, annual, and semi-annual Fourier modes with process noise on the state parameters as described in Table 1. 2.3. Measurement Update During the measurement update the state forecasts (i.e. predictions) obtained through Eq. (6) are adjusted by minimizing the misfits between the physical model and the station position and EOP observations (cf Equations 4). As mentioned in Section 2.1, the observation equations involved in the measurement update have to establish a one-to-one relationship between frame inputs in the measurement space and the state parameters. The linearized similarity transformation of Eq. (5) may be conveniently used, as outlined in Altamimi et al. (2002), to construct the rank-deficient linear system central to the measurement update: 8 s;P P c;P P P c;P > x ¼ xc;P i þ bi þ k xi þ eijk rj xk > > i > > < xs;R ¼ xc;R þ bR þ kR xc;R þ eijk rR xc;R i i i i j k ð9Þ s;L c;L L c;L L L c;L > ¼ xi þ bi þ k xi þ eijk rj xk > > xi > > : xs;D ¼ xc;D þ bD þ kD xc;D þ e rD xc;D ijk j k i i i i

Please cite this article as: C. Abbondanza, T. M. Chin, R. S. Gross et al., A sequential estimation approach to terrestrial reference frame determination, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.11.016

C. Abbondanza et al. / Advances in Space Research xxx (2019) xxx–xxx

5

Table 1 Set-up adopted for the filter initialization: initial values of the state parameters as well as of the uncertainties adopted in the first step of the filter are reported in Column 3 (a priori state values) and Column 4 (a priori state uncertainties). Process noise for the state parameters is selected to be random walk driven by white noise, whose variance square roots are reported in Column 5. Column 6 defines the units adopted in current JPL frame reduction. Parameter xk vk sk ; qk xp ; y p x_ p ; y_ p UT LOD T k R

Description Non-periodic Coordinates Station Velocities Periodic Displacements Polar Motions Polar Motion Rates Universal Time Excess Length-of-Day Translations Scale Rotations

a priori Value

Uncertainty

Process Noise r

Units

2 2 5 400 400 400 400 2 20 200

b

m m/yr cm mas mas/day ms ms/day m 0:1ppm 0:01  as

a

0 0 0 0 0 0 0 0 0

0 0 100 100 100 100 1 1:6 3:2

a

A priori station coordinates are derived from the first observed station coordinates in the time series. Station-dependent standard deviations of the position process noise are derived from the analysis of geophysical models of loading displacements caused by terrestrial mass transport as described in Wu et al. (2015). b

where  xs;t and xc;t are vectors whose elements are the techniquespecific total coordinates described in Eq. (8) (the superscript ~ has been here dropped for clarity of representation). The index t represents the SG technique where P stands for GNSS, R for VLBI, L for SLR and D for DORIS;  xs;t relates to the observed weekly/daily station positions reported in the frame input files;  xc;t to the weekly combined (outputted) station positions;  bi denotes the vector of the translations mapping the combined frame to the technique-specific xs;t frame;  kt denotes the scale factor between the combined and the technique-specific input frames;  rtj is the infinitesimal angle about the coordinate axis j clockwise rotating the combined frame to the technique-specific input frame. Rather than adopt extended Kalman filter formulations (see e.g. Chapter 13 of Simon, 2006) for the system of Eq. (9) –essentially non-linear in the frame state parameters, JPL chooses to apply standard Kalman filtering to the linearized measurement equations. As a result of the linearization, the optimal state estimates are to be regarded as perturbations about the approximate values of the similarity transformation parameters and the station positions. It should also be noted that, since small station position motions at sub-weekly time scale are being ignored, the time dependency in Eq. (9), otherwise crucial in ITRF reductions (see e.g. Eq. (A9) in Altamimi et al. (2002) and Eq. (1) in Altamimi et al. (2007)) is here dropped. The system of Eq. (9) is complemented by EOP measurement equations relating the combined and observed Earth rotation to the rotational parameters, as documented in Altamimi et al. (2011):

8 s;t xp > > > > > y s;t > p > > > < UT s;t > x_ s;t > p > > > > > y_ s;t p > > : LODs;t

¼ ¼

xcp þ rt;2 y cp þ rt;1

¼ UT c  f1 rt;3 ¼

x_ cp y_ cp

¼

LODc

¼

ð10Þ

where  the subscript t denotes the Earth Orientation (Helmert rotational) parameters as observed by (determined through) the SG technique t;  the superscript c denotes the combined EOP parameters included   in the state vector;  xp ; y p denotes the daily polar motion;    x_ p ; y_ p denotes the daily polar motion rate;  UT denotes the difference ðUT 1  UTC Þ;  LOD is the excess length of day representing discrete variations of UT over one day;  f ¼ 1:002737909350795 is the conversion factor from UT into sidereal time. During the measurement update the technique-specific EOPs are simultaneously calibrated with the TRF rotational parameters so as to compensate for EOP-related systematic biases. It has to be observed that, since sub-weekly variations of the transformation parameters are being in JPL implementation, the dependency of neglected  x_ p ; y_ p and LOD in Eq. (10) upon r_ t;i , as described in Eq. (2) of Altamimi et al. (2007), vanishes. In JPL reductions, polar motion observations from VLBI, SLR, GNSS and DORIS, polar motion rates from VLBI and GNSS, UT and LOD from VLBI are typically assimilated. The satellite-derived LOD, though available, has not been utilized in order to avoid potential corruption of VLBI-

Please cite this article as: C. Abbondanza, T. M. Chin, R. S. Gross et al., A sequential estimation approach to terrestrial reference frame determination, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.11.016

6

C. Abbondanza et al. / Advances in Space Research xxx (2019) xxx–xxx

determined LOD (Ray, 1996). The observations of positions and EOPs fed into the sequential estimator are weighted through the inverse of the measurement (full) covariance matrices reported in the input SINEX files of the four SG solutions. Eqs. (9) and (10) constitute the basis for the observation equations adopted in CATREF (see e.g. Altamimi et al., 2002; Altamimi et al., 2005; Altamimi et al., 2005), the least-square batch estimator through which ITRF products are computed. 2.4. Frame definition To define the origin of the output frame, the SLR translation parameters bLi in the rank-deficient system of Eq. (9) are set to zero during the measurement update. The conditions bLi ¼ 0 guarantee that no translation gets applied to the SLR inputs at each step of the measurement update. Since translations for the remaining SG techniques are not set to zero, the origin of VLBI, GNSS and DORIS networks are conversely being adjusted to SLR quasiinstantaneous center of mass. The output frame scale can be defined by setting the kt parameters; For instance, the conditions kR ¼ 0 ensure that the VLBI scale is not being adjusted during the weekly measurement updates, while scale parameters are estimated for GNSS, SLR, and DORIS, thus implying their intrinsic scale information is removed and adjusted to the instantaneous VLBI scale. In such case, the scale combined frame is solely determined by VLBI. The scale for JTRF frame products was determined by setting both kL ¼ 0 and kR ¼ 0. In principle valid were there no bias between VLBI and SLR scale, such definition would realize a combined TRF whose scale is given by the quasi instantaneous average of the interferometric (VLBI) and optical (SLR) scales weighted through the inverse of the covariance measurement matrices. To define orientation, rotation parameters rti featured in Eq. (9) can be estimated for all of the four SG techniques while applying rotational minimal constraints to an external reference frame as additional observations during the measurement update. Following the framework outlined in Chapter 4 of the IERS Conventions (Petit and Luzum, 2010), soft (to a given r) minimal rotation constraints are formed and applied at each step of the weekly filter:   d ij;k xcj;k  xrj;k ¼ 0 ð11Þ where xcj;k denotes the j-th Cartesian component of the output frame coordinates at time index k, and xrj;k is the reference coordinate extrapolated at time index k. Note that, when the reference coordinates are drawn from predetermined frame realizations, xcj;k is said to inherit, as a result of the rotational constraints, the orientation from the pre-existing frame (Blewitt, 2015). Alternatively, the reference coordinates xrj;k can be constructed by adopting

horizontal velocities inferred from absolute plate motion models (PMMs). xcj;k is therefore made, as a result of the application of Eq. (11), rotationally consistent with the underlying PMM. The time-varying d ij;k terms in Eq. (11) are elements of the projector h 1 i d ij;k ¼ Atk Ak Atk ð12Þ ij

where

2

.. .. 6. . 61 0 6 6 Ak ¼ 6 60 1 6 60 0 4 .. .. . .

.. .

.. .

.. .

0 0

0 xr3;k

xr3;k 0

1

xr2;k

xr1;k

.. .

.. .

.. .

.. .

3

7 xr2;k 7 7 7 r x1;k 7 7 7 0 7 5 .. .

ð13Þ

An alternate approach, based upon the notion of Tisserand reference axes (see e.g. Chapter 3 of Moritz and Mu¨ller, 1987), can be envisaged when defining the orientation of a reference frame. It relies on the definition of internal angular momentum hi (see e.g. Eq. (2) in Altamimi et al., 2003) computed with respect to Earth’s geocenter Z   hi ¼ q eijk xj x_ k dX ð14Þ X

where xj denotes the j-th component of the Cartesian geocentric coordinates, x_ k the k-th component of the geocentric velocities and X, the integration domain, is over Earth’s volume. Nullifying hi in Eq. (14) leads to select the orientation of a set of Tisserand axes. The practical implementation of an explicit no-netrotation (NNR) condition entails defining a NNRcompliant velocity field x_ Ni which is functionally related to the SG-derived velocity field x_ i through the following transformation: x_ Ni ¼ x_ i þ eijk xj xk

ð15Þ

where the angular velocity vector xi of Eq. (15) is estimated in such a way that Z   ð16Þ hi ¼ q eijk xj x_ Nk dX ¼ 0 X

Incidentally, we note that, in TRF literature, the application of Eq. (16) in which the volume integral is substituted by a surface integral computed over Earth’s crust is also referred to as No-Net-Rotation Condition over the Crust (NNRC) (see e.g. Eq. (5) in Altamimi et al., 2002). The implementation of an explicit NNR condition as featured in Eq. (16) may differ based upon the way in which the integral in Eqs. (14) and (16) is represented and discretized. The reader interested on a review on the subject is addressed to the comprehensive work of Legrand (2007). It is important to emphasize that the adoption of the minimal rotational constraints featured in Eq. (11) will lead to a frame product xcj;k satisfying the NNR condition inso-

Please cite this article as: C. Abbondanza, T. M. Chin, R. S. Gross et al., A sequential estimation approach to terrestrial reference frame determination, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.11.016

C. Abbondanza et al. / Advances in Space Research xxx (2019) xxx–xxx

far as the pre-existing frame realization xrj;k is NNRcompliant. Finally, we observe that, notwithstanding minimal rotational constraints and NNRC should both lead in theory to the same result, the two approaches are analytically distinct: the former (Eq. (11)) aims to determine and nullify in a least-square fashion the rotational differences between two frame realizations; The latter determines a velocity field explicitly nullifying the internal angular momentum (see Eq. (16)). 2.5. Local ties and co-motion constraints Local ties are assimilated during the measurement update. Tie vector observations provide linear measurement equations in the form c;m xc;l ¼ Dxlm i  xi i

ð17Þ

where the observables Dxlm i represent differential Cartesian coordinates between co-located SG stations l and m and are extracted from local tie SINEX files. The terms xc;l i and xc;m in Eq. (17) are frame outputs related to the state i parameters through the linear combination of Eq. (8). Tie vectors are applied only once at the epoch of their measurements, as reported in their input SINEX files and are weighted through the inverse of the measurement covariance matrices scaled by properly selected variance factors. Incidentally, we observe that the algorithm, due to its sequential nature, allows multiple local tie estimates for the same site, surveyed and determined at different epochs, to be seamlessly assimilated in the frame determination. In principle the sequential estimator could even assimilate, were they available, time series of local ties. Following standard ITRF combination practice (see e.g. Altamimi et al., 2007, 2011), the optimal set of tie variance factors in JPL frame products is empirically determined as a function of the agreement between SG and terrestrial observations. In order to limit local frame distortions potentially caused when assimilating local ties, tie vectors rather than SG observations are down-weighted. By adopting the state-space model illustrated in Section 2.2, the process noise variances of Table 1 and the additive model of Eq. (8), the total displacement at station l from time k to k þ 1 can be described as ~xli;kþ1  ~xli;k ¼ vli;k  Dt þ sli;kþ1  sli;k þ qli;kþ1  qli;k þ wlxi;k

ð18Þ

Co-motion constraints dictate the equivalence of the displacements of Eq. (18) when observed by two or more nearby co-located stations, and require  The trend variable (velocity) of co-located SG stations l and m be equivalent

7

 The seasonal oscillator state parameters relevant to colocated stations be the same.  The station position process noise wlxi;k and wmxi;k relevant to co-located stations l and m be characterized by a 1:0 correlation coefficient. The equivalence of the station velocities and of the seasonal oscillators for two generic co-located stations l and m is formalized through soft (stochastic) constraints associated with pre-assigned uncertainty: 8 l v  vmi;k ¼ 0; ðrv Þ > < i;k sli;k  smi;k ¼ 0; ðrs Þ ð19Þ >   : l m qi;k  qi;k ¼ 0; rq The constraints of Eq. (19) provide additional measurement equations to be applied during the measurement update. The seasonal oscillator terms (cf Eq. (7)) for colocated stations are constrained to be equivalent with uncertainty value rs;q uniformly set to 10 lm, based on the assumption that, conceivably, nearby stations are subject to and detect the same seasonal displacements. As to the velocity of co-moving stations, the uncertainties rv associated with the equivalence constraints are, unlike what happens for the seasonal oscillators, stationdependent. In JPL frame products, velocity constraintrelated uncertainties are empirically determined by analyzing the time series of station position residuals at colocated sites, based on the assumption that inaccurately constrained velocities tend to produce systematically biased (non-zero mean) residual position time series. Table 2 gathers statistics computed on the values of the stochastic velocity constraints adopted in JTRF2014, thus showing the degree of variability of rv across the entire JTRF network. The statistics are computed by grouping the stochastic constraints by SG technique. Only pairs to GNSS colocated stations have been included in this analysis which does not account for the large number of intra-technique (i.e. GNSS-to-GNSS, DORIS-to-DORIS, and, to a lesser extent, VLBI-to-VLBI and SLR-to-SLR) velocity constraints. For co-located stations showing a remarkable agreement in velocity, rv is in the order of 10 lm/yr. In few exceptional cases where large discrepancies in the velocities of GNSS and VLBI stations are found, rv can be made as large as 10 cm/yr. The largest dispersions for rv relate to VLBI, the smallest to SLR, followed by DORIS. DORIS-to-GNSS median rv is at the mm-level, unlike VLBI and SLR which are characterized by submm median values. The constraints on the station position processes wlxk and wjxk are realized by imposing the correlation term qlm between stations l and m to be 1:0 in the process noise covariance matrix at each step of the forward filter.

Please cite this article as: C. Abbondanza, T. M. Chin, R. S. Gross et al., A sequential estimation approach to terrestrial reference frame determination, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.11.016

8

C. Abbondanza et al. / Advances in Space Research xxx (2019) xxx–xxx

Table 2 Statistics expressing the variability of the stochastic equivalence constraints applied to the velocities of co-located stations. Statistics are computed by grouping the stochastic constraints by SG technique. Only pairs to GNSS co-located stations are being considered. The first row relates to the pairs VLBI (R)-GNSS(P); The second to SLR(L)-GNSS(P); The third to DORIS(D)-GNSS(P). Units for the statistics are [mm/yr]. The Column (N) reports the number of pairs utilised to compute the minimum (min), maximum (max), median values (med), and standard deviations (std). The Column DOMES (Site) indicates the DOMES number of the ITRF site wherein the velocity constraints are applied at a sigma level as reported in Column (max). T

min

max

med

std

N

DOMES (Site)

R L D

0.01 0.01 0.10

100.0 5.0 10.0

0.1 0.5 1.1

15.8 1.1 2.3

39 35 55

13407S010 (Madrid-Robledo, Spain) 13504M002 (Kootwijk, Netherlands) 41705S007 (Santiago, Chile)

3. Properties The main features of JPL frame products are here listed and discussed in relation to ITRF and DTRF (Seitz et al., 2012) products (cf Table 3). 3.0.1. (1) JTRF is sub-secular The qualifier sub-secular, used, for instance, in Earth rotation studies to indicate the broad-range spectral band including decadal, seasonal, inter/intra annual, and higher-frequency portions of the polar motion spectrum, here denotes the timescale of JPL products. Such finer timescale (i.e. sub-secular) set JTRF apart from ITRF and DTRF which, secular by construction, are designed to capture and represent the long-term (i.e. secular) mean physical properties of the frame. JPL products are constructed in such a way that the frame origin is set at the quasi-instantaneous center of mass as sensed by SLR and the scale is the quasi-instantaneous average of VLBI and SLR scale. JPL sequential approach to TRF determination makes use of local tie observations to spatially transfer the SLR

quasi-instantaneous center of mass to the SG techniques either insensitive to l ¼ GM  , such as VLBI, (where G is the universal Gravitational constant and M  is Earth’s mass) or for which geocenter motion is poorly determined (DORIS and GNSS) as well as the averaged VLBI/SLR instantaneous scale to GNSS and DORIS. Likewise, the adoption of co-motion constraints along with the estimation at each step of the filter of the similarity parameters (cf Eq. (9)) ensures that the instantaneous origin and scale information are properly transferred in time. As a result of the sub-secular timescale of JPL frames, JTRF gives users access to the quasi-instantaneous (SLRderived) center of mass as well as to the quasiinstantaneous scale as realized by VLBI and SLR systems. 3.1. (2) JTRF is non-linear Non-linearity should be here understood in a broader sense: JPL frame products can be best described as nonlinear by virtue of station position process noise based upon random walk. This property is peculiar to JTRF, in that neither ITRF nor DTRF incorporates randomness in their station motion models (cf Table 3).

Table 3 Synoptic outline of the TRF2014 combined products released by the official IERS ITRS Combination centres hosted at LAREG-IGN (ITRF), JPLCaltech (JTRF), DGFI-TUM (DTRF). CM in the Row Origin stands for Center of Mass. Additional pieces of information on the TRS-related activities promoted under the aegis of the IERS can be found on https://www.iers.org/IERS/EN/Organization/ITRSCombinationCentres/ITRSCC.html.

Combination Center Computational Code Frame Type Estimator Process Noise Origin Scale Orientation Station Motion Model a b c d e f (1) (2) (3) (4)

JTRF

ITRF

DTRF

JPL-Caltech KALREF(1) Time Series Kalman Filter Random Walk Instantaneous SLR CM Instantaneous VLBI/SLR RMCd to ITRF2008(4) Trend, Annual, Semi-Annual

LAREG-IGN CATREF(2) Parametric Least-square None Long-term Mean SLR CM VLBI/SLRb RMC to ITRF2008 Trend, Periodic, PSDf

DGFI-TUM DOGS(3) Parametric Least-squarea None Long-term Mean SLR CM VLBI/SLRc NNRe to ITRF2008 Trend, Non-Tidal Corrections

Based on the inversion of accumulated Normal Equations. Simple Average of VLBI and SLR Scales. Weighted Average of VLBI and SLR Scales. RMC stands for Rotational Minimal Constraints (see Eq. (11)). NNR stands for No-Net-Rotation (see Eq. (16)). PSD stands for Post-Seismic-Displacement parametric model. Abbondanza et al. (2017). Altamimi et al. (2007) Seitz et al. (2012). Altamimi et al. (2011).

Please cite this article as: C. Abbondanza, T. M. Chin, R. S. Gross et al., A sequential estimation approach to terrestrial reference frame determination, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.11.016

C. Abbondanza et al. / Advances in Space Research xxx (2019) xxx–xxx

The adoption of station position process noise imparts to JPL frame products a distinctive character. In allowing for deviations from the station motion model, the term wxi;k  N ð0; r2 Þ in Eq. (8) captures non-secular and nonseasonal components of the displacements. This enables, in the absence of technique-specific systematic effects that might affect the smoothed coordinates, a more realistic description of the time-variable nature of the deformable solid Earth. In truth, localized/site-dependent phenomena such as thermoelastic deformations of the bedrock to which geodetic sensors are anchored (see e.g. Dong et al., 2002; Xu et al., 2017), thermal tilting of GNSS monuments (Yan et al., 2009), gravity induced structural deformations of VLBI radio-telescopes (Sarti et al., 2011), SLR stationdependent range biases (Appleby et al., 2016), instabilities in the definition of phase center for DORIS beacons (Willis et al., 2007; Tourain et al., 2016; Saunier et al., 2016) if unaccounted for within the data reductions, are able to perturb, to a different degree and often with distinctive patterns, the SG positioning and might, in turn, affect the quality of random-walk based non-linear frames. Since random walk in JPL products can be controlled by varying the amount of white noise injected during the filtering, two diametrically opposing setups are in principle possible and worth illustrating. On the one hand, the algorithm produces the largest smoothing when no process noise is adopted. When no stochastic perturbations to the state parameters are allowed, the algorithm degenerates, as mentioned in Section 2.2, into a recursive least square estimator and the final combined positions are constrained to follow the dynamical model. On the other, the adoption of large (virtually infinite) values of station position process noise would lead to recover temporally unconstrained positions close to the observed values. In standard practice, JPL frame solutions correspond to intermediate cases wherein the site-dependent variances r2 are determined analysing time series of loading-induced elastic displacements from the three fluid components, i.e. non-tidal atmosphere, oceans and continental hydrology (see e.g. Wu et al., 2015; Abbondanza et al., 2017). The elastic displacements are separately de-trended, de-seasonalized and subsequently aggregated in weekly averages. The residuals from each separate fluid component are added up and first differences of the weekly global displacements are computed. The variances of the differenced time series at each site characterize the site-dependent position process noise wxi;k of Eq. (7). For JTRF2014, median values of the square root of the position process noise variances are found to be 1.20, 0.97, 2.22 mm with associated standard deviations of 0.10, 0.13, 0.87 mm for the North, East and Height component of the displacements, respectively. Although none of the loading models used so far in JPL products accounts for either cryospheric (e.g. ice mass loss) or seismic processes, we found that the loading-derived variances are large enough to allow for an adequate representation of

9

some of the rapidly varying, non-linear and non-seasonal geophysical signals in the position time series. 3.1.1. (3) Is JTRF a regularized frame product? According to the set of conventional definitions reported in Chapter 4 of Petit and Luzum (2010), frame products should rely on the notion of regularized coordinates in their relation to the instantaneous position of points (i.e. geodetic stations) anchored to Earth’s crust: X DX i ðtÞ ð20Þ X ðt Þ ¼ X R ðt Þ þ i

where the decomposition in Eq. (20) dissects the instantaneous –and therefore time-variable– station position X ðtÞ into regularized coordinates X R ðtÞ superposed on highfrequency (mostly) tidal-related displacements DX i ðtÞ. The former are easily determinable from the frame realization by applying a Galilean transformation: X R ðtÞ ¼ X ðt0 Þ þ V  ðt  t0 Þ

ð21Þ

The latter are accessible through conventional geophysical models (cf Chapter 7 of Petit and Luzum, 2010) which are removed at the observation level within the SG reduction procedures. Both ITRF and DTRF provide the users with regularized coordinates linked to a frame whose origin is given by the long-term mean center of mass as observed by SLR. While essentially providing regularized products, JPL departs from the linear representation implicit in Eq. (21) and constructs time series frames based on the additive model of Eq. (8). JPL frame products are to be considered tide-free, regularized frames that, contrary to secular frames, contain loading-induced non-tidal displacements. It is nevertheless important to observe that, insofar that frame inputs include residual tidal effects on the station displacements attributable to body-tide and ocean-tide model deficiencies (see e.g. Section 5 of Kang et al., 2015), JPL frame products, if constructed with a larger amount of process noise, might carry as well such tidal residual signatures. In this particular instance, JPL products may not be thought of as a purely regularized frame. Moreover, JTRF being a sub-secular frame, its coordinates embed the quasi-instantaneous variability due to smoothed SLR-sensed geocenter motion as well as the time-dependent fluctuations of the quasi instantaneous VLBI/SLR scale (cf Fig. 5 in Abbondanza et al., 2017). Users are therefore warned that a larger degree of variability is inherently present in JPL products, when compared to secular products. We note that both ITRF and DTRF are extending their frame products beyond the purely linear representation of the station motions (cf Table 3). In acknowledging the role that non-linear phenomena have on the time-varying deformable Earth and their impact on the frame quality, IGN chose to enhance, for

Please cite this article as: C. Abbondanza, T. M. Chin, R. S. Gross et al., A sequential estimation approach to terrestrial reference frame determination, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.11.016

10

C. Abbondanza et al. / Advances in Space Research xxx (2019) xxx–xxx

ITRF2014 (Altamimi et al., 2016), its stacking model by adjusting periodic and post-seismic displacement (PSD) terms. Of the three IERS-TRF products currently available, ITRF is the only one providing PSD corrections. The introduction of additional parametric functions to describe time-variable non-linear phenomena leads to the formulation of so-called augmented parametric reference frames (Altamimi et al., 2019). As was the case for ITRF2014, Eq. (20) can therefore be expanded as follows: X p X X ðt Þ ¼ X R ðt Þ þ DX j ðtÞ þ DX i ðtÞ ð22Þ j

i

P where j DX pj ðtÞ is a sum of parametric functions of time describing post-seismic displacements of the stations (Altamimi et al., 2019). An opposing strategy for handling non-linearities is utilized by DGFI in DTRF2014 (Seitz et al., 2016). In analogy with standard processing of GRACE data where high-frequency non-tidal mass variability is removed from the range-rate observations and can be optionally restored by the users (Dobslaw et al., 2017), DGFI applies a remove-restore approach to its frame reduction. Rather than fit harmonic terms to position time series, DGFI removes from the SG coordinates geophysical models of non-tidal elastic displacements whose temporal and spatial sampling has been made consistent via spatial interpolation and numerical integration with SG data, and produces after-the-fact piecewise linear estimates of station positions. DTRF users, provided with the removed models, may optionally restore the non-tidal loading variability into the frame product. 3.1.2. (4) JTRF adopts a time-series representation By virtue of the sequential nature of the algorithm, JPL frames are represented through time-dependent parameters at a fixed (weekly) time step, e.g. smoothed time-variable station positions, similarity transformation parameters, data assimilation residuals. JPL time-dependent parameters can be gathered and output in time series of SINEX files. The time-series representation characterizes and distinguishes JTRF from ITRF and DTRF, both adopting a parameterized representation. ITRF and DTRF products are typically released through single monolithic SINEX files containing geocentric station positions expressed at a reference epoch, velocities, and EOPs along with their full covariance. The ultimate goal of the JPL frame products is to unify all of the SG inputs into a self-consistent frame. Bridged together in a unified TRF, the station position and Earth rotation time series form longer and more robust SG records which can be valuably used to infer highly accurate time-variable signals of the Earth deformation, its rotation, and the geocenter motion. An example of the station position time series outputted by the algorithm when analysing ITRF2014 frame inputs is given in Fig. 1 which illustrates the co-located GNSS and VLBI stations at Tsukuba (Japan). Table 4 recalls the main

features of the SG inputs adopted to determine JTRF2014 and the time series of Fig. 1 and lists recently published references wherein a thorough discussion on the nature of the solutions can be found. The black dots in the plots represent the observed positions of the two stations whereas the red solid lines indicate the estimated positions derived from JTRF2014. The coseismic displacement of the two observing stations during the 9.1 Mw T ohoku megaquake of March 11, 2011 (Koketsu et al., 2011; Thorne, 2018) is readily apparent, particularly in the East (see top panels of Fig. 1) and Height (bottom panels) components. By resetting the variance of the station position process noise to 25.0 m2 at the epoch of the seismic event, the algorithm is able to correctly represent the sizable co-seismic dislocation. At the same time, a discontinuity in the velocity before and after the quake is introduced both for VLBI and GNSS. This is, likewise, achieved by activating the velocity process noise and by instantaneously setting its variance to the 2 value of 25.0 ðm=sÞ . After the quake, the velocity process noise is deactivated. In so doing, the sequential algorithm is forced to estimate a new (constant) value of velocity after the seismic event. Breaking the velocity continuity in this fashion is equivalent to producing piecewise constant estimates of the station velocities. In between discontinuities, the nominal value of the process noise at this site causes the optimal station positions to be a smoothed version of the observed position. The dynamical model of Eq. (7) is used to interpolate across gaps in the observations and to extrapolate the station position to epochs before the station itself started observing (before 1994.0) or after it stopped observing (not shown in the plot). Since the VLBI station (7345) at Tsukuba is co-located with a GNSS station that started observing before it and because of the co-motion constraints applied to co-located stations, the estimated position of the GNSS station is transferred to the VLBI station at the epochs in which 7345 was not yet operating. This is most clearly seen in the Height component (right bottom panel of Fig. 1) in the time window 1994-1998. It has to be observed that, in the case of Tsukuba, comotion constraints between GNSS and VLBI have been applied both before and after the seismic event. 4. Perspectives JPL products are sub-secular frames optimally estimated by adopting a Kalman filter and RTS smoother. JPL frames offer to the geodetic community interesting features. With its sub-secular time scale, JTRF gives access to the quasi-instantaneous center of mass as sensed by SLR as well as to a quasi-instantaneous scale as realized by VLBI and SLR which, if accurately realized by data, can be used by geodetic systems with inaccurate scale information for monitoring time-variable station coordinates. If adopted by the SG community, sub-secular frame products allows the origin duality between secular center

Please cite this article as: C. Abbondanza, T. M. Chin, R. S. Gross et al., A sequential estimation approach to terrestrial reference frame determination, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.11.016

C. Abbondanza et al. / Advances in Space Research xxx (2019) xxx–xxx

11

Table 4 SG inputs adopted for JTRF2014. T in the first column designates the SG techniques (P identifies GNSS, R VLBI, L SLR whereas D is DORIS). TC in the second column identifies the Technique Center. TS in the third column designates the time span for each technique. TR in the fourth column relates to the temporal resolution, ST in the fifth column designates the number of stations adopted in JTRF2014 (Parenthesized is the total number of stations provided in the unedited data by the Technique Centers). The last column, Source, provides references. T

TC

P R L D

IGS IVS ILRS IDS

TS

TR

ST

Source

1994.0–2015.1 1979.5–2015.0 1983.0–2015.0 1993.0–2015.0

Daily Session-Wise Fortnightly (83–93)/ Weekly Weekly

671 (1845) 71 (158) 71 (138) 159 (138)

Rebischung et al. (2016) Bachmann et al. (2016) Luceri and Pavlis (2016), Pearlman et al. (2019) Moreaux et al. (2016)

Fig. 1. Position Time series for the co-located GNSS and VLBI stations at Tsukuba (Japan). The plots in the left panel relate to the GNSS station, whereas those on the right side to VLBI. Black dotted lines represent the station position observations assimilated by sequential estimator, the red solid lines are the weekly smoothed coordinates obtained as a result of the combination. Solid green vertical lines mark the time tags of station position discontinuities. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

of mass and sub-secular center of figure/network affecting GNSS solutions Helmert-transformed into ITRF/DTRF to be overcome. Such an aspect is critically important to those observing systems, such as radar altimeters, aimed at quantifying global mean sea level (Melachroinos et al., 2013; Couhert et al., 2015). As efforts to improve JPL frame products continue, positive feedback on JTRF2014 has been received. When adopted in precise orbit determination of high and low Earth orbiters JTRF2014-derived time series are reported to make valuable contributions in terms of an improved geocenter and reduction of post-fit residuals (Zelensky et al., 2018; Rudenko et al., 2018).

State predictions beyond the data assimilation time span and based on Eq. (3), albeit available, were not included in JTRF2014 release. Predictions are all of the more crucial in time-series products, in that they provide the only way to access the frame beyond its assimilation time span. JPL frame predictions are considered to be, at the current stage of development, experimental products still exhibiting deficiencies, particularly in relation to those sites affected by post-seismic relaxation. Unlike ITRF, JTRF does not implement yet PSD terms in the state-space model of Eq. (7). As a result, predictions for seismic sites are largely inaccurate and not adoptable in current space-geodetic reduction procedures.

Please cite this article as: C. Abbondanza, T. M. Chin, R. S. Gross et al., A sequential estimation approach to terrestrial reference frame determination, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.11.016

12

C. Abbondanza et al. / Advances in Space Research xxx (2019) xxx–xxx

Offsets of the Helmert transformation parameters relating JPL products to ITRFs are all well within 1 mm-level, with the exception of the scale whose larger offset is attributable to the different strategies adopted in averaging the VLBI/SLR information. Drifts of the Helmert transformation parameters are found to be smaller than 0.18 mm/yr, thus proving the secular parts of the two frames are highly consistent. (see e.g. Wu et al. (2015) for JTRF to ITRF2008 and Abbondanza et al. (2017) for JTRF to ITRF2014 comparisons). In this instance, the requirements formulated in Gross et al. (2009) demand TRF realizations be accurate to the level of 1 mm and stable to within the level of 0.1 mm/yr. Although similarity transformation parameters between TRFs do not provide a metric valid for assessing TRF accuracy (i.e. degree of closeness to some extrinsic physical truth), they establish the level of agreement between JTRF and ITRF products to be acceptably compliant with the GGOS specifications. Frame accuracy is an extremely arduous notion to pinpoint and quantify in that it entails evaluating and comparing TRF products to exogeneous data types, be they geophysical models or geodetic observations not directly adopted in frame reductions (see e.g. Collilieux et al., 2014). Wu et al. (2011) inverted ITRF2008 station velocities (Altamimi et al., 2011), linear trends computed from space gravimetry data and an ocean bottom pressure (OBP) model and found that the ITRF2008 origin stability is in the order of 0:5 mm/yr. Since such metric describes the level of consistency between the geocentric velocity fields as observed by SG and those implied by degree-1 surface mass variation and load-induced deformation extracted from trends in relative deformation, OBP, and space gravity, it can be used as a proxy for TRF accuracy. It is worth observing that, because the origin drift between ITRF2014 and ITRF2008 is practically negligible (Altamimi et al., 2016), the results described by Wu et al. (2011) are directly applicable to and hold true for ITRF2014. The level of accuracy for current frame products, then, turns out to be five times larger than the value of 0:1 mm/yr set by Gross et al. (2009) and Blewitt et al. (2010). Furthermore, TRF sequential approaches provide an ideal framework wherein the frame can be recursively updated as new observations become available. ITRF realizations are currently determined every 3 to 5 years. Such a long interval between ITRF releases makes the frame products vulnerable to temporal degradation and limits their ability to accurately extrapolate station positions and Earth Orientation Parameters into the future. In a TRF sequential algorithm, because the pair ðX k ; P k Þ can be saved at the last step of the full blown frame reduction, as additional measurements become available the filter can be restarted from the saved state and run forward in time to assimilate the additional measurements. The updates, produced with no need to regenerate the full solution, will therefore improve the predictions through the constraints offered by the newer observations. This is potentially appealing for all those scientific applications in which

frame predictability is crucial, such as sea level change, vertical land motion determination, geodynamics and precise orbit determination. JPL is actively researching the feasibility of producing timely updates to TRFs and issues related to the stability, accuracy and precision of updated products are currently being evaluated. With the ever-increasing amount of global spacegeodetic data sets (Blewitt et al., 2018), processing time for JPL products might become a limiting factor. For JTRF2014, limited by the high computational burden our sequential estimator entails, we opted for reducing the number of stations adopted to 972 (cf Table 4), against the 1499 stations included in ITRF2014. In an effort to modernize the computational code, the TRF team at JPL has selected lines of inquiries which are being actively researched as part of JTRF development/ improvement plan:  Algorithm enhancement.The computational code will be upgraded to a time-variable sequential algorithm based on the assimilation of the information pair ðX k ; Kk Þ (where Kk ¼ P 1 is the information matrix, (see e.g. k Chin, 2001)) and the adoption of a square-root formalism (see e.g. Bierman, 2006; Simon, 2006). By transitioning to a square-root information filter, the new algorithm will take full advantage of (i) the increased numerical stability achieved as a result of the squareroot representation, and (ii) the flexibility ensured by the propagation of the information pair during the initialization. The information filter can easily handle those instances when the initial uncertainty on the state parameters is undetermined (i.e. P0 ¼ 1) by setting K0 ¼ 0.  State vector augmentation.The state representation will be expanded to account for extra-parameters such as post-seismic displacements and radio-source coordinates that can be assimilated along with VLBI observations of precession and nutation to establish a direct link to the celestial reference frame (Kwak et al., 2018). Also, to allow for a more sophisticated representation of the stochastic state parameters (via e.g. high-order autoregressive processes), the state vector will be augmented by including additional noise process parameters.  Dynamical model enhancement.The state-space model of Eq. (6) will be complemented with EOP dynamics and a more appropriate formulation of their underlying stochastic processes as documented in Chin et al. (2009). Potential areas of further improvement also relate to e.g. strategies for providing accurate predictions of the station motions, data-editing pipelines, break detection, etc. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Please cite this article as: C. Abbondanza, T. M. Chin, R. S. Gross et al., A sequential estimation approach to terrestrial reference frame determination, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.11.016

C. Abbondanza et al. / Advances in Space Research xxx (2019) xxx–xxx

Acknowledgements This work was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (NASA). We thankfully acknowledge funding support from NASA’s Space Geodesy Program. Comments and feedback from the Editor-in-Chief, Pascal Willis, helped improving the presentation of the material. Also, we are highly appreciative of the thourough assessment of the manuscript from two anonymous reviewers. JTRF data (input SINEX files) and products (time series of station positions, Earth Orientation Parameters, Geocenter Motion, data assimilation residuals) are available upon request to [email protected]. B Soja’s research was supported by an appointment to the NASA Postdoctoral Program at the NASA Jet Propulsion Laboratory, administered by Universities Space Research Association under contract with NASA. Maps in this manuscript have been produced with Carto Py, freely available under the terms of the GNU License. Ó 2019 California Institute of Technology. All rights reserved. U.S. Government sponsorship acknowledged. References Abbondanza, C., Chin, T.M., Gross, R.S., Heflin, M.B., Parker, J.W., Soja, B.S., van Dam, T., Wu, X., 2017. JTRF2014, the JPL Kalman filter and smoother realization of the International Terrestrial Reference System. J. Geophys. Res.-Sol. Ea. 122, 8474–8510. https://doi. org/10.1002/2017JB014360. Altamimi, Z., Boucher, C., Gambis, D., 2005. Long-term stability of the terrestrial reference frame. Adv. Space Res. 36, 342–349. https://doi. org/10.1016/j.asr.2005.03.068, satellite Dynamics in the Era of Interdisciplinary Space Geodesy. Altamimi, Z., Boucher, C., Willis, P., 2005. Terrestrial reference frame requirements within GGOS perspective. J. Geodyn. 40, 363–374. https://doi.org/10.1016/j.jog.2005.06.002. 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.-Sol. Ea. 112, B09401. https://doi.org/ 10.1029/2007JB004949. Altamimi, Z., Collilieux, X., Me´tivier, L., 2011. ITRF2008: an improved solution of the International Terrestrial Reference Frame. J. Geodesy 85, 457–473. https://doi.org/10.1007/s00190-011-0444-4. Altamimi, Z., Rebischung, P., Collilieux, X., Me´tivier, L., Chanard, K., 2019. Review of Reference Frame Representations for a Deformable Earth. International Association of Geodesy Symposia. Springer, Berlin, Heidelberg. https://doi.org/10.1007/1345_2019_66. Altamimi, Z., Rebischung, P., Me´tivier, L., Collilieux, X., 2016. ITRF2014: a new release of the International Terrestrial Reference Frame modeling nonlinear station motions. J. Geophys. Res.-Sol. Ea. 121, 6109–6131. https://doi.org/10.1002/2016JB013098. Altamimi, Z., Sillard, P., Boucher, C., 2002. ITRF2000: a new release of the International Terrestrial Reference Frame for Earth science applications. J. Geophys. Res.-Sol. Ea. 107 (B10), 2114–2133.

13

Altamimi, Z., Sillard, P., Boucher, C., 2003. The impact of a No-NetRotation Condition on ITRF2000. Geophys. Res. Lett. 30, 36-1–36-4. https://doi.org/10.1029/2002GL016279. Appleby, G., Rodrı´guez, J., Altamimi, Z., 2016. Assessment of the accuracy of global geodetic satellite laser ranging observations and estimated impact on ITRF scale: estimation of systematic errors in LAGEOS observations 1993–2014. J. Geodesy 90, 1371–1388. https:// doi.org/10.1007/s00190-016-0929-2. Bachmann, S., Thaller, D., Roggenbuck, O., Lo¨sler, M., Messerschmitt, L., 2016. IVS contribution to ITRF2014. J. Geodesy 90, 631–654. https://doi.org/10.1007/s00190-016-0899-4. Bierman, G.J., 2006. Factorization Methods for Discrete Sequential Estimation. Dover Books on Mathematics Series. Dover Publications. (unabridged republication of the original volume published in 1977 by Academic Press, Inc., N.Y.). Blewitt, G., 2015. Terrestrial Reference Frame Requirements for Studies of Geodynamics and Climate Change. In: Springer International Publishing Switzerland (Ed.), Reference Frames For Applications in Geosciences and Technology. doi: https://doi.org/10.1007/1345_2015_ 142. Blewitt, G., Hammond, W.C., Kreemer, C., 2018. Harnessing the GPS data explosion for interdisciplinary science. EOS 99. https://doi.org/ 10.1029/2018EO104623. Chin, T.M., 2001. On Kalman filter solution of space-time interpolation. IEEE Trans. Image Proc. 10, 663–666. https://doi.org/10.1109/ 83.913601. Blewitt, G., Altamimi, Z., Davis, J., Gross, R., Kuo, C.Y., Lemoine, F.G., Moore, A.W., Neilan, R.E., Plag, H.P., Rothacher, M., Shum, C.K., Sideris, M.G., Scho¨ne, T., Tregoning, P., Zerbini, S., 2010. Geodetic Observations and Global Reference Frame Contributions to Understanding Sea-Level Rise and Variability. Wiley-Blackwell, pp. 256–284. https://doi.org/10.1002/9781444323276.ch9. Chin, T.M., Gross, R.S., Boggs, D.H., Ratcliff, J.T., 2009. Dynamical and Observation Models in the Kalman Earth Orientation Filter. Technical Report DSMS No. 821–104 Rev B, JPL D-17235. Jet Propulsion Laboratory. Chin, T.M., Gross, R.S., Dickey, J.O., 2004. Modeling and forecast of the polar motion excitation functions for short-term polar motion prediction. J. Geodesy 78, 343–353. https://doi.org/10.1007/s00190004-0411-4. Collilieux, X., Altamimi, Z., Argus, D.F., Boucher, C., Dermanis, A., Haines, B.J., Herring, T.A., Kreemer, C.W., Lemoine, F.G., Ma, C., MacMillan, D.S., Ma¨kinen, J., Me´tivier, L., Ries, J., Teferle, F.N., Wu, X., 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. Springer, Berlin, Heidelberg, pp. 197–202. Collilieux, X., Altamimi, Z., Coulot, D., van Dam, T., Ray, J., 2010. Impact of loading effects on determination of the International Terrestrial Reference Frame. Adv. Space Res. 45, 144–154. https:// doi.org/10.1016/j.asr.2009.08.024. Collilieux, X., Altamimi, Z., Coulot, D., Ray, J., Sillard, P., 2007. Comparison of very long baseline interferometry, GPS and satellite laser ranging height residuals from ITRF2005 using spectral and correlation methods. J. Geophys. Res.-Sol. Ea. 112, B12403. https:// doi.org/10.1029/2007JB004933. Couhert, A., Cerri, L., Legeais, J.F., Ablain, M., Zelensky, N.P., Haines, B.J., Lemoine, F.G., Bertiger, W.I., Desai, S.D., Otten, M., 2015. Towards the 1mm/y stability of the radial orbit error at regional scales. Adv. Space Res. 55, 2–23. https://doi.org/10.1016/j.asr.2014.06.041. Davis, J.L., Wernicke, B.P., Tamisiea, M.E., 2012. On seasonal signals in geodetic time series. J. Geophys. Res.-Sol. Ea. 117. https://doi.org/ 10.1029/2011JB008690. Dehant, V., Mathews, P.M., 2015. Earth Rotation Variations. In: Schubert, G. (Ed.), Treatise on Geophysics. second ed. Elsevier, Oxford, pp. 263 – 305. doi: https://doi.org/10.1016/B978-0-444-538024.00066-X.

Please cite this article as: C. Abbondanza, T. M. Chin, R. S. Gross et al., A sequential estimation approach to terrestrial reference frame determination, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.11.016

14

C. Abbondanza et al. / Advances in Space Research xxx (2019) xxx–xxx

Dobslaw, H., Bergmann-Wolf, I., Dill, R., Poropat, L., Thomas, M., Dahle, C., Esselborn, S., Koenig, R., Flechtner, F., 2017. A new highresolution model of non-tidal atmosphere and ocean mass variability for de-aliasing of satellite gravity observations: AOD1B RL06. Geophys. J. Int. 211, 263–269. https://doi.org/10.1093/gji/ggx302. Dong, D., Fang, P., Bock, Y., Cheng, M.K., Miyazaki, S., 2002. Anatomy of apparent seasonal variations from GPS-derived site position time series. J. Geophys. Res.-Sol. Ea 107. https://doi.org/10.1029/ 2001JB000573, ETG 9-1–ETG 9-16. Dong, D., Herring, T.A., King, R.W., 1998. Estimating regional deformation from a combination of space and terrestrial geodetic data. J. Geodesy 72, 200–214. https://doi.org/10.1007/s001900050161. Fukumori, I., Malanotte-Rizzoli, P., 1995. An approximate Kalman filter for ocean data assimilation: an example with an idealized Gulf Stream model. J. Geophys. Res.-Oceans 100, 6777–6793. https://doi.org/ 10.1029/94JC03084. Gelb, A., 1974. Applied Optimal Estimation. MIT Press https://books.google.com/books?id=KlFrn8lpPP0C. Gross, R., Beutler, G., Plag, H.P., 2009. Integrated scientific and societal user requirements and functional specifications for the GGOS. In: Plag, H.P., Pearlman, M. (Eds.), Global Geodetic Observing System, Meeting the Requirements of a Global Society on a Changing Planet in 2020. Springer-Verlag, Berlin Heidelberg. pp. 209–223. doi: https://doi. org/10.1007/978-3-642-02687-4 (chapter 7). Gross, R.S., 2015. Earth Rotation Variations - Long Period. In: Schubert, G. (Ed.), Treatise on Geophysics, second ed. Elsevier, Oxford, pp. 215–261. doi: https://doi.org/10.1016/B978-0-444-53802-4.00059-2. Gross, R.S., Eubanks, T.M., Steppe, J.A., Freedman, A.P., Dickey, J.O., Runge, T.F., 1998. A Kalman-filter-based approach to combining independent Earth-orientation series. J. Geodesy 72, 215–235. https:// doi.org/10.1007/s001900050162. Hay, C.C., Morrow, E., Kopp, R.E., Mitrovica, J.X., 2013. Estimating the sources of global sea level rise with data assimilation techniques. Proc. Natl. Acad. Sci. USA 110, 3692–3699. https://doi.org/10.1073/ pnas.1117683109. Herring, T.A., Davis, J.L., Shapiro, I.I., 1990. Geodesy by radio interferometry: the application of Kalman filtering to the analysis of very long baseline interferometry data. J. Geophys. Res.-Sol. Ea. 95, 12561–12581. https://doi.org/10.1029/JB095iB08p12561. Ka´lma´n, R.E., 1960. A new approach to linear filtering and prediction problems. J. Basic Eng. 82, 35–45. https://doi.org/10.1115/1.3662552. Kang, K., Wahr, J., Heflin, M., Desai, S., 2015. Stacking global GPS verticals and horizontals to solve for the fortnightly and monthly body tides: Implications for mantle anelasticity. J. Geophys. Res.-Sol. Ea. 120, 1787–1803. https://doi.org/10.1002/2014JB011572. Koketsu, K., Yokota, Y., Nishimura, N., Yagi, Y., Miyazaki, S., Satake, K., Fujii, Y., Miyake, H., Sakai, S., Yamanaka, Y., Okada, T., 2011. A unified source model for the 2011 Tohoku earthquake. Earth Planet Sci. Lett. 310, 480–487. https://doi.org/10.1016/j.epsl.2011.09.009. Kwak, Y., Bloßfeld, M., Schmid, R., Angermann, D., Gerstl, M., Seitz, M., 2018. Consistent realization of celestial and terrestrial reference frames. J. Geod. 92, 1047–1061. https://doi.org/10.1007/s00190-0181130-6. Legrand, J., 2007. Champ de vitesses de l’ITRF. Proprie´te´s cine´matiques de la crouˆte terrestre et Condition de non rotation globale. Ph.D. thesis. Observatoire de Paris. URL https://hal.archives-ouvertes.fr/tel02071422. Lichten, S.M., Border, J.S., 1987. Strategies for high-precision Global Positioning System orbit determination. J. Geophys. Res.-Sol. Ea. 92, 12751–12762. https://doi.org/10.1029/JB092iB12p12751. Luceri, V., Pavlis, E., 2016. The ILRS contribution to ITRF2014. Int Laser Rang Serv, Greenbelt, MD, http://itrf.ign.fr/ITRF_solutions/ 2014/doc/ILRS-ITRF2014-description.pdf (on-line Accessed 4 August 2017). Melachroinos, S.A., Lemoine, F.G., Zelensky, N.P., Rowlands, D.D., Luthcke, S.B., Bordyugov, O., 2013. The effect of geocenter motion on Jason-2 orbits and the mean sea level. Adv. Space Res. 51, 1323–1334. https://doi.org/10.1016/j.asr.2012.06.004.

Moreaux, G., Lemoine, F.G., Capdeville, H., Kuzin, S., Otten, M., Sˇteˇfa´nek, P., Willis, P., Ferrage, P., 2016. The International DORIS Service contribution to the 2014 realization of the International Terrestrial Reference Frame. Adv. Space Res. 58, 2479–2504. https:// doi.org/10.1016/j.asr.2015.12.021. Moritz, H., Mu¨ller, I.I., 1987. Earth Rotation: Theory and Observation. Ungar New York. Papoulis, A., Pillai, S., 2001. Probability, Random Variables and Stochastic Processes, Fourth ed. McGraw-Hill Education https://books.google.com/books?id=ZIAMnwEACAAJ. Pearlman, M., Arnold, D., Davis, M., Barlier, F., Biancale, R., Vasiliev, V., Ciufolini, I., Paolozzi, A., Pavlis, E.C., Sos´nica, K., Bloßfeld, M., 2019. Laser geodetic satellites: a high-accuracy scientific tool. J. Geodesy. https://doi.org/10.1007/s00190-019-01228-y. Petit, G., Luzum, B., 2010. IERS Conventions (2010). International Earth Rotation and Reference Systems Service Central Bureau Bundesamt fu¨r Kartographie und Geoda¨sie Richard-Strauss-Allee 11 60598 Frankfurt am Main Germany. IERS Technical Note, No. 36. Ray, J., Altamimi, Z., Collilieux, X., van Dam, T., 2008. Anomalous harmonics in the spectra of GPS position estimates. GPS Solut. 12, 55– 64. https://doi.org/10.1007/s10291-007-0067-7. Ray, J.R., 1996. Measurements of length of day using the global positioning system. J. Geophys. Res.-Sol. Ea. 101, 20141–20149. https://doi.org/10.1029/96JB01889. Rebischung, P., Altamimi, Z., Ray, J., Garayt, B., 2016. The IGS contribution to ITRF2014. J. Geodesy 90, 611–630. https://doi.org/ 10.1007/s00190-016-0897-6. Rudenko, S., Bloßfeld, M., Mu¨ller, H., Dettmering, D., Angermann, D., Seitz, M., 2018. Evaluation of DFTR2014, ITRF2014 and JTRF2014 by Precise Orbit Determination of SLR Satellites. IEEE T. Geosci. Remote 56, 3148–3158. https://doi.org/10.1109/ TGRS.2018.2793358. Sarti, P., Abbondanza, C., Petrov, L., Negusini, M., 2011. Height bias and scale effect induced by antenna gravitational deformations in geodetic VLBI data analysis. J. Geodesy 85, 1–8. https://doi.org/10.1007/ s00190-010-0410-6. Saunier, J., Auriol, A., Tourain, C., 2016. Initiating an error budget of the DORIS ground antenna position: genesis of the Starec antenna type C. Adv. Space Res. 58, 2717–2724. https://doi.org/10.1016/j. asr.2016.02.013. Seitz, M., Angermann, D., Bloßfeld, M., Drewes, H., Gerstl, M., 2012. The 2008 DGFI realization of the ITRS: DTRF2008. J. Geodesy 86, 1097–1123. https://doi.org/10.1007/s00190-012-0567-2. Seitz, M., Bloßfeld, M., Angermann, D., Schmid, R., Gerstl, M., Seitz, F., 2016. The new DGFI-TUM realization of the ITRS: DTRF2014 (data). doi: https://doi.org/10.1594/PANGAEA.864046. Sillard, P., Boucher, C., 2001. A review of algebraic constraints in terrestrial reference frame datum definition. J. Geodesy 75, 63–73. https://doi.org/10.1007/s001900100166. Simon, D., 2006. Optimal State Estimation: Kalman, H Infinity, and Nonlinear Approaches. Wiley https://books.google.com/books?id= UiMVoP_7TZkC. Stengel, R.F., 1994. Optimal Control and Estimation. Dover Books on Mathematics. Dover Publications, Inc., New York. Tapley, B.D., Bettadpur, S., Watkins, M., Reigber, C., 2004. The Gravity Recovery And Climate Experiment: mission overview and early results. Geophys Res Lett 31. https://doi.org/10.1029/2004GL019920. Thorne, L., 2018. A review of the rupture characteristics of the 2011 Tohoku-oki Mw 9.1 earthquake. Tectonophysics 733, 4–36. https:// doi.org/10.1016/j.tecto.2017.09.022. Tourain, C., Moreaux, G., Auriol, A., Saunier, J., 2016. DORIS Starec ground antenna characterization and impact on positioning. Adv. Space Res. 58, 2707–2716. https://doi.org/10.1016/j.asr.2016.05.013, Scientific Applications of DORIS in Space Geodesy. Willis, P., Haines, B., Kuang, D., 2007. DORIS satellite phase center determination and consequences on the derived scale of the Terrestrial Reference Frame. Adv. Space Res. 39, 1589–1596. https://doi.org/ 10.1016/j.asr.2007.01.007.

Please cite this article as: C. Abbondanza, T. M. Chin, R. S. Gross et al., A sequential estimation approach to terrestrial reference frame determination, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.11.016

C. Abbondanza et al. / Advances in Space Research xxx (2019) xxx–xxx Wu, X., Abbondanza, C., Altamimi, Z., Chin, T.M., Collilieux, X., Gross, R.S., Heflin, M.B., Jiang, Y., Parker, J.W., 2015. KALREF - a Kalman filter and time series approach to the International Terrestrial Reference Frame realization. J. Geophys. Res.-Sol. Ea. 120, 3775– 3802. https://doi.org/10.1002/2014JB011622. Wu, X., Collilieux, X., Altamimi, Z., Vermeersen, B.L.A., Gross, R.S., Fukumori, I., 2011. Accuracy of the International Terrestrial Reference Frame origin and earth expansion. Geophys. Res. Lett. 38. https://doi.org/10.1029/2011GL047450. l13304. Xu, X., Dong, D., Fang, M., Zhou, Y., Wei, N., Zhou, F., 2017. Contributions of thermoelastic deformation to seasonal variations in GPS station position. GPS Solut. 21, 1265–1274. https://doi.org/ 10.1007/s10291-017-0609-6.

15

Yan, H., Chen, W., Zhu, Y., Zhang, W., Zhong, M., 2009. Contributions of thermal expansion of monuments and nearby bedrock to observed GPS height changes. Geophys. Res. Lett. 36. https://doi.org/10.1029/ 2009GL038152. Zelensky, N.P., Lemoine, F.G., Beckley, D.B., Chinn, D.S., Pavlis, D., 2018. Impact of ITRS 2014 realizations on altimeter satellite precise orbit determination. Adv. Space Res. 61, 45–73. https://doi.org/ 10.1016/j.asr.2017.07.044. Zumberge, J.F., Heflin, M.B., Jefferson, D.C., Watkins, M.M., Webb, F. H., 1997. Precise point positioning for the efficient and robust analysis of GPS data from large networks. J. Geophys. Res.-Sol. Ea. 102, 5005– 5017. https://doi.org/10.1029/96JB03860.

Please cite this article as: C. Abbondanza, T. M. Chin, R. S. Gross et al., A sequential estimation approach to terrestrial reference frame determination, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.11.016