C E Gadda

Abstract

The accurate and precise implementation of the International Terrestrial Reference Frame (ITRF) is a fundamental requirement for the development of Earth System Sciences. The actual realization of the reference frame, in fact, directly impacts a number of different tasks ranging from precise satellite orbit determination to altimeter calibration, satellite antenna offset assessment for Global Navigation Satellite System (GNSS) and validation of center of mass corrections for spacecrafts carrying on board retro-reflectors for Satellite Laser Ranging (SLR). As a consequence, all the studies investigating motions of the Earth's surface, including oceans and ice-sheets, strictly depend on the availability of a reliable TRF that is fundamental for geo-referencing the relevant measurements.

ITRF realizations must then be periodically updated, in order to account for newly acquired observations and for upgrades in data analysis procedures and/or combination methods. Any innovative computation strategy should ameliorate the realization of the frame physical parameters, namely the origin and the scale, upon which a number of scientific applications critically rely. This work addresses the potential of combining GNSS and SLR observations via their co-location on board GPS/GLONASS satellites. GNSS vehicles equipped with retro-reflector arrays can be tracked by SLR ground stations, which allows determining the spacecraft orbits by means of both optical and microwave signals.

In principle, the inter-technique connection so achieved could be exploited for the computation of the ITRF in place of terrestrial ties. These lasts are known to be currently a limiting factor of the frame accuracy because of their inhomogeneous distribution and of their discrepancies with space geodesy estimates due to technique systematic errors. In this study, the strength of the alternative link in orbit has been thoroughly investigated in order to evaluate the performances of the selected space tie approach under the available operational conditions. The analysis focuses on the characterization of the precision, the accuracy and the pertinence of the combined frame parameters. Keywords: ITRF, GNSS, SLR, space ties, terrestrial ties, origin, geocenter, scale Sintesi L'implementazione accurata e precisa del quadro di riferimento terrestre internazionale (International Terrestrial Reference Frame, ITRF) è un requisito fondamentale per lo sviluppo delle Scienze del Sistema Terra. Dalla realizzazione di tale quadro dipendono infatti la precisa determinazione delle orbite satellitari, la calibrazione degli altimetri e la determinazione della posizione dei dispositivi di localiz-zazione a bordo di un satellite rispetto al centro di massa dello stesso. Solo una conoscenza affidabile di queste quantità permette una georeferenziazione accurata delle osservazioni su cui si basano gli studi di deformazione crostale o di variazioni del livello del mare e dello spessore dei ghiacci.

La realizzazione dell'ITRF deve pertanto essere periodicamente aggiornata in modo tale da assimilare le nuove osservazioni acquisite ed i miglioramenti nelle tecniche di analisi dati e/o nelle procedure di combinazione. Lo sviluppo di nuovi metodi per il calcolo dell'ITRF deve avere come scopo il perfezionamento della materializzazione dei parametri fisici del quadro, ovvero l'origine e la scala, da cui dipendono criticamente applicazioni scientifiche, civili ed economiche.

Questo lavoro si pone come obiettivo la valutazione dell'efficacia di una strategia di combinazione basata sulla co-localizzazione in orbita delle tecniche coinvolte. In particolare, è stato investigato il potenziale della connessione tra i sistemi Global Navigation Satellite System (GNSS) e Satellite Laser Ranging (SLR) a bordo dei satelliti GNSS. Grazie ad appositi pannelli di specchi retro-riflettori, é infatti possibile tracciare otticamente l'intera costellazione GLONASS e due veicoli GPS. L'orbita di questi satelliti può quindi essere determinata dall'insieme delle osservazioni delle due tecniche, creando una connessione che, in linea di principio, potrebbe rimpiazzare i "vincoli terrestri" attualmente necessari per il calcolo dell'ITRF. Prima di affidarsi ad un approccio alternativo, è tuttavia necessario valutarne l'efficacia nelle reali condizioni operative. In questo studio, i vincoli selezionati sono stati valutati caratterizzando la precisione, l'accuratezza e la stabilità del quadro di riferimento che ne può essere derivato.

Parole chiave : ITRF, GNSS, SLR, vincoli in orbita, vincoli terrestri, origine, geocentro, scala

Resumé

La mise en oeuvre exacte et précise du repère international de référence terrestre (ITRF) est une exigence fondamentale pour le développement des Sciences du Système Terre. La réalisation du référentiel mondial, en fait, concerne directement de nombreux domaines allant de la détermination précise des orbites des satellites, à la calibration des altimètres, à l'évaluation des étalonnages absolus d'antennes satellites pour le Global Navigation Satellite System (GNSS) et la validation des corrections du vecteur du centre de masse pour les véhicules spatiaux portant à bord des rétro-réflecteurs pour la technique de télémétrie laser sur satellite (SLR). En conséquence, toutes les études portant sur les mouvements de la surface de la Terre, y compris les océans et les calottes glaciaires, dépendent étroitement de la disponibilité d'un repère de référence fiable qui est fondamental pour référencer les mesures pertinentes.

La réalisation de l'ITRF doit alors être périodiquement mise à jour, afin d'intégrer des nouvelles observations et progrès dans les procédures d'analyse des données et/ou des méthodes de combinaison. Toutes les nouvelles stratégies de calcul doivent viser l'amélioration de la réalisation des paramètres physiques du repère, à savoir l'origine et l'échelle, sur lesquels se fondent de façon critique un grand nombre d'études scientifiques et d'applications civiles.

Ce travail se concentre sur le potentiel de combiner les observations GNSS et SLR par leur liens à bord de satellites GPS / GLONASS. En fait, les satellites GNSS équipés de rétro-réflecteurs peuvent être observés par les stations SLR, ce qui permet de déterminer les orbites des satellites à travers les deux signaux : optiques et à micro-ondes. En principe, la connexion inter-technique si réalisée pourrait être exploitée pour le calcul de l'ITRF en place des liens terrestres actuellement utilisés. Ces derniers sont connus pour être aujourd'hui un facteur limitant de la précision du repère en raison de leur distribution inhomogène et de leurs divergences avec les estimations de la géodésie spatiale en conséquence des erreurs systématiques dans les observations. Dans cette étude, la force du lien alternatif en orbite a été soigneusement analysée afin d'évaluer les performances de l'approche de combinaison sélectionnée dans les conditions opérationnelles disponibles. L'investigation porte sur la caractérisation de la précision, de la fiabilité et de la pertinence des paramètres combinés du repère de référence.

Mots-clés : ITRF, GNSS, SLR, liens spatiaux,liens terrestres, origine, geocenter, échelle

Sintesi sostanziale

Lo scopo della geodesia si concretizza in ultima analisi nella caratterizzazione della geometria, del campo di gravità e della rotazione della Terra e delle loro variazioni temporali. Gli studi geodetici riguardano l'osservazione e la modellizzazione dei movimenti crostali, della ridistribuzione di massa tra le diverse geosfere e degli scambi di momento angolare al di sotto e al di sopra della superficie. L'evoluzione di questi tre campi è strettamente interconnessa e richiede la realizzazione di un'unica struttura che permetta di studiare il Sistema Terra nel suo insieme. Una tale cornice è fornita dal sistema di riferimento terrestre (Terrestrial Reference System -TRS) e dalla sua materializzazione fisica, il quadro di riferimento terrestre (Terrestrial Reference Frame -TRF), che permettono di posizionare organicamente gli oggetti situati sulla superficie della Terra o in orbita attorno ad essa. Le implicazioni scientifiche, sociali ed economiche della disponibilità di un TRF preciso e stabile sono state riconosciute dall'Assemblea Generale delle Nazioni Unite nella sua prima risoluzione geospaziale, adottata il 26 febbraio 2015.

La determinazione delle coordinate dell'ITRF, ovvero del particolare quadro adottato dall'International Association of Geodesy (IAG) e dall'International Union of Geodesy and Geophysics (IUGG), si basa attualmente sulle osservazioni di quattro tecniche di geodesia spaziale : Very Long Baseline Interpherometry (VLBI), Satellite Laser Ranging (SLR), Global Navigation Satellite System (GNSS) e Doppler Orbitography and Radio-positioning Integrated by Satellite (DORIS). Le diverse reti sono connesse rilevando la distanza a Terra tra i punti di riferimento di strumenti di tecniche differenti. Le misure possono essere realizzate solo presso siti multitecnica, la cui distribuzione spaziale è altamente disomogenea. Inoltre, a causa della natura manuale del rilevamento, le campagne non vengono ripetute frequentemente. A causa degli errori sistematici che affliggono potenzialmente tutte le misure coinvolte, poi, discrepanze significative sono spesso osservate tra le stime da tecniche spaziali ed i vincoli valutati a Terra. L'indagine di metodi di combinazione alternativi risulta dunque di primario interesse per lo sviluppo di quadri di riferimento di sempre maggior accuratezza e precisione.

In questo lavoro è stato studiato il potenziale della combinazione di misure GNSS ed SLR attraverso la loro co-locazione a bordo di satelliti GNSS equipaggiati con specchi retro-riflettori. L'analisi è stata effettuata sia per mezzo di realizzazioni quasi istantanee del quadro di riferimento che tramite stime di lungo periodo. É inoltre stata condotta una serie di simulazioni per stabilire l'effetto apportato alla precisione dei parametri fondamentali del quadro di riferimento da miglioramenti tecnologici ed infrastrutturali. Nell'interpretazione dei risultati, occorre tener presente che tutte le nuove strategie di combinazione devono avere come obiettivo il miglioramento della realizzazione dei parametri fisici del quadro, ovvero origine e scala, da cui dipendono criticamente le applicazioni scientifiche e civili. I risultati ottenuti mostrano che, nelle condizioni operative attuali, i vincoli spa-ziali a bordo dei satelliti GNSS non possono essere utilizzati per rimpiazzare le connessioni a Terra. Nella combinazione basata sulla co-localizzazione in orbita, infatti, l'informazione sulla scala e sull'origine non viene mai trasferita da una rete all'altra. Confronti settimanali hanno mostrato che i due sistemi indipendenti, GNSS ed SLR, mantengono inalterata la propria realizzazione del quadro all'interno della soluzione combinata. In particolare, l'origine e la scala del GNSS non beneficiano della combinazione con le osservazioni LAGEOS capaci di fornire una realizzazione più accurata del quadro.

Un'ulteriore indicazione che la co-localizzazione a bordo dei satelliti GNSS è insufficiente per sintetizzare un quadro organico a partire dai dati GNSS e dalle osservazioni SLR ai satelliti GNSS è fornita dall'analisi dei risultati della loro soluzione di lungo periodo. Tutti i confronti effettuati rispetto a tale soluzione, infatti, conducono a risposte incoerenti tra le due reti.

Infine, differenti scenari sono stati simulati aumentando il volume e la precisione delle osservazioni SLR ai satelliti GNSS. Persino nelle più favorevoli condizioni implementate, il guadagno osservato nella precisione dei parametri di riferimento è risultato marginale, soprattutto considerando lo scarto tra lo stato dell'arte del tracciamento laser e gli scenari implementati.

Alla luce dei risultati presentati in questa tesi, la co-localizzazione delle tecniche GNSS ed SLR a bordo dei satelliti GNSS non fornisce un solido approccio di combinazione per la realizzazione del quadro di riferimento terrestre. In questo contesto, i miglioramenti nella qualità e nella quantità delle osservazioni SLR sembrano in grado di produrre solo progressi limitati nel breve termine. La contribuzione più rilevante di tali osservazioni resta, pertanto, la validazione indipendente delle orbite stimate sulla base dei soli dati GNSS.

Résumé substantiel

Le but ultime de la géodésie est la détermination de la géométrie, du champ de gravité et de la rotation de la Terre ainsi que leurs variations dans le temps. Les études géodésiques concernent la surveillance et la modélisation des mouvements de la croûte, de la redistribution de masse entre les différents géosphères et des échanges de moment angulaire au sein et au-dessus de la surface de la Terre. L'évolution de ces trois champs est strictement inter-connectée et exige la réalisation d'un repère unique qui permet l'étude du système Terre dans son ensemble. Un tel repère est fourni par la définition d'un Système de Référence Terrestre (Terrestrail Reference System, TRS) et de sa matérialisation physique, un Repère de Référence Terrestre (Terrestrial, Reference Frame, TRF), qui permettent le positionnement uniforme des objets situé sur la surface de la Terre ou en orbite autour d'elle. La disponibilité d'un TRF précis et stable est si fondamentale pour les études scientifiques, les applications sociétales et les activités économiques que son importance ait été soulignée par l'Assemblée Générale des Nations Unies dans sa première résolution géospatiale adoptée le 26 Février 2015.

La détermination des coordonnées ITRF, qui est le repère officiel adopté par l'Association Internationale de Géodésie (IAG) et par l'Union Internationale de Géodésie et Géophysique (IUGG), repose actuellement sur les observations de quatre techniques de géodésie spatiale : Very Long Baseline Interpherometry (VLBI), Satellite Laser Ranging (SLR), Global Navigation Satellite System (GNSS) et Doppler Orbitography and Radio-positioning Integrated by Satellite (DORIS). Les différents réseaux sont connectés par des rattachements géodésiques qui font le lien entre les points de référence des instruments de mesures co-localisés dans un certain nombre de sites ITRF. La répartition spatiale de ces sites co-localisés est très inhomogène et leurs rattachements géodésiques sont rarement mis à jour. En outre, à cause des erreurs systématiques qui affectent potentiellement toutes les observations pertinentes, des écarts sont observés entre les liens terrestres et les estimations correspondantes des techniques de géodésie spatiale .

Dans ce travail, on a étudié une approche de lien alternatif, basée sur la colocalisation des techniques dans l'espace. En particulier, le potentiel de la combinaison des mesures GNSS et SLR à travers leur lien à bord des satellites GNSS munis de réflecteurs laser a été examiné. L'analyse a été effectuée par réalisations quasiinstantanée et à long terme du repère de référence. En outre, des perspectives futures ont été étudiées par simulations, pour quantifier ce que des améliorations technologiques et d'infrastructure pourraient apporter à la précision des paramètres fondamentaux du repère. Dans l'interprétation des résultats, il faut considérer que toutes les nouvelles stratégies de combinaison doivent viser l'amélioration de la réalisation des paramètres physiques du repère, à savoir l'origine et l'échelle, sur lesquels se fondent de façon critique un grand nombre d'études scientifiques et d'applications civiles.

Les résultats présentés montrent que, dans les conditions opérationnelles actuelles, les liens spatiaux à bord des satellites GNSS ne peuvent pas remplacer les rattachements terrestres dans la dérivation des repères de référence. Le lien spatial est, en fait, incapable de transmettre l'information de l'échelle ou l'origine d'un sousréseau à l'autre.

Une comparaison au niveau hebdomadaire a montré que les deux sous-réseaux indépendants, GNSS et SLR, conservent la réalisation du repère spécifique de leur technique au sein de la solution "combinée". En particulier, l'origine et l'échelle GNSS ne bénéficient pas du tout de la combinaison avec les observations LAGEOS qui fournissent une matérialisation plus pertinente des paramètres du repère.

Une autre indication que les liens spatiaux à bord des satellites GNSS sont insuffisants pour synthétiser un repère unique à partir de données GNSS et des observations SLR sur les satellites GNSS vient de l'analyse des résultats de leur cumul à long-terme. Toute les comparaisons avec cette solution conduisent à des résultats incohérents entre les sous-réseaux spécifiques de chaque technique.

En ce qui concerne les simulations, plusieurs scénarios basés sur l'incrément du volume et de la qualité des observations SLR sur les satellites GNSS ont été effectués. Même dans les meilleures conditions étudiées, le gain observé dans la précision des paramètres du repère de référence est marginal, surtout en tenant compte de l'écart technologique et d'infrastructure entre l'état de l'art de la télémétrie SLR et les scénarios envisagés.

Considérant les résultats présentés dans cette thèse, la co-localisation GNSS et SLR à bord des satellites GNSS ne fournisse pas une approche de lien solid pour la réalisation du repère de référence terrestre. Dans ce contexte, même des améliorations dans la qualité et la quantité des observations SLR aux satellites GNSS sont susceptibles d'avoir juste un impact limité dans un avenir proche. La contribution la plus pertinente de ces mesures SLR reste alors le processus de validation indépendante des orbites GNSS.

List of Tables

2.1 Ambiguity resolution strategies adopted according to baseline length.

Applicability to GLONASS system is also indicated. . . . . . . . . . . Chapter 1

Average values

Theoretical notes on Terrestrial Reference Systems and Frames

General context and outline

The ultimate goal of geodesy is the determination of the geometry, gravity field and rotation of the Earth and their changes in time [START_REF] Blewitt | Geodetic Observations and Global Reference Frame Contributions to Understanding Sea-Level Rise and Variability[END_REF]. Geodetic studies concern monitoring and modeling crustal movements, angular momentum exchanges within and above the Earth surface and mass redistribution between the different geospheres. The evolution of these three fields is strictly interconnected and demands for the realization of a unique framework enabling the study of the Earth system as a whole. Such framework is provided by the definition of a Terrestrial Reference System (TRS) and its associated physical materialization, a Terrestrial Reference Frame (TRF), which allow the consistent positioning of objects located anywhere on the Earth surface or orbiting around it.

In principle, a TRF provides "the foundation for nearly all ground-based and space-based observations" [START_REF] Minster | Precise geodetic infrastructure: national requirements for a shared resource[END_REF], fully enabling a multi-technique approach in the monitoring of geophysical phenomena. In addition, long-term reference frames allow tying together data acquired in successive measuring campaigns, which is a fundamental requirement for the reliable detection of slow movements and subtle crustal deformations [START_REF] Blewitt | Reference frame requirements for studies of geodynamics and climate change[END_REF]. As a consequence, advances in Earth System sciences and global change studies critically depend on the accuracy, stability and availability of the TRF. Interested processes encompass sea level change, variations in ice sheet thickness, ground subsidence, hydrological forcing, pre-seismic strain accumulation and visco-elastic characterization of Earth interior from postseismic deformation.

Moreover, TRFs are necessary to describe the evolution of Earth motion in space, since such movement can be generally expressed as a set of rotations between geocentric terrestrial and celestial frames. The accurate description of Earth rotation serves a number of purposes ranging from satellite tracking to supporting studies of mass movements and angular momentum exchanges between different components of the System Earth. Besides its support to scientific activities, a reliable TRF serves as basis for societal applications such as navigation, agriculture, land-use planning and natural hazard mitigation [START_REF] Blewitt | Geodetic Observations and Global Reference Frame Contributions to Understanding Sea-Level Rise and Variability[END_REF].

The scientific, social and economic importance of an accurate and stable TRF has been emphasized by the United Nation General Assembly within its first geospatial resolution adopted on February 26th, 20151 . The Assembly recognized the need of international cooperation for the realization of a Global Geodetic Reference Frame (GGRF) and urged Member States to implement open sharing of geodetic data, standards and conventions. The need for official national commitment to improving and maintaining appropriate geodetic infrastructure was stressed in the resolution. A dedicated working group, the United Nations Global Geospatial Information Management (UN-GGIM) Working Group, has also been formed to address the creation of a roadmap for the enhancement of the GGRF. The effective use of available resources should benefit from a strengthened international coordination addressing present day infrastructural limitations.

Together with the maintenance and development of observational facilities, progresses in data analysis and combination techniques are major contributors in the realization of TRF of increasing precision and stability. In particular, any innovative computation strategy should ameliorate the definition of the frame physical parameters, namely the origin and the scale, upon which a number of scientific applications critically rely. [START_REF] Kierulf | Precise point positioning requires consistent global products[END_REF], for example, showed that an error of 2 mm/y in the origin rate would produce an artifact mean sea level variation in the order of 0.4 mm/y, with local effects up to 3 mm/y. A frame scale drift error of 0.1 ppb/y would instead produce a drift up to 0.6 mm/y in mean sea level as determined by tide gauge records corrected using GPS data expressed in ITRF [START_REF] Beckley | A reassessment of global and regional mean sea level trends from TOPEX and Jason-1 altimetry based on revised reference frame and orbits[END_REF][START_REF] Collilieux | Global sea-level rise and its relation to the terrestrial reference frame[END_REF]. The geodetic community evaluated that the most demanding challenges require an accuracy of 1 mm and a stability of 0.1 mm/y in the realization of the frame [START_REF] Blewitt | Geodetic Observations and Global Reference Frame Contributions to Understanding Sea-Level Rise and Variability[END_REF]. As the ITRF2008 [START_REF] Altamimi | an improved solution of the international terrestrial reference frame[END_REF]] was still about one order of magnitude far from this goal [Altamimi et al., 2013a;[START_REF] Wu | Accuracy of the international terrestrial reference frame origin and earth expansion[END_REF][START_REF] Collilieux | External evaluation of the terrestrial reference frame: report of the task force of the iag sub-commission 1.2[END_REF], alternative/complementary computational schemes are worth investigating.

The present study aims at characterizing the frame realized by a combination of Satellite Laser Ranging (SLR) and Global Navigation Satellite Systems (GNSS) observations exploiting their co-location on-board GNSS spacecrafts (space tie).

The first chapter provides a general introduction on terrestrial reference systems and frames. The main features relevant for their definition are detailed and the mathematical formalism needed for their practical implementation is illustrated. The dissertation concentrates primarily on techniques and combination strategies currently exploited in the derivation of the International Terrestrial Reference Frame (ITRF), which is the reference frame officially adopted by the International Earth Rotation and Reference System Service (IERS) and by the Internal Union of Geodesy and Geophysics (IUGG). The space tie approach investigated in this study is also introduced.

The second chapter focuses on the GNSS and the SLR techniques and introduces the dataset selected to address the evaluation of the frame realized on the basis of their co-location on-board GNSS satellites. Technique-specific TRFs are presented and validated against the results that the official IERS services, namely the International GNSS Service (IGS) and the International ILRS Service (ILRS), submitted for the computation of the ITRF2014.

The derivation of accurate TRFs highlights the necessity of identifying discontinuities in geodetic time series in order to reliably assess linear trends. An automated blind procedure appointed for the detection of jumps in GNSS datasets is presented in App. C. The proposed methodology tailors to GNSS time series the application of the STARS algorithm [START_REF] Rodionov | A sequential algorithm for testing climate regime shifts[END_REF], originally conceived in the field of marine sciences. The performances of the devised screening strategy are tested against a synthetic data set and the application to a real case study is also presented [START_REF] Bruni | Detecting discontinuities in GNSS coordinate time series with STARS: case study, the Bologna and Medicina GPS sites[END_REF].

The third chapter, eventually, addresses the evaluation of a reference frame realized exploiting the space ties on-board GNSS vehicles. In particular, it is established how such a frame compares to the traditional ITRF implementation and what is the impact on the realization of the frame origin and scale. The study encompasses the computation of both long-term and quasi-instantaneous frames based on geodetic data collected by a global network of stations in the period 2011 -2014. Possible perspectives are also investigated on the basis of simulations contemplating a number of ideal observing scenarios.

Terrestrial Reference Systems and Frames

As "position" and "motion" are intrinsically relative concepts, a reference system must be introduced over a given domain in order to determine location and displacements of objects moving within it. From a mathematical point of view, the choice of a reference system is completely arbitrary and its determination is accomplished once its origin, scale and orientation are defined. To make the system accessible to users, however, it is convenient that there exists an observational relationship between its ideal points and a set of physical objects [START_REF] Kovalevsky | Reference frames: in astronomy and geophysics[END_REF]. The materialization of a defined reference system is provided by its associated reference frame which inherits the mathematical properties of the system and is realized by a set of physical markers with precisely determined coordinates [Petit and Luzum, 2010, Chapt. 4]. Once these coordinates are available, the system origin, orientation and scale are inherently defined. Locating an object in a particular reference system is then achieved by positioning such object relatively to the points which realize the frame.

According to the fourth chapter of the International Earth Rotation and Refrence System Service (IERS) Conventions [START_REF] Petit | IERS conventions[END_REF], "a Terrestrial Reference System (TRS) is a spatial reference system co-rotating with the Earth in its diurnal motion in space" . This statement implies that no rotations or collective translations of its points should be observable. If the Earth's crust is chosen as the physical structure supporting the definition of the system, the founding condition of a TRS is formalized by the system of equations:

C v dm = 0 C x × v dm = 0 (1.1)
where the integral domain C refers to the Earth's crust and x and v are, respectively, the position and the velocity of the infinitesimal mass element dm. When it comes to practice, however, it is impossible to survey the crust as a continuum. Therefore, the integrals of Eq. 1.1 must be replaced by summations over a limited number N of points distributed over the Earth's surface:

           N i=1 v i dm i = 0 N i=1 x i × v i dm i = 0. (1.2)
In order to verify the general property of the TRS, Eq. 1.2, the crust must be provided with a modeling describing the relationship between the structure configuration and the actual value of the coordinates associated to its points [START_REF] Kovalevsky | Reference frames: in astronomy and geophysics[END_REF]. The characterization of the reference system thus require the definition of a number of constants, conventions and theoretical models allowing the description of crustal physical evolution with time. Once this modeling is designated, point coordinates in the specific TRS are fully defined and can be computed for the fundamental points realizing the materialization of the system. Users can then retrieve the position of any other object within the relevant domain by evaluating its ties to such frame markers, or by exploiting ITRF-based products for datum definition purposes (e.g. using IGS GNSS orbits when performing PPP analysis).

Both the TRS and the TRF computed according to the presented procedure are known as "conventional". On the one hand, the definition reflects the arbitrariness in the choice of the structure modeling, which is far from being unique; on the other side, the means exploited in the practical realization of the frame (e.g. observing techniques, combination strategies...) are also selected among different possible options. As will be detailed in the following sections, the reference system officially adopted by the IERS is based on the resolutions of international scientific unions which endorse specific models, procedures and constants to be used in the data analysis and combination. The full reference to such conventions are found in [START_REF] Petit | IERS conventions[END_REF].

IERS Conventions on Terrestrial Reference Systems and Frames

The International Earth Rotation and Reference Systems Service (IERS) was established in 1987 by the International Astronomical Union (IAU) and the International Union of Geodesy and Geophysics (IUGG) with the aim of computing and maintaining the International Terrestrial and Celestial Reference Frames and their relative rotations as a function of time. On the basis of constantly updated scientific contributions, the IERS provides data, models and products which serve as a standard for the activities of the astronomical, geodetic and geophysical communities.

TRS Specifications

In 2010, the latest edition of the IERS Conventions was published collecting findings and recommendations of the different international scientific unions. Concerning the definition of the Terrestrial Reference System [Petit and Luzum, 2010, Chapt. 4], it is established that a TRS is a tridimensional Euclidean space (O, E), which is therefore fully characterized by the identification of its origin, scale and orientation. According to the convention, the origin O is close to the Earth's center where δ ij is the Kronecker symbol. The norm of each basis versor represents the scale of the space and it is isotropically close to the S.I. meter. Concerning the orientation, the Z-axis is aligned with the Earth's rotation axis, while X-and Y-axes lie in the Equatorial plane.

The transformation between different reference systems follows the 3-dimensional similarity transformation:

X 2 = T + kRX 1 (1.4)
where X 2 and X 1 are the coordinate vectors of the same point in the two TRS, T = [T x , T y , T z ] T represents the translational vector between the origin of the two systems, k expresses the ratio between the relevant scales and R describes the 3-dimensional rotation needed to align the bases of the two spaces. A visual representation of these transformation parameters is presented in Fig. 1.1, from [START_REF] Pollet | Combinaison des techniques de géodésie spatiale: contributions aux réalisations des systèmes de référence et à la détermination de la rotation de la Terre[END_REF].

In general, R is given by the product of the rotational matrices around the x-, y-and z-axis:

R = R x R y R z (1.5) with R x =   1 0 0 0 cos r x -sin r x 0 sin r x cos r x   R y =   cos r y 0 sin r y 0 1 0 -sin r y 0 cos r y   R z =   cos r z -sin r z 0 sin r z cos r z 0 0 0 1   (1.6)
According to [START_REF] Altamimi | Systeme de reference terrestre: definition, realisation, application a l'ITRF[END_REF], second and higher order terms can be neglected in Eq. 1.4 as long as the rotation angles between the two systems are smaller than 10 -5 rad and the scale factor ratio is smaller than 10 -5 . In addition, under such conditions it holds:

lim θ->0 cos θ = 1 lim θ->0 sin θ = θ (1.7)
and the rotational matrix R can be approximated by:

R = I + R = I +   0 -r z r y r z 0 -r x -r y r x 0   (1.8)
Eq. 1.4 can then be linearized into

X 2 = X 1 + T + dX 1 + RX 1 (1.9)
where the scale factor between the two TRSs is rewritten in the form k = d + 1.

In general, the exploited transformation parameters are function of time in order to account for possible deviations between the temporal evolution of the two frames.

TRF Specifications

As previously introduced, a reference system is a purely mathematical concept; it can be accurately defined, but the location of its axes is not directly measurable. Actual access to a given TRS is provided by its associated frame (see section 1.3.1 for details concerning the derivation of frame coordinates).

Terrestrial reference frame realizations can be divided in two main categories according to the time span covered by the data exploited in their computation. Quasi-instantaneous reference frames are realized over short (up to one week) periods, and give direct access to point coordinates at the representative epoch of each interval. Long-term frames, on the other hand, typically span several years and provide mean coordinates at epoch t 0 together with the set of parameters necessary to propagate them to any generic epoch t. In the simplest case, only linear displacements Ẋ are considered and frame point coordinates can be propagated according to X(t) = X 0 (t) + Ẋ • (t -t 0 ) [Altamimi et al., 2013a]. Since the long-term approach is followed in the derivation of the official IERS product, this dissertation privileges a presentation of such frames which are of prominent interest for Earth Science applications as they allow monitoring the evolution of geophysical phenomena over time. Nevertheless, it is worth mentioning that short-time implementations naturally account for non-linear motions and, at every epoch, are referred to the quasi-instantaneous realization of the Earth's center of mass.

For long-term crust-based TRFs, a piece-wise linear parametrization of point coordinates has been historically implemented. Such modeling accounts for instantaneous changes in point position or velocity (e.g. due to seismic events, monument instabilities, instrumentation changes...) and for most of the effects of litosphere plate motions and secular geophysical signals, but it is obviously not suitable for a comprehensive description of the complex crustal movements. On the other hand, the reference coordinates of TRF points shall be as representative as possible of the actual marker positions in order to ensure a consistent realization of the frame parameters over time. For this reason, the effects of well described geophysical non linear signals are usually hived off in the adjustment of geodetic observations, leading to the definition of a set of coordinates which vary more regularly with time. At any given epoch t, the link between the conventional X(t) and the regularized X R (t) coordinates of a TRF point is given by:

X(t) = X R (t) + i ∆X i (t) (1.10)
where the corrective terms ∆X i (t) account for the so-called conventional displacements and are calculated on the basis of well established models. The IERS Conventions [Petit and Luzum, 2010, Chapter 7, and references there in] currently recommend to account for:

-solid Earth tides: elastic response of the crust to the external tide generating potential (TGP) produced by the Moon, the Sun and other planets of the Solar System; the tides act mainly on a daily and sub-daily scale and can reach several tens of cm -tidal ocean loading deformations: periodic deformations of the Earth's crust induced by temporal variation of the load associated to the tidal redistribution of oceanic mass; act mainly on a daily and sub-daily scale and can reach several cm -S 1 and S 2 atmospheric pressure loading deformations: crustal periodic motions caused by surface pressure oscillations which follow the diurnal heating cycle of the atmosphere; S 1 and S 2 refer to the diurnal and semidiurnal frequencies respectively and vertical effects are up to a few mm -rotational deformations due to polar motion: local displacements generated by the changes in the gravitational potential produced by variations in the direction of the Earth's rotation axis; station coordinates might be affected in the order a few centimeters. The stated changes in gravitational potential also affect the distribution of oceanic masses and, consequently, the associated loading. Models available to account for the ocean pole tide loading deformations are also recommended in geodetic data analysis.

Accounting for the recommended corrective terms (i.e. switching to regularized coordinates) improves the agreement between the assumed piece-wise linear modeling and the temporal evolution of frame point coordinates. Nevertheless, the presented parametrization cannot, under any circumstances, accommodate nonconventional non-linear motions. Discrepancies between actual marker locations and their associated regularized coordinates are therefore expected and relevant deviations may be observed in presence of strong non-linear signals (e.g. post seismic relaxation effects). The most accurate realizations of terrestrial reference frames have already gone beyond the piece-wise linear parametrization described in the IERS Conventions and account for post-seismic displacements and residual seasonal signals induced by atmospheric, oceanic and hydrological non-conventional loadings [Altamimi et al., 2013b]. Altamimi et al. [2013a] and [START_REF] Collilieux | Impact of loading effects on determination of the international terrestrial reference frame[END_REF], however, demonstrated that estimating non-linear seasonal displacements does not impact the values and the associated linear trends of the realized ITRF parameters provided that several years of data are available.

International Terrestrial Reference System and Frame

Since 1991 and 2007 respectively, the International Association of Geodesy (IAG) and the International Union of Geodesy and Geophysics (IUGG) formally adopt a particular TRS and TRF, namely the International Terrestrial Reference System (ITRS) and Frame (ITRF). Quoting from the Conventions,the ITRS definition fulfills the following conditions:

1. it is geocentric, its origin being the center of mass for the whole Earth, including oceans and atmosphere;

2. the unit of length is the meter (S.I.) [...];

3. its orientation was initially given by the BIH orientation at 1984.0 (meaning that the ITRS orientation at epoch 1984.0 is that of the Bureau International de l'Heure (BIH) reference system, Ed); 4. the time evolution of the orientation is ensured by using a no-netrotation condition with regards to horizontal tectonic motions over the whole Earth.

The system so defined is made accessible to users through the computation of the associated ITRF. Due to the continuous deformations undergone by the Earth's crust, the accuracy of any specific ITRF is entailed to degrade with time. Its realization must then be periodically updated in order to account for newly acquired observations and for upgrades in data analysis procedures and/or combination techniques. Thirteen editions have been issued by the IERS since its establishment in 1988; the last one, ITRF2014, was released in January, 2016. The ITRF is a long-term frame, which, until the penultimate realization ITRF2008, used to be characterized by a piece-wise linear parametrization of coordinates. Together with station positions and velocities, time series of daily Earth Orientation Parameters (EOPs) are also made available to users since ITRF2005. ITRF2014 additionally encompasses the modeling of post-seismic relaxations [START_REF] Lercier | Evaluation of parametric post-seismic models and application in reference frame determination[END_REF]; in the preparation of the frame, station-specific seasonal signals have also been estimated in order to improve the reliable estimation of long-term linear trends, but such signals are not included in the official model for coordinate propagation.

A multi-technique approach is necessary to achieve present day accuracy and stability of the frame. Combining different types of measurements, in fact, allows compensating for system-specific biases and drawbacks with the strengths of the other observing methods. Starting from ITRF94, the computation of frame positions, velocities and of the relevant full covariance matrix depends on the observations of four space geodetic techniques: Very Long Baseline Interferometry (VLBI), Satellite Laser Ranging (SLR) to LAGEOS and ETALON satellites, Global Positioning System (GPS) and Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS). ITRF points represent the reference coordinates of the stations participating to the technique-specific observing networks. A brief introduction on the contributing systems is provided in the next chapter, with a particular focus on the SLR and GNSS techniques which are exploited for the derivation of the results presented in this thesis. The role of the different techniques in the definition of the TRF datum is instead detailed in Sec. 1.3.3.

Together with space geodetic observations, the availability of a set of terrestrial ties is required in the computation of the ITRF. Expressing the distance and the angles between the reference markers of the different techniques at those sites where multiple observing systems are co-located, these vectors allow connecting the diverse system-specific networks to each others. More details on local ties evaluation and their role in the computation of TRF are provided in Sec.1.4.1. Eventually,next Sec. 1.3.1 illustrates the computation strategy followed in the preparation of ITRF2008 and exploited in the derivation of the long-term TRF calculated for this work.

ITRF computation

Since ITRF2005, the ITRF computation stems from the combination of time series of space geodetic solutions. In order to construct an optimal set of coordinates and velocities, input information must be consistently exploited. Following [START_REF] Altamimi | The choice of reference system in ITRF formulation[END_REF], it is considered that each input technique-specific solution at epoch t k materializes its own realization of the reference system. Therefore, assuming a piece-wise linear parametrization for the definition of the frame, the observed coordinates x T,i (t k ) of the i -th station from technique T can be expressed as a function of the relevant unknown TRF products, (x 0,i , v i ), and of the transformation parameters p T,k (see Sec. 1.2.1.2) linking the combined reference frame to the one materialized for technique T at epoch t k :

x T,i (t k ) = f T,i (p T,k , x 0,i , v i ) .
(1.11)

In order to achieve the highest possible accuracy, the ITRF is formulated as the result of the simultaneous stacking of the four contributing techniques, where the term stacking designates "the determination of initial coordinates and velocities with simultaneous transformation of every epoch coordinates to a common reference system" [START_REF] Altamimi | The choice of reference system in ITRF formulation[END_REF]. The procedure is, however, computationally prohibitive, and a two-step process is generally implemented [Altamimi et al., 2007a].

In a first phase, long-term reference frames are computed for each individual technique; subsequently these frames are combined into a unified solution taking advantage of the terrestrial ties information. Provided there is no dependency between the observations from different techniques, the two steps approach is mathematically equivalent to the completely simultaneous stacking [Altamimi et al., 2013a]. The proposed combination model as implemented in the Combination and Analysis of Terrestrial Reference Frames (CATREF) Software package [START_REF] Altamimi | A new release of the international terrestrial reference frame for earth science applications[END_REF][Altamimi et al., , 2005a,b] ,b] has been exploited for achieving the results presented in this thesis.

Step 1: long-term stacking of technique-specific time series

In the first step, technique-specific time series, typically derived from daily/weekly observations, are accumulated into their relative long-term frame. For the computation of the ITRF, input datasets result from the combination of different Analysis Centers (ACs) solutions. Each IAG Service in charge of coordinating a specific technique contribution is responsible to appoint Combination Centers (CCs) that synthesize the processing efforts of the official ACs into a unique adjustment. Once this preliminary combination is finalized, the resulting time series are submitted to the ITRS Center and to the IERS Central Bureau. Details concerning combination strategies and main findings of the different CCs in the framework of the ITRF2014 call for participation can be found in [START_REF] Bachmann | Ivs contribution to the next ITRF[END_REF], [START_REF] Rebischung | Combination and stacking of the igs repro2 terrestrial frames[END_REF], [START_REF] Luceri | Extending the ilrs terrestrial reference frame development contribution to ITRF2014[END_REF], [START_REF] Moreaux | Ids combined solution improvements between ITRF2008 and ITRF2013[END_REF] for the International VLBI Service (IVS), the International GNSS Service (IGS), the International Laser Ranging Service (ILRS) and the International DORIS Service (IDS) respectively.

For each selected sampling interval (i.e. daily for GNSS, weekly for SLR and DORIS and session-wise for VLBI), input time series must provide the set of parameters and their associated variance-covariance matrix best fitting the data collected by ground stations. The reference frame information inherent to the available observations can be accessed only with the simultaneous estimation of all the quantities affecting the technique-specific observables. Failing to account for all of them, would mean fixing some of the parameters to a priori values expressed in an external reference frame which would then be aliased, to some extent, in the new realization. Data adjustment should then encompass:

-station coordinates -satellite orbit parameters for GNSS, SLR and DORIS/ quasar positions for VLBI -Earth Orientation Parameters (EOPs) describing the Earth's orientation in space, see below -delays associated to the signal propagation through the atmosphere -technique-specific biases absorbing systematic errors and mismodeling

As only the official ITRF products, namely station coordinates and EOPs, are actually combined in the stacking, the other quantities may be pre-eliminated [START_REF] Dach | Manual of Bernese GPS Software Version 5.0[END_REF] for computational convenience. Pre-eliminated parameters are implicitly solved for at the level of each time series sample, but no long-term combined values will be provided at the end of the TRF adjustment.

The long-term stacking of technique-specific time series is then realized with the combination model [Altamimi et al., 2007a]:

X i s = X i c + (t s -t 0 ) Ẋi c + + T k + D k X i c + R k X i c (1.12)
which refers to a total network of N points, i = 1...N , and assumes the availability of a set of S solutions, s = 1...S, encompassing station positions X i s at epoch t s expressed in a TRF k. The long-term accumulation aims at estimating:

-station positions X i c at a given epoch t 0 and velocities Ẋi c expressed in the combined frame c -transformation parameters P k between the combined TRF and each input frame k. In the most general formulation, P k encompass a translation vector T k , a scale factor D k and a rotation matrix R k written according to Eq. 1.8. If each solution s is expressed in its own reference frame, it is implied that S = K.

The combination of the available EOPs is performed according to:

                   x p s = x p c + r yk y p s = y p c + r xk U T s = U T c -1 f r zk ẋp s = ẋp c ẏp s = ẏp c LOD s = LOD c (1.13)
where:

-r xk , r yk , r zk are the non zero elements of the rotation matrix R k -f is the conventional rate of advance of the Earth Rotation Angle (ERA) with respect to the Universal Time U T 1 -x p s , y p s , ẋp s and ẏp s represent pole coordinates and their daily associated drifts, while U T and LOD account for deviations and relative rates between the Universal Time U T 1 and the Coordinated Universal Time U T C. The equations related to polar motion rate are introduced in order to handle input files where pole coordinates are expressed in the form of offset and drift; if a continuous piecewise linear function is used, fourth and fifth lines of Eq. 1.13 are discarded.

Parameters listed in the last point, together with celestial pole coordinates, are needed to characterize the Earth's motion in space. In fact, if the space is provided with a reference system centered in the ITRS origin and with a fixed orientation towards distant celestial objects, the motion of the Earth can be described as the set of rotations demanded to align the ITRS to such system, commonly referred to as the Geocentric Celestial Reference System (GCRS). According to [START_REF] Petit | IERS conventions[END_REF], the necessary rotations can be parametrized through the previously introduced quantities, as

X GCRS (t) = Q(t)R(U T, t)W(x p , y p , t)X IT RS (t) (1.14)
where:

-Q(t) is a matrix describing the motion of the Earth's rotation axis w.r.t. the GCRS (precession and nutation)

-R(t) represents changes in the Earth's angular velocity in response to mass redistribution within the system. This quantity can be parametrized through the sum of a modeled component, deriving from tides and libration, plus the non-conventional contribution U T -W(t) accounts for the motion of the Earth's rotation axis with respect to the crust. As for the previous case, the effects of tides and libration can be suitably modeled, while the non-conventional components x p and y p must be estimated.

As a final remark on Eq. 1.13, it can be noted that the connection between the stacked frame and the EOP time series is realized by the three rotation parameters only. The combination of EOPs can therefore be neglected in investigations focusing just on the frame physical parameters (origin and scale). By so doing, the computation burden is considerably reduced.

Step 2: multi-technique combination of long-term solutions

When stacking the long-term solutions of the different techniques, the combination model previously introduced should be expanded in order to account for station velocities Ẋi s available in the input files [START_REF] Altamimi | The choice of reference system in ITRF formulation[END_REF]Altamimi et al., 2007a]:

           X i s = X i c + (t s -t 0 ) Ẋi c + + T k + D k X i c + R k X i c + (t s -t k ) Ṫk + Ḋk X i c + Ṙk X i c Ẋi s = Ẋi c + Ṫk + Ḋk X i c + Ṙk X i c (1.15)
where Ṫ, Ḋ and Ṙ represent the estimated temporal drifts of the transformation parameters between each input frame and the combined solution (ITRF). The combination scheme for EOPs, Eq. 1.13, holds unaltered.

In this second phase, the different input frames are totally distinct, meaning that they do not share any common points. Some sort of supplementary information must then be introduced in order to connect the various system-specific subnetworks. Rigorously speaking, then, the ITRF is computed on the basis of space geodetic observations complemented by the set of exploited inter-technique ties. Two alternative approaches for the realization of the network connection are presented in Sec. 1.4. 1.3.1.3 The role of transformation parameters Equations 1.12, 1.13 and 1.15 are written in the most general way, i.e. encompassing all possible transformation parameters between the different input frames and the combined one. Their simultaneous estimation is, however, non compulsory; according to specific analysis choices, limited subsets (or even none) of these parameters might be actually set up in the combination.

When performing the long-term stacking of technique-specific time series, estimating transformation offsets enables a critical evaluation of the quality of the input information [Altamimi et al., 2013a]. This approach allows evaluating the Weighted Root Mean Square (WRMS) of epoch results with reference to the combination solution, thus providing a robust check on the temporal variation of the internal precision of the original contributions. In addition, this formulation enables imposing internal minimal constraints (see Sec. 1.3.2.2) during the long-term stacking, ensuring that the resulting frame realization preserves the mean information expressed by the input time series. The analysis of the estimated transformation offsets, then, allows characterizing the temporal stability of the frame parameters. Particular attention should be devoted to the identification of discontinuities and unexpected trends which, if not properly handled, might impact the accuracy of the derived fully combined TRF. The datum definition selected for the ITRF, see Sec. 1.3.3, follows from such considerations on the quality of technique-specific frame realizations [START_REF] Altamimi | A new release of the international terrestrial reference frame for earth science applications[END_REF][Altamimi et al., , 2007a[START_REF] Altamimi | an improved solution of the international terrestrial reference frame[END_REF].

During the second step of the ITRF computation, setting up transformation offsets allows specifying which techniques should be responsible for the definition of the combined frame parameters. Each input long-term solution provides its own materialization of the reference system origin and scale. As a consequence of characteristic systematic errors, all these realizations are a biased version of the actual geocenter location and of the conventional S.I. meter. However, according to the sensitivity of the different techniques to the distinct degrees of freedom of the frame, some realizations are likely to be more accurate than others and should be used as a reference in the combined solution.

To exclude a certain technique from the definition of (a subset of) the combined frame parameters, it is sufficient to estimate its relevant similarity biases during the final stacking. By so doing, a correlation is introduced between the estimation of sta-tion positions and transformation offsets, see Sec. 1.3.2, and the frame information associated to that observing system is no longer uniquely defined. The realization of the corresponding combined frame parameters rely, therefore, on the mean of the contributions provided by the remaining techniques. For discarded systems, the estimated offsets absorb the inconsistencies between their original input solution and the resulting combined frame, computed according to the available ties. The uncertainties associated to transformation biases define the precision with which the reference information is transferred to the other techniques [Altamimi et al., 2013a].

Constraints

For each transformation parameter estimated during the stacking, a singularity is introduced. Translating all station coordinates of a constant amount along the X axis, for example, can be exactly compensated by countershifting all the input frames of an equal amount [START_REF] Rebischung | Can GNSS Contribute to Improving the ITRF Definition?[END_REF]. The impossibility to discriminate between these two explanations of the observed data results in a singularity of the combination system when X-translation offsets are set up for all the input solutions. Analogous considerations hold, obviously, for every transformation parameter.

To cope with these singularities, suitable constraints have to be provided to the combination system, Eq. 1.12, 1.13 and 1.15. One possibility consists in forcing station positions to some values known a priori. The strength of these absolute constraints is evaluated in terms of the maximum allowed discrepancy between a priori information and estimated results; tight constraints require an agreement at the level of 10 -10 m, while a 1 m accordance defines loose bounds. Such a conditioning heavily depends on the quality of the available a priori values, which is, however, difficult to assess at the accuracy level currently required to support geophysical applications. Moreover, relying on external solutions for station positions might distort the network adjustment, perturbing the information content expressed by the available set of geodetic observations. Alternatively, constraints might be imposed concerning the realized frame. If suitably implemented, this conditioning provides the minimum information necessary to overcome the system singularities, while preserving the optimal network shape defined by a weighted least square adjustment of the original observations [START_REF] Altamimi | The choice of reference system in ITRF formulation[END_REF]. Parameters estimated under such conditioning, in fact, still satisfy the unconstrained equations [Rebischung, 2014, appendix B.2.3], meaning that minimal constraints remedy exactly the rank deficiencies of the original system, while leading to a non distorted solution.

In the ITRF combination, the last approach is preferred. The CATREF Software package allows to set up both minimum (MC) and internal (IC) constraints (see Sec.1.3.2.3 and 1.3.2.2 respectively), depending on the choice/necessity to align the results to an external reference frame or to the mean information expressed by the input solutions.

Mathematical formulation

Mathematically speaking, constraints are expressed as relationships to which the estimated parameters have to obey. In case of linear constraints, as for the ITRF computation, they can be written in the form:

C T x = 0 (1.16)
where C T is full rank and its lines specify the constraining conditions on the parameters [Rebischung, 2014, appendix B.2]. A possible approach to account for these constraints during the combination, consists in adding them to the original system as properly weighted pseudo-observations. This means that the linear(ized) system Ax = l (1.17) derived from the original observations l associated to the weight matrix P l , is transformed into the system

A c x = l c (1.18) with A c = A C T and l c = l 0 (1.19)
associated to the weight matrix

P lc = P l 0 0 P c (1.20)
The weight matrices P l and P c are built from the uncertainties of actual and pseudo observations respectively. According to the metrics defined by P lc , the least square adjustment of Eq. 1.18 results in

A T P l A N + CP c C T Nc x = A T P l l b (1.21)
where N and N c are named normal matrix and normal matrix of the constraints respectively. If Im(N) + Im(N c ) = R p , then (N + N c ) is invertible and the solution of Eq. 1.18 is unique. According to [START_REF] Dermanis | Geodetic linear estimation technique and the norm choice problem[END_REF], the suggested metric ensures that the resulting parameter values are associated to minimum variances. The full covariance matrix associated to the solution is computed from:

Q = (N + N c ) -1 (1.22)
For the reader's convenience, an example of Eq. 1.17 is provided, concerning the observation equation for the long-term stacking of technique-specific time series. The system built from Eq. 1.12 and Eq. 1.13 can be written in matrix form as:

       x 1 . . . x S        =        J 1 (t 1 -t 0 ) J 1 v A 1 . . . . . . . . . J S (t S -t 0 ) J S v A S                   x c ẋc θ 1 . . . θ S            (1.23)
where:

-x s collects the estimated parameters of the s -th input solution

x s = . . . X i s x p s y p s U T s ẋp s ẏp s LOD s . . . T with i = 1...N and s = 1...S.
-the unknown combined parameters are organized in

x c = . . . X i c x p c y p c U T c ẋp c ẏp c LOD c . . . T and ẋc = . . . Ẋi c . . . T with i = 1...N .
θ s lists the transformation parameters between the s -th input solution and the combined one -J s = ∂x s /∂x c and J s v = ∂x s /∂ ẋc are matrices of one and zeros which associate each parameter x i in the s -th input solution to the corresponding combined parameter x c or ẋc respictively [START_REF] Rebischung | Can GNSS Contribute to Improving the ITRF Definition?[END_REF] -A s = ∂x s /∂θ s is the design matrix linking input data of a specific solution with the transformation parameters set up between that solution and the combined one:

A s =   A s XV A s EOP   (1.24)
where

A s XV =            1 0 0 X 1 s 0 Z 1 s -Y 1 s 0 1 0 Y 1 s -Z 1 s 0 X 1 s 0 0 1 Z 1 s Y 1 s -X 1 s 0 . . . . . . . . . . . . . . . . . . . . . 1 0 0 X N s 0 Z N s -Y N s 0 1 0 Y N s -Z N s 0 X N s 0 0 1 Z N s Y N s -X N s 0            A s EOP =         0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 -1/f 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0        
From a least square adjustment of Eq. 1.23 according to the metrics introduced by the weight matrix of the input observations, the normal matrix of the stacking can be computed. This matrix will present a rank deficiency equal to the number of transformation parameters set up in θ. The introduced singularities need to be compensated by introducing a suitable normal matrix of the constraints. For the long-term stacking of technique-specific time series, internal constraints are typically imposed to the definition of the frame physical parameters; the frame orientation, on the other hand, needs necessarily to be referred to an external frame realization.

Internal constraints (IC)

"The aim of using internal constraints to define the combined frame resulting from the stacking of time series is to preserve the intrinsic physical parameters, namely the origin and the scale of satellite techniques and the scale of VLBI, without involving external frames" [Altamimi et al., 2007a]. In practice, internal constraints are those that allow selecting the set of values which minimizes the sum of the squares of all unknown parameters [START_REF] Altamimi | The choice of reference system in ITRF formulation[END_REF]. The actual constraining condition can be derived from the consideration that a general transformation T with parameters θ changes the unknown p into p ′ = T (θ) p. If T can be linearized, the stated relation can be written as

p ′ = p + Eθ (1.25)
The desired condition p T p = min is then achieved when E T p = 0.

In the computation of the ITRF, the seven transformation parameters θ are assumed to vary linearly over time; at any epoch t k , the values of the transformation parameters can then be retrieved through a linear extrapolation in form of Eq. 1.25 where

E = I (t k -t 0 )I
The internal constraints conditions, then, result:

           K k=1 θ k = 0 K k=1 (t k -t 0 )θ k = 0 (1.26) 1.3.2.3 Minimum constraints (MC)
The combined frame can be aligned to an external reference known apriori. For some applications, relying on external references cannot be avoided. The orientation of the network of all the geodetic techniques concurring to the realization of the ITRF, for example, cannot be determined unambiguously when EOPs are estimated. As pointed out by [START_REF] Rebischung | Can GNSS Contribute to Improving the ITRF Definition?[END_REF], "any rotation of the station network can be compensated by variations of the EOPs without any impact on the geodetic observations". For the VLBI system, then, also determining the origin of the frame is intrinsically impossible given the differential nature of the basic observable. In all these cases, if the corresponding transformation parameters are set up in the combination model, the estimated values must necessarily be referred to a pre-existing frame. In addition, minimum constraints can be useful to compare different frame realizations, as the estimated parameters would automatically highlight the effects produced by the different options/data exploited in the derivation of the input solutions.

When the relative transformation parameters are assumed to evolve linearly over time, the standard relation between two TRFs is given by [START_REF] Altamimi | CATREF software: Combination and analysis of terrestrial reference frames[END_REF]:

X 2 = X 1 + Aθ (1.27) where θ = T x , T y , T z , d, r x , r y , r z , Ṫx , Ṫy , Ṫz , ḋ, ṙx , ṙy , ṙz and A =             • • • • • • • • • • • • • • 1 0 0 x 0 i 0 z 0 i -y 0 i 0 1 0 y 0 i -z 0 i 0 x 0 i 0 0 0 1 z 0 i y 0 i -x 0 i 0 1 0 0 x 0 i 0 z 0 i -y 0 i ≈ 0 0 1 0 y 0 i -z 0 i 0 x 0 i 0 0 1 z 0 i y 0 i -x 0 i 0 • • • • • • • • • • • • • •             (1.28)
Forcing one TRF onto the other is equivalent to ask for zero transformation offsets between the two frames. As a least square adjustment of Eq. 1.27 leads to:

θ = A T A -1 A T B (X 2 -X 1 ) (1.29)
the desired constraints are realized by imposing:

B (X 2 -X 1 ) = 0 (1.30)
In terms of normal matrix, the same condition reads:

B T Σ -1 θ B (X 2 -X 1 ) = 0 (1.31)
where the weight matrix Σ θ , typically diagonal, sets the strength of the desired datum alignment.

Datum definition

The strategy described for the ITRF computation aims at realizing the best possible materialization of the ITRS defined by the IERS Conventions (Sec. 1.3). The proposed combination model preserves the quality characterizing space geodetic observations and the consistency of their associated technique-specific frame. In addition, it allows selecting the techniques to exploit in the implementation of the frame origin, scale and orientation (datum definition). A careful assessment of the quality of the technique-specific materializations is therefore necessary to select the best combination strategy. Choices adopted for the realization of ITRF2014 are here briefly illustrated.

According to the official ITRF website [IERS, http://itrf.ensg.ign.fr/ITRF_ solutions/2014/frame_ITRF2014.php], the origin of ITRF2014 is inherited from the SLR observations to LAGEOS and ETALON satellites. From the computational point of view, this has been achieved estimating translation parameters and rates for all input frames contributing to the multi-technique stacking, but SLR. The VLBI contribution has not been taken into account as its basic observable is intrinsically insensitive to displacements of the frame origin. The choice of rejecting the other satellite techniques is, instead, less obvious. In principle, all of them are sensitive to the geocenter position determination because the equations of motion governing the orbit determination process hold only in a frame centered on the Earth's center of mass. Concerning the GNSS, however, Rebischung et al. [2014] has shown that the simultaneous estimation of station and satellite clock offsets together with tropospheric parameters leads to collinearity issues that strongly degrade the quality of the sensed geocenter motion. The IDS report 2014 [Soudarin, 2015, Session 3.d.iv -Fig. 2], then, clearly showed that the DORIS system is still providing a poorer origin (and scale) information compared to the SLR.

The ITRF2014 scale and relative rates are realized by the averaged contributions of SLR and VLBI observations. The DORIS technique is once again disregarded because of its problematic temporal behavior. A clear discontinuity of about 1 cm is still observed from mid 2012, and a previously detected offset in 2002 was just recently recovered thanks to beacon frequency offset estimation [Soudarin, 2015, Session 3.d.iv]. GNSS, on the other hand, is long known to be weakly sensitive to the scale due to the correlation between mean station heights, clock offsets, tropospheric parameters and satellite antenna phase center offsets (z-PCOs) for which no pre-launch characterization is available [START_REF] Zhu | Detecting discontinuities in GNSS coordinate time series with STARS: case study, the Bologna and Medicina GPS sites[END_REF][START_REF] Rebischung | Can GNSS Contribute to Improving the ITRF Definition?[END_REF].

The orientation of the frame is, instead, purely conventional since the estimation of EOPs correlates with that of rotational parameters [START_REF] Rebischung | Can GNSS Contribute to Improving the ITRF Definition?[END_REF]. Rotation offsets and rates are, therefore, set up for all the four contributing solutions and the generated singularity is recovered by setting up minimum constraints with reference to the previous ITRF official realization [START_REF] Altamimi | an improved solution of the international terrestrial reference frame[END_REF]. Constraints are applied over a core network of 127 stations distributed among 125 ITRF sites.

As a general remark, it should be noted that the chosen coordinate parametrization directly impacts the agreement level between the materialized frame parameters and the ITRS specifications issued by the IERS. In particular, if a linear evolution model is assumed for station coordinates, all non linear variations of the Earth's center of mass cannot be reproduced [START_REF] Rebischung | Can GNSS Contribute to Improving the ITRF Definition?[END_REF]. This misalignment might pose some problems in the unambiguous detection of subtle geophysical signals which would benefit from the direct access to the quasi-instantaneous location of the geocenter. Impacted studies concerns, for example, the quantification of the sea-level rise trends or the identification of present-day ice melting fingerprints [START_REF] Beckley | A reassessment of global and regional mean sea level trends from TOPEX and Jason-1 altimetry based on revised reference frame and orbits[END_REF][START_REF] Collilieux | Global sea-level rise and its relation to the terrestrial reference frame[END_REF]. For the purposes of this thesis, however, a linear parametrization is still acceptable as the main conclusions will be derived by the relative comparison of reference frames realized exploiting different tying approaches. Provided that the same parametrization is exploited in both combinations, the conclusions drawn from the comparison are trustworthy and can be used to address the problem under study.

Inter-technique ties

In order to establish a connection between distinct technique-specific networks, it is necessary that different geodetic systems are co-located within a distance that can be measured with independent high precision strategies. On ground, co-locations are realized at multi-technique stations, where different instruments are installed within an area which can be surveyed with classical geodesy approaches; alternatively, satellites trackable by different systems can be exploited to realize the link in space.

Terrestrial ties

In the traditional ITRF computation, the four contributing techniques are linked through the so called local ties, which measure the distance and angles between the different observing instruments installed at those ITRF sites hosting more than one system. These three-dimensional vectors and their associated uncertainties are evaluated with terrestrial surveying campaigns performed by geodetic agencies. For this reason, in the context of this thesis, local ties will be rather referred to as terestrial ties, to avoid confusion with the alternative possible link, see Sec. 1.4.2. In the computation of the ITRF, equality constraints are also imposed on the velocities of co-located points. The static connection of the diverse networks realized by introducing the available ties, is then propagated over time bounding the displacements of co-located stations.

Quoting from Altamimi et al. [2007a], "the very existence of the ITRF relies on the availability and quality of local ties in co-location sites as well as the number and distribution of these ties over the globe." Nevertheless, the location of multitechnique ground sites is still far from being optimal, as can be seen in Fig. 1.2 that maps co-located stations active in 2014. Most of the connections are realized in the northern emisphere and tend to be clustered in dense groups rather than being homogeneously distributed over different countries. The situation is particularly unbalanced for the SLR and the VLBI networks. These techniques play a fundamental role in the definition of the ITRF parameters, see Sec. 1.3.3, but they are long known to lack adequate global coverage. In addition, Fig. 1.2 clearly shows the importance of GPS/GNSS in tying together the different technique-specific frames. GPS stations, in fact, are involved in about 80% of the total available ground colocations: of 56 co-locations involving a DORIS beacon, 43 encompass a GPS/GNSS antenna and the ratio increases to 39/48 and 27/32 for VLBI and SLR respectively.

Besides their suboptimal geographical distribution, terrestrial ties are also rarely updated because site surveying is performed in episodic dedicated campaigns rather than through automated and continuous measurements. On the other hand, it is known that monument instability [START_REF] Langbein | Noise in GPS displacement measurements from Southern California and Southern Nevada[END_REF][START_REF] Haas | Evaluation of GNSS monument stability[END_REF], and even intrasite movements [Sarti et al., 2013b] might affect the internal geometry of co-located sites, leading to a progressive degradation of the tie quality over time.

Further issues arise from the actual measurement strategy followed during the evaluation of the ties. According to Sarti et al. [2013a], a major problem is posed by the lack of standard approaches and critical comparisons of the different sitedependent methods. The whole procedure typically encompasses four stages, all potentially subjected to systematic errors impacting the accuracy of the final results [Sarti et al., 2013a, and references therein]: 1. surveying: in this phase, measurements of the distance and angles between the reference points (RPs) of the different instruments are acquired, typically by means of high-precision total stations. The subtlest aspect to consider is that terrestrial surveys can only access physical structures. As a consequence, they can just be referred to a Conventional RP, determined by the instrument geometry. Space geodetic techniques, on the other hand, acquire their observable in the Electronic RP, so that a gap arise between the surveyed vectors and the actual quantity required in the computation of the ITRF. A further assessment of the stability of the tie between the Conventional and the Electronic RP is needed to fully characterize terrestrial ties. In the case of the GNSS network, for example, this implied a considerable calibration effort to evaluate absolute antenna phase center corrections [START_REF] Schmid | Absolute phase center corrections of satellite and receiver antennas[END_REF]. In addition, survey observations might get corrupted by atmospheric artifacts or calibration issues.

2. statistical reduction of collected observations and data analysis: this step is considered the least critical as it is usually performed with well-established softwares performing in comparable ways [START_REF] Dawson | Indirect approach to invariant point determination for SLR and VLBI systems: an assessment[END_REF]. It is important to stress that the full variance/covariance matrix shall be provided in order to ensure proper weighting of the tie information in the ITRF computation.

3. data conditioning: when surveying complex structures such as SLR and VLBI antennas, measurements are directed to a series of targets expressly attached to the instrument and the location of the Conventional RP is retrieved adjusting the observations while imposing suitable geometric conditions. Depending on the different exploited observation geometries, however, biases up to 3 mm can be observed [START_REF] Dawson | Indirect approach to invariant point determination for SLR and VLBI systems: an assessment[END_REF].

4. alignment of the terrestrial ties to the global frame: at the end of phase 3, terrestrial ties are completely estimated, but they are expressed in the local topocentric frame. In order to be consistent with the results of space geodetic techniques, and useful for the computation of the ITRF, ties must be transformed into a global frame. [START_REF] Ray | Evaluation of co-location ties relating the VLBI and GPS reference frames[END_REF] showed that different aligning approaches commonly used in the survey practice might lead to discrepancies up to several millimeters in the transformed vectors.

Terrestrial ties are introduced in the ITRF combination as independent measurements, accounting for their full variance/covariance matrix (when available) and for the surveying epoch. Weights are then adjusted in an iterative process based on the agreement between each tie and the corresponding space geodetic estimates [START_REF] Altamimi | A new release of the international terrestrial reference frame for earth science applications[END_REF][Altamimi et al., , 2007a]]. Statistically significant inconsistencies exceeding the relevant formal errors are often observed. [START_REF] Altamimi | an improved solution of the international terrestrial reference frame[END_REF] showed that half of the VLBI/SLR-to-GPS ties exploited in the computation of ITRF2008 disagreed with space geodesy estimates up to 6 mm and about 30% of the cases presented discrepancies at the centimeter level. It is important to stress that systematic errors in the results of spatial techniques [START_REF] Ray | New global positioning system reference station in Brazil[END_REF]Seeber, 2003, e.g.] are likely to contribute to the observed dissimilarities and might even surmount possible errors in the derivation of the terrestrial ties. Nevertheless, the inhomogeneous spatial distribution and the rare revisit process of ground ties provide the rationale for the investigation of the alternative link available in space.

Space ties

Satellites that can be tracked by multiple systems provide an inter-technique linking opportunity alternative to the exploitation of terrestrial ties. This connection approach is usually referred to as space tie, since the different observing strategies are co-located on-board flying spacecrafts whose orbit parameters can be inferred from a joint combination of the diverse independent measurements. Technique co-location in space automatically overcomes some of the most problematic aspects of terrestrial ties, namely their spatial distribution and the frequency of their updates. Notwithstanding the technique with which they are acquired, all observations to a certain satellite can be gathered in the same adjustment encompassing orbit parameters and ITRF products. The inter-technique link is then realized continuously and every station of the observing networks may contribute in its realization, provided that it manages to track the satellites hosting the redundant positioning devices. Some critical aspects concerning the actual implementation of the space tie approach in the realization of the ITRF, however, have to be investigated. First of all, it should be noted that satellite orbit parameters describe the motion of the center of mass (CoM) of the spacecraft [Dach et al., 2007, Sec. 2.2], while the observations from space geodetic techniques are referred to the eccentric reference points of the relevant positioning payloads. Space ties are therefore realized by the vectors connecting these reference points to the satellite's center of mass (see Fig. 1.3 for a sketch representation). Accurate pre-launch evaluations of these quantities are however lacking for multi-technique spacecrafts currently available, posing a critical limitation on the link performances. In addition, abundance and quality of the system-specific observations play a major role in the effectiveness of the realized tie. Sparse and noisy observations might be insufficient to perform a link accurate enough to transfer frame information between the different techniques. Despite its conceptual elegance, the reliability of the inter-system connection via space ties depends on the assessment of the stated practical issues. In this thesis, the co-location between the GNSS and SLR techniques on-board GNSS satellites is taken into account for thorough evaluation. The investigated link is mainly provided by GLONASS satellites which are operationally equipped with Laser Retroreflector Arrays (LRA) and can therefore be tracked by SLR ground stations; concerning the GPS constellation, instead, only GPS-36 and GPS-35 can be exploited. In order to take advantage of this particular group of space ties, SLR observations to GNSS satellites have to be added to SLR and GNSS data currently contributing to the realization of the ITRF. As explained in Sec. 1.3.1.1, the reference frame information inherent to a certain set of observations can only be accessed adjusting the data over all the parameters impacting the technique-specific observables. In the adjustment of ranging measurements to GNSS satellites, it is necessary to set up, among others, orbital parameters and station-dependent quantities. The same elements must be estimated also while fitting GNSS observations and SLR shots to LAGEOS and ETALON satellites respectively. In a simultaneous adjustment of these three groups of observations, then, the joint estimation of common quantities provides a link between SLR and GNSS contributions without introducing relative inter-technique constraints on station positions (namely the terrestrial ties).

These considerations show that, in principle, the proposed space ties should be able to replace terrestrial constraints in the realization of the ITRF. The thesis aims at evaluating the actual strength of the link under the available operational conditions. In particular, it is evaluated how such an implementation of the reference frame compares to the traditional computation and it is established to which extent the provided tie is able to transfer to the GNSS network the origin and scale information as realized by LAGEOS observations. These aspects will be addressed from the perspective of both long-term and quasi-instantaneous reference frames. Diverse combination scenarios will be investigated encompassing various options for the datum definition. The comparisons between differently combined solutions and/or with technique-specific results (computed for reference purposes) allow evaluating the impact of the space tie approach on the definition of the frame origin and scale.

Within the literature, the technique co-location on-board GNSS satellites has already been proposed as a possible link for the combination of GNSS and SLR observations [START_REF] Thaller | Combination of GNSS and SLR observations using satellite co-locations[END_REF][Thaller et al., , 2012, among others], among others]. A critical and extensive assessment of the accuracy, precision and stability of the frame parameters derived with this approach was, however, still lacking. Earlier studies focused primarily on the impact that SLR observations to GNSS satellites brought to the derivation of station coordinates and on the possibility of validating the available set of terrestrial ties with the independent estimates provided by space geodetic techniques combined in orbit. A partial evaluation of the ability of the link to transfer origin and scale information is presented in the conference proceedings [START_REF] Thaller | Geocenter Coordinates from GNSS and Combined GNSS-SLR Solutions Using Satellite Co-locations[END_REF] and [START_REF] Thaller | Precombined GNSS-SLR Solutions: What Could be the Benefit for the ITRF[END_REF] respectively. Specific published results will be recalled in the following chapters and directly compared to the findings achieved in this study; it is already worth mentioning, however, that relevant space tie vectors, namely the GNSS Satellite Antenna Offset (SAOs) and the Laser Retroflector Array (LRA) offsets, are systematically recognized as poorly calibrated and in need of a dedicated re-evaluation within observation adjustment. Such estimation is, however, rather delicate due to the partial correlation between tie vectors and frame parameters.

Besides combination studies, SLR observations to GNSS satellites have been long exploited for orbit validation purposes [START_REF] Urschl | Calibrating GNSS orbits with SLR tracking data[END_REF][START_REF] Montenbruck | GNSS satellite orbit validation using satellite laser ranging[END_REF][START_REF] Sośnica | Satellite laser ranging to GPS and GLONASS[END_REF]. SLR tracking data are completely independent from microwave GNSS observations and are therefore entitled to provide a quality assessment of orbital products derived from the latters. The analysis of laser residuals is intended to reveal orbit modeling deficiencies especially in the radial component directly surveyed by ranging measurements [START_REF] Urschl | Calibrating GNSS orbits with SLR tracking data[END_REF]. According to the latest estimates [Sośnica et al., 2015, Table 6], SLR residuals to GPS and GLONASS orbits show average biases at the centimeter level when the newECOM model [START_REF] Arnold | CODEś new solar radiation pressure model for GNSS orbit determination[END_REF] is used in the microwave Precise Orbit Determination (POD). Offsets of this size are, however, difficult to be unambiguously interpreted as indicators of improper satellite positioning. Systematic SLR errors (mostly range biases and uncertainties in the values of LRA offsets) might, in fact, easily explain the observed results. In other words, orbit models are likely to have attained the maximum accuracy which can be validated with the state of the art knowledge of SLR characteristic response to GPS/GLONASS tracking. In addition, it should be observed that the statistical significance of the reported biases is also questionable, given that the RMS of SLR residuals is about 2 cm for GPS satellites and 4 -4.5 cm for GLONASS. Nevertheless, SLR observations have been valuable in the assessment of the newECOM model performances. Such parametrization allowed, in fact, to reduce the previously observed time variability of the residuals which, depending on the Sun-satellite elongation, was attributed to inaccurate modeling of solar radiation pressure.

The study presented in this thesis aims at evaluating whether the space ties between the SLR and the GNSS techniques realized on-board GNSS satellites provide a sufficiently strong link to substitute ground ties in the realization of Terrestrial Reference Frames. The ability of the link to transfer LAGEOS frame origin and scale information to a GNSS network is critically assessed.

The potential of the space tie approach has long been recognized by the scientific community, even though technical aspects concerning its actual implementation are known to pose serious concerns about the quality of the performed link. Issues highlighted in this study shall be addressed in the development of geodetic missions co-locating several positioning systems on-board the same spacecraft. It should be noted that the additional payloads required to allow redundant tracking of a satellite might pose mass, and consequently budget, issues in the realization of new space vehicles. A precise characterization of the actual benefits that can be brought by these supplementary links is therefore essential to provide sound recommendations to decision makers.

In the context of future geodetic missions implementing technique co-location in orbit, the Geodetic Reference Antenna in Space (GRASP) shall be mentioned. The mission concept consists of an oversimplified platform hosting the sensors of the four techniques currently exploited in the realization of the ITRF. The satellite design and an exhaustive pre-launch calibration shall ensure stable co-locations characterized at the millimeter-level. The positioning of the center of mass of the satellite in the inertial space shall also be accessible with an accuracy of 1 mm RMS, and stability of 0.1 mm/year [START_REF] Bar-Sever | Geodetic Reference Antenna in Space (GRASP) -a mission to enhance space-based geodesy[END_REF]. Just with such an optimized knowledge of the space tie vectors, the full value of technique co-location in orbit can be exploited. The great interest of the geodetic community is testified by the fact that the GRASP concept is currently submitted for evaluation both at NASA, within the framework of the "Earth Venture Mission" Call for Launch in 2020, and at ESA, possibly as a follow-up mission if the first proposal will be accepted.

Chapter 2 Techniques

The International Terrestrial Reference Frame is currently computed on the basis of the data of four space geodetic techniques and a set of terrestrial ties connecting the distinct networks. Each technique depends on specific observables which, together with the systematic errors affecting the measurements, determine the frame parameters to which observations are more sensitive. These considerations are reflected in the datum definition selected for the computation of the ITRF, see Sec. 1.3.3.

Space techniques currently exploited encompass:

-SLR: optical technique measuring the two-way travel time of an ultra-short laser pulse shot from a ground telescope to an orbiting retroflector and collected back at the original site. The measured time of flight is then translated into range records from which satellite orbits, station positions, EOPs and signal specific parameters can be estimated.

-GPS/GLONASS: space-based navigation system that provides location and time information in all weather conditions. Positioning is achieved triangulating the signals broadcast from, at least, four satellites of the constellation and simultaneously collected at a ground receiver. Despite the geometrical problem would require observations from three satellites only, a fourth signal must be acquired in order to synchronize the clocks on-board different spacecrafts among them and with the one installed at the ground station. Requiring the cheapest terrestrial equipment and considered the amount of scientific and civil applications of its data, the GNSS network is currently the densest in the framework of ITRF contributing techniques. As a consequence, it plays a major role in the distribution of ITRF products to users and in providing tying information between the different networks.

-DORIS: radio technique relying on the signal emitted by ground beacons and received on-board satellites. The concept is somewhat equivalent to the GNSS one, but the configuration is reversed. Since all measurements are recorded on-board the satellites, the only requirement for the installation of a terrestrial beacon is the availability of electrical power supply, facilitating the set-up of new stations even in remote areas. Despite encompassing just about 60 sites, the DORIS network is then remarkably distributed what makes it suitable for supporting orbit determination and ITRF dissemination.

-VLBI: system based on a network of telescopes observing astronomical radiosources, such as QUASARs (QUASi stellAR objects). The time lag between the arrival of the same signal at independent sites can be used to determine their relative distance (baseline). Due to the differential nature of the basic observable, the VLBI technique is insensitive to translation in the origin of the terrestrial system and cannot, therefore, contribute in the determination of the frame origin. On the other hand, VLBI is the only method, among those contributing to the ITRF computation, which does not rely on satellite measurements. As a consequence, it is the only system which can provide absolute information about the difference UT1-UTC: for satellite techniques, in fact, variations of such quantity could be completely compensated by a consistent rotation of the whole constellation and therefore absorbed in the estimation of orbit parameters.

In this study, TRFs realized with GPS/GLONASS and SLR measurements are computed and analyzed in order to assess the relevance of the contribution of space ties on-board GNSS satellites. A brief introduction on these techniques is provided in the following sections. The discussion is focused on the characterization of the system-specific observation equation and highlights the critical issues impacting the accuracy of the derived TRF parameters. A detailed description of the functional principles is, however, well established in the literature and therefore out of the scope of this introduction.

In addition, the chapter presents the datasets considered in this project and the strategies designed for the analysis of raw observations. The derived solutions are then validated against the official products provided by IAG Services for the computation of ITRF2014. Technique-specific long term frames are also computed and a comparison with ITRF2008 is carried out.

Results presented in this chapter are intended to confront the following questions:

1. Does the analysis strategy devised for this study provide results meeting the state of the art precision and accuracy of IAG Services' official solutions?

2. Does the analysis of data collected over the period 2011-2014, i.e. after the release of ITRF2008, reveal unprecedented features in the realization of frame parameters?

Least-Square Adjustment theory

Each technique contributing to the ITRF realization relies on a specific observable which depends on the status of a series of quantities. Therefore, variations in the measured observations can be used to infer changes in the independent variables which are of interest for the user. The mathematical expression detailing such dependency is called observation equation. This is in general a non-linear relation f of the form

f (x) = ô (2.1)
where the vector of parameters x lists all the relevant independent variables and ô contains the recorded observations. x includes station positions, satellite orbit parameters (except for VLBI for which quasar coordinates are pertinent), EOPs, delays associated with the signal propagation in atmosphere and system-specific quantities which can absorb unmodeled effects and technical biases.

Recorded observations ô are representative of the actual state o of the observable just up to the effects of measurement and modeling errors:

ô = o + ν (2.2)
where ν is called vector of residuals. If a set of a priori values x 0 is available for the unknown parameters x, so that

x = x 0 + δx (2.3)
with small corrective terms δx, Eq. 2.1 can be linearized with a first-order Taylor series expansion:

o + ν = f (x 0 ) + ∂f ∂x x=x 0 • δx (2.4)
In order to simplify the notation, the reduced observation vector l and the design matrix A can be introduced

l = o -f (x 0 ) (2.5) A = ∂f ∂x x=x 0 (2.6)
leading to the equation of residuals:

ν = A • δx -l (2.7)
Under the assumption that no systematic deviations affect the observations in the framework of the selected parametrization, ν contains only random errors normally distributed around 0. It is therefore reasonable to proceed in the resolution of eq. 2.1 by searching for the set of x values minimizing the norm of the vector of residuals. This approach is known as Least Square Adjustment. The space of the residuals can be provided with the metrics defined by the weight matrix P built on the basis of observation noise. The problem at hand consists, then, in finding

x : ν T Pν -→ min (2.8)
Substituting Eq. 2.7 in Eq. 2.8 leads to:

A T PA N •δx = A T Pl b (2.9)
where N is the normal matrix derived from the observation equation and b is typically referred to as the vector of the right-hand side.

Stacking of normal equations systems

If some parameters appear in different normal equation systems, a unique estimation can be provided exploiting all the available information. This merging procedure is called stacking and, provided that no correlations exist between the different groups of contributing observations, can be easily realized starting from the original independent systems. In fact, if two uncorrelated sets of observations l 1 and l 2 associated to the weight matrices P 1 and P 2 are available for a common set of parameters x, Eq. 2.9 prescribes that the estimation of x depends on the design matrix A and on the weight matrix P given respectively by

A = A 1 A 2 and P = P 1 0 0 P 2 .
(2.10)

A least square adjustment of the system so achieved leads to:

A T 1 P 1 A 1 + A T 2 P 2 A 2 • x = A T 1 P 1 l 1 + A T 2 P 2 l 2 .
(2.11)

On the other hand, if the two autonomous datasets are handled separately, the two estimates of x would result from:

A T 1 P 1 A 1 N 1 •x = A T 1 P 1 l 1 b 1 A T 2 P 2 A 2 N 2 •x = A T 2 P 2 l 2 b 2 .
(2.12)

The comparison of Eq. 2.11 and 2.12 shows that the stacking of two normal equations systems derived from uncorrelated observations simply results in a new normal equation system characterized by a normal matrix and a right-hand side vector provided by the sum of the original corresponding quantities:

(N 1 + N 2 ) • x = b 1 + b 2 (2.13)
In the outlined presentation it has been assumed that the different groups of observations depended exactly on the same set of parameters. In the most general case, however, such an assumption is not justified and Eq. 2.11 cannot be consistently written. In order to enable the computation of the proposed stacking, a new parameter vector should be introduced encompassing both common and systemspecific variables and each normal equation system shall be expanded accordingly. For each contributor, this is achieved by adding lines and columns of zeros in the original normal matrix and zero values in the original right-hand side vector in correspondence of the parameters peculiar to the other systems [START_REF] Seitz | Comparison of different combination strategies applied for the computation of terrestrial reference frames and geodetic parameter series[END_REF].

As a final comment, it should be noted that when stacking different normal equation systems, it might be desirable to differently weight the distinct contributions. This can be achieved by scaling the corresponding normal matrix and right-hand side vector of a factor s. If a scaling is applied, the stacked system would read

(sN 1 + N 2 ) • x = sb 1 + b 2 (2.14)
The presented combination scheme has been exploited in this study for the computation of GNSS+SLRtoGNSS weekly solutions and for the derivation of fully combined quasi-instantaneous frames which will be presented and detailed in Chapter 3. As all the technique-specific normal equation matrices have been generated with the same software, the highest possible consistency in data handling is guaranteed. This approach, in fact, ensures that the same modeling and parametrization are exploited for all the available observations [START_REF] Gambis | Why Combining at the Observation Level? In Reference Frames for Applications in Geosciences[END_REF].

The illustrated stacking strategy is also known in the literature as combination at the normal equation level. This approach is typically opposed to the combinations that could be carried out at other stages of the least square adjustment, namely at the level of the observation equations or at the level of parameters (solutions). Theoretically speaking, differences do exist between these possibilities. The combination at the observation level is the only one that ensures a consistent screening of the available data as all the measurements are composed into one system and jointly pre-processed. On the other hand, the combination at the parameter level differs from the other approaches because, if not properly handled, constraints applied to derive the contributing solutions are reflected in the combination product. For practical applications, however, these three alternatives are actually equivalent. The effects of performing data pre-processing technique-wise or in a common adjustment are believed to be "unverifiable, particularly in the case of the ITRF computation" [START_REF] Seitz | Comparison of different combination strategies applied for the computation of terrestrial reference frames and geodetic parameter series[END_REF]. Moreover, SINEX (Solution (Software/technique) INdependent EXchange) files, which are currently the standard format for solutions presented by IAG Services, can accommodate the normal matrix of the constraints exploited in the derivation of the solution, so that the original unconstrained information can be easily recovered according to 1.22. In addition, if the theoretical singularities of the unconstrained normal equation matrix are known, even unreported minimal constraints can be removed [Rebischung, 2014, App. B.2.4]. As already stated, then, the three approaches can be considered as equivalent for all practical purposes.

GNSS

The GNSS technique encompasses nowadays different constellations operating on the same basic positioning principle. The first deployed project has been the American NAVSTAR GPS (NAVigation Satellite Timing And Ranging Global Positioning System), developed since 1973 by the US Department of Defence for military use and later made available to support civil applications. Russia, China, Japan and the European Union have subsequently developed their own systems, namely the GLONASS, Beidou, QZSS and GALILEO constellations respectively, which have attained, or will attain in the close future, the operational phase. Among the main differences between the diverse systems, it should be noted that GPS, GLONASS and GALILEO provide global coverage, while the trackability of Beidou and QZSS satellites is/will be limited to regional areas.

For the ITRF2014 computation, most IGS ACs provided solutions based exclusively on GPS data, with the only exceptions of CODE, ESA and GRGS-CNES 1 [START_REF] Rebischung | The IGS contribution to ITRF[END_REF] which included also GLONASS observations. Considering that all GLONASS satellites are operationally equipped with laser retroreflectors arrays (LRA) against just two GPS vehicles, namely GPS-35 and GPS-36, data from both constellations are included in this study. The analysis strategy devised in this work to properly handle GLONASS observations is presented in Sec. 2.2.2. As none of the other GNSS systems is relevant in the context of this thesis, the following discussion will focus just on GPS and GLONASS. The acronym GNSS will be used to jointly address both of them and them only, notwithstanding the fact that, properly speaking, also the other constellations should be included.

The GPS constellation is fully deployed since 1994 and encompasses 24 satellites 1 CODE Center for Orbit Determination in Europe ESA European Space Agency GRGS-CNES Groupe de Recherches de Géodésie Spatiale -Centre National d' Études Spatiales distributed over six orbital planes inclined at 55➦ and separated by 60➦ in the equatorial plane. The selected orbits are quasi-circular and characterized by a semi-major axis of about 26600 km and a period of 11 hours and 58 minutes, i.e. about half a sidereal day. Such configuration ensures the continuous coverage of the whole globe, as it guarantees that at least four satellites can be tracked by any point at any time. All the satellites of the constellation transmit on the same two carrier frequencies (a third one is additionally available, but only for last generation spacecrafts) and can be distinguished by ground stations thanks to a message modulated upon them (Code Division Multiple Access, CDMA, approach).

The GLONASS project has been characterized by two operational phases: the first one attained in 1995, followed by a severe constellation degradation due to lack of funds, and a more recent one began in 2003 and culminated in 2011 when the nominal extent of the orbiting fleet was restored. The encompassed 24 satellites are evenly assigned to three orbital planes inclined at 64.8➦ and separated of 120➦ in the equatorial plane. As for the GPS case, quasi-circular orbits have been chosen, with a semi-major axis of about 25510 km and a period of 11 hours and 15 minutes. Moreover, also GLONASS satellites broadcast their signal on two frequencies, but these are shifted from vehicle to vehicle (Frequency Division Multiple Access, FDMA, approach). This peculiarity is exploited by receivers to identify the emitting spacecraft, but poses some critical issues on the data analysis (see Sec. 2.2.2). A more comprehensive review of GPS and GLONASS systems can be found in Dach et al. [2007, Section 2.1].

Principle

GNSS systems rely on the transmission of an electromagnetic signal between GNSS satellites and active receivers. For the purposes of this study, only ground receivers are here taken into account, but a number of space missions (e.g. Jason 1-2-3, ICESat, GRACE, TerraSAR-X...) are equipped with GNSS positioning payloads in order to support precise orbit determination (POD). For GPS and GLONASS constellations, each satellite transmits on two carriers, L1 and L2, which are integer multiples of the main frequency provided by the accurate atomic clocks carried onboard GNSS satellites (at least two caesium and as many as two rubidium clocks in the case of GPS). The broadcast signal is then additionally modulated according to specific codes also derived from the fundamental frequency, Fig. 2.1.

Receivers are capable to measure the phase of the incoming signal and are provided with an oscillator on the basis of which the broadcast information can be reproduced. Assuming that all the involved clocks are perfectly synchronized, the phase of the wave generated within the receiver at the time of recording t differs from the acquired phase by an amount which depends on the signal time of flight τ . Considering the finite velocity c of electromagnetic waves, in fact, measurements In the case of GPS two codes are employed. The C/A (clear acquisition) code has a frequency of 1/10 of the fundamental one maintained by the clocks and is modulated on L1 only; the P (precise) code, instead, provides a ten times higher accuracy and is superimposed to both L1 and L2. The navigation message, identifying the satellite that broadcast the signal, its health status and other operative information, is also modulated on both carriers. carried out at time t actually sense the signal emitted by the satellite at time t -τ . The basic observable ψ i F k (t), then, consists in the phase difference [START_REF] Dach | Manual of Bernese GPS Software Version 5.0[END_REF]:

ψ i F k (t) = φ F k (t) -φ i F (t -τ ) + n i F k (2.15)
where In operational conditions, Eq. 2.15 can be linearized with a first order Taylor series expansion

ψ i F k (t)
ψ i F k (t) = φ F k (t) -φ i F (t) + τ f F + n i F k . (2.16)
For ideal, perfectly synchronized clocks, the difference φ F k (t) -φ i F (t) is equal to zero, but in operational conditions station-and satellite-clock errors, δ k and δ i respectively, must be taken into account so that

φ F k (t) -φ i F (t) = δ k -δ i f F .
(2.17)

If the relation expressed in Eq. 2.17 is substituted in Eq. 2.15 and this last is multiplied by the carrier wavelength λ F , the basic form of the phase observation equation in meters is obtained

L i F = ρ i k + cδ k -cδ i + λ F n i F k (2.18)
where ρ i k = cτ represents the geometric distance between satellite i at time t -τ and the receiver k at time t. For users interested in positioning applications, ρ i k can be explicitly related to the 3-dimensional coordinates of the involved satellite-receiver pair through

ρ i k = R EOP • X k (t) -X i (t -τ ) (2.19)
where X k is the coordinate vector of receiver k at the time of signal acquisition t expressed in a Terrestrial Reference Frame X i is the coordinate vector locating the satellite i at the time of signal emission t -τ in a Geocentric Celestial Reference Frame R EOP is the rotation matrix of Earth Orientation parameters necessary to transform coordinates from a Terrestrial to a Geocentric Celestial Reference Frame, see Sec. 1.3.1.1. The application of such operator allows expressing station and satellite positions in the same system, so that it is possible to differentiate the two vectors.

The electromagnetic signal broadcast by the satellite, however, interacts with the environment through which it propagates causing delays that need to be accounted for in the observation equation in order to achieve the best possible accuracy. In addition, systematic errors do also affect the measurements. Equation 2.15, then, should be expanded to

L i k = R EOP • X k (t) -X i (t -τ ) + + λ • N i k + c • δt k -c • δt i + + δρ trop + δρ ion + δρ phas + δρ rel + δρ mult + ǫ i k (2.20)
where the added corrective terms accounts for -ρ trop delays accumulated for the signal propagation through troposphere -ρ ion delays accumulated for the signal propagation through ionosphere -ρ phas corrections for instabilities/inaccurate characterization of phase center offsets and variations at satellite and receivers antennas -ρ rel corrections for relativistic effects affecting satellite clocks and orbits as well as signal propagation (Shapiro effect) -ρ mult multipath effects

-ǫ i k measurement errors
An observation of this type can be formulated for each of the two fundamental carriers and for any of their linear compositions. In GNSS analysis practice, several combinations are often employed according to the different parameters which need to be addressed. Solutions computed for this study, are obtained with the quasi ionosphere-free linear combination

L 3 = 1 f 2 1 -f 2 2 (f 2 1 L 1 -f 2 2 L 2 )
for which the observation equation is, in first approximation, independent from ionospheric path delays. Apart from second and higher order terms, in fact, ionospheric refraction is given by

I i k ≈ 1 f 2 (2.21)
where f is the carrier frequency, and its relevant contributions at f 1 and f 2 cancel out in the L 3 equation.

Additionally, differences of the original observation equation can be exploited in order to eliminate or, at least, reduce some biases. Results presented in this study derive from a double-difference adjustment between the data collected by pairs of receivers, kl, simultaneously tracking pairs of satellites, ij:

L ij F kl = L i F kl -L j F kl (2.22)
where

L i F kl = L i F k -L i F l
(2.23) Equation 2.22 is independent from satellite and receiver clock errors and is therefore often used in GNSS data analysis.

Analysis Strategy

Results presented in this study are obtained from GPS and GLONASS observations processed in double-difference mode exploiting the quasi ionosphere-free linear combination of L1 and L2. Conventional movements have been modeled according to IERS Conventions [Petit and Luzum, 2010, Chapter 7]. Raw GNSS observations have been treated with the Bernese GNSS Software v. 5.2 [START_REF] Dach | Manual of Bernese GPS Software Version 5.0[END_REF]. A dedicated automated routine (Bernese Processing Engine, BPE) has been developed to compute daily normal matrices from the records of 104 globally distributed stations, see Fig. 2.4. In the analysis, data sampling rate is set to 180 s, the elevation angle cutoff is fixed to 3➦ and an elevation-dependent weighting in form of 1/cos2 z is imposed.

Following the adjustment approach presented in Sec. 2.1, accurate a priori values must be provided for the estimated parameters. All the stations selected for the analysis belong to the IGS permanent network, a priori coordinates have then been extrapolated from IGS cumulative solution 2 aligned with ITRF2008. Concerning orbit parameters, CODE products have been preferred to official IGS ones because the firsts are computed from a simultaneous and fully consistent adjustment of both GPS and GLONASS data [START_REF] Dach | GNSS processing at CODE: status report[END_REF]. In addition, it should be considered that precise satellite orbits are distributed as a series of spacecraft positions tabulated every 15 minutes. These must be interpolated in order to retrieve the instantaneous position of the satellite at any time. In the Bernese Software, developed and used at CODE for operational GNSS analysis, the interpolating orbit is the solution of a physical model of the forces acting on the satellite, as will be briefly introduced later in this section. The resulting orbital product can therefore be accurately reproduced by other Bernese users, provided that the same parametrization is exploited. IGS combined orbits, on the other hand, benefits from the contribution of all the official IGS ACs, but do not follow any motion equation and can hardly be exactly reproduced during the interpolation. Even if effects are likely to be within the orbital uncertainties, CODE products are preferred in this work. A priori values for polar motion and length of the day parameters, i.e. for EOPs, are provided by the IERS C04 series aligned to ITRF2008 [START_REF] Bizouard | The combined solution C04 for earth orientation parameters consistent with international terrestrial reference frame[END_REF]. Official IGS values based on ITRF2008 have been adopted as a priori for Satellite Antenna Offsets (SAOs) and Satellite Antenna Phase Center Variations (PCVs). As far as available, absolute elevation-and azimuth-dependent PCVs and L1/L2 offsets have been used for ground antennas. Antenna radome calibrations and GLONASS-specific models have also been accounted for.

The pre-processing phase of the analysis encompasses several different steps:

-clock syncronization: receiver clocks need to be synchronized to GPS time.

Even if receiver clock errors δ k do not explicitly appear in the double-differentiated observation equation, they still impact the determination of the geometric term. Synchronization errors are, in fact, proportional to the uncertainty in the radial positioning of the satellite with reference to the receiver, as shown by dρ i k = -ρi k dδ k which is the time derivative of the relation determining the impact of δ k over the radial distance between satellite i and receiver k

ρ i k (t) = ρ i k (t k -δ k )
.

-baseline generation: the pre-processing further encompasses the definition of the set of receiver pairs (baselines) that will constitute the basis for the formation of double differences. If N receivers are encompassed in the tracking network, only N -1 independent baselines can be formed. The selection criterion adopted in this work maximizes the number of simultaneous observations from the associated stations to the same satellite. In other words, it is identified the set of ground baselines from which the maximum number of double-differences can be successively formed.

-cycle-slip detection and handling: each selected baseline is checked for cycle-slips, i.e. for unexpected signal losses while tracking a satellite. These events have to be addressed and either corrected or absorbed in additional parameters set up in the observation equation of the affected data. In fact, the unknown number of integer cycles n i k does not change as long as a continuous receiver-satellite lock is ensured. The same parameter can then be set up in the observation equations referring to successive epochs, reducing the number of variables to estimate. Physical obstructions (trees, buildings...), low signalto-noise ratio or equipment failures, however, can cut the lock off and, if it is impossible to estimate the number of cycles lost during the interruption interval, a new ambiguity parameter needs to be set up.

-outlier screening: raw data are screened for outliers in an iterative process based on the comparison of observation residuals with progressively smaller thresholds. Residuals are obtained as the difference between actual records and the corresponding quantities computed from the results of a preliminary solution encompassing only station coordinates, constrained at 1 cm level in each component, and tropospheric parameters. All the other unknowns are fixed to their a priori values and data are down-sampled to three minutes to reduce computational time.

-ambiguity resolution: before computing the final solution, a baseline-bybaseline estimation of the initial ambiguities n i kl is attempted. Strictly speaking, this step is not part of the pre-processing, as n i kl are actual parameters of the observation equation. However, considering their high number (at least one per baseline-satellite pass), the computational burden associated with their estimation in the final session-wise computation is just prohibitive. For this reason, the different algorithms for ambiguity resolution implemented in the Bernese Software are applied in baseline mode and the estimated values are introduced as known quantities in the subsequent processing. The diverse resolution strategies are based on specific assumptions that can be verified just within a limited range of baseline length. According to the CODE Analysis Strategy Summary for IGS repro2 campaign, [START_REF] Code | CODE Analysis Strategy Summary for IGS repro2 campaign[END_REF], the choices adopted in this study are reported in Tab. 2.1. The interested reader is addressed to the Bernese Software manual [START_REF] Dach | Manual of Bernese GPS Software Version 5.0[END_REF], Chapter 8 and references therein] for a detailed description of the different algorithms.

BL length Strategy GLONASS?

< 20 km Direct L1/L2 method yes < 200 km Phase-based widelane/narrowlane method yes < 2000 km Quasi-Ionosphere-Free (QIF) approach yes, same f < 6000 km Melburne-Wubbena approach no Table 2.1 -Ambiguity resolution strategies adopted according to baseline length. Applicability to GLONASS system is also indicated.

Ambiguity resolution is particularly troublesome for GLONASS systems which rely on the FDMA principle for the identification of the transmitting satellite. Different spacecrafts broadcast over slightly separate frequencies (apart from antipodal vehicles which share the same channel) so that the double-difference observation equation reads

L ij kl = ρ ij kl + nij kl λ i + ∆λ ij n j kl + other biases (2.24)
where λ i is the carrier wavelength of satellite i and ∆λ ij is the carrier wavelength difference between satellites j and i. Because of the presence of ∆λ ij , resolution strategies exploiting the integer nature of double-difference ambiguities cannot be applied for baselines encompassing GLONASS receivers.

The percentages of ambiguities resolved with each strategy for both GPS and GLONASS observations are displayed in App. A. About 80% of GPS initial phase ambiguities have been fixed, while a satisfactory solution was found in just 40% of the cases when GLONASS baselines were screened. Unresolved ambiguities were kept in the normal equation system and estimated within the final adjustment.

Once the data pre-processing is completed, the normal equation matrix associated to the observation system in double-difference mode can be derived. A suitable parametrization must be introduced for the different quantities:

-station coordinates are set up as constant offsets over the time span covered by the observations at hand. Mean position are considered to be representative of the actual station position up to about 10 days; for longer periods, at least linear drifts should be taken into account. In the analysis, non-linear conventional displacements are modeled according to IERS Conventions [Petit and Luzum, 2010, Chapter 7], so that the actual unknowns set up in the observation equations are the regularized station positions.

-satellite positions are computed as the solution of the differential motion equation:

       Ẍ = a i t, Ẋi (t), X i (t), p i 1 , ..., p i m X i (t 0 ) = X i 0 Ẋi (t 0 ) = Ẋi 0 (2.25)
where a i is the function describing the acceleration fields acting on the satellite i. On the one hand, such term includes well modeled conventional contributions, as gravitational acceleration. In addition, insufficiently well modeled effects, such as those caused by solar radiation pressure, must also be estimated. In the Bernese Software framework, the complete characterization of each orbital arc is then achieved through the identification of the initial six obsculating Keplerian elements and a series of parameters p i 1 , ..., p i m accommodating for empirical accelerations. In this study, the lasts have been parametrized exploiting the Extended CODE Orbit Model (ECOM; [START_REF] Springer | A new solar radiation pressure model for GPS satellites[END_REF][START_REF] Beutler | Extended orbit modeling techniques at the CODE processing center of the international GPS service for geodynamics (IGS): theory and initial results[END_REF]), which, in its complete formulation, encompasses a constant offset and two once-per-revolution periodic accelerations for each of the three axes D, Y and B of a satellite-Sun oriented frame, see Fig. The adopted model for empirical accelerations is the one exploited by a large number of official ACs within the framework of the repro2 campaign which constituted the IGS contribution to the realization of IRTF2014. Nevertheless, recent studies [START_REF] Rodriguez-Solano | Reducing the draconitic errors in GNSS geodetic products[END_REF][START_REF] Meindl | Geocenter coordinates estimated from GNSS data as viewed by perturbation theory[END_REF] have shown that the ECOModel was responsible for spurious draconitic lines in geophysical parameters derived from combined GPS/GLONASS analysis, such as Earth Orientation Parameters (ERPs) and geocenter estimates. These effects have been attributed to the lack of twice-per-revolution signals in the D-direction, see the animations in [START_REF] Dach | Updating the CODE GNSS Orbit Model[END_REF] for a direct representation. Signals are more evident for GLONASS satellites since the cylindrical shape of these spacecrafts is characterized by a strong variation of the structural portions exposed to direct solar irradiation during an orbital cycle. As a consequence of these findings, the traditional ECOModel has been revisited and an extended version, encompassing twice-per-revolution terms, is now available [START_REF] Arnold | CODEś new solar radiation pressure model for GNSS orbit determination[END_REF]. Alternatively to the empirical approach developed at CODE, attempts are being made to actually model the solar radiation pressure impacting GPS and GLONASS satellites. In these studies, the box-wing approach for the characterization of the spacecraft surfaces and their interaction with solar radiation is typically followed [START_REF] Rodriguez-Solano | Reducing the draconitic errors in GNSS geodetic products[END_REF]ESA, 2014].

-delays caused by the signal propagation through the atmosphere are composed of an ionospheric and a tropospheric contribution. In first approximation, ionospheric effects can be neglected when running an analysis exploiting the L3 linear combination of the carrier frequencies. The so-called higher order ionospheric terms, however, are known to be responsible for spurious variations up to 1 cm in station height component during period of maximal solar activity and network distortions may also be observed [START_REF] Petrie | Higher-order ionospheric effects on the GPS reference frame and velocities[END_REF]; Hernandez-Pajares et al. [2007]. CODE's estimate for second and third order terms have been applied in this study.

Concerning the troposphere, IERS Conventions [Petit and Luzum, 2010, Chapt. 9] recommend to parametrize the relevant correcting terms accounting for zenith delays and North/South, East/West gradients. These lasts are typically estimated as daily linear functions; zenith delays, instead, can be inferred from numerical weather models, except for the component depending on the water vapor content that needs to be evaluated during the adjustment. With state of the art products, in fact, troposphere humidity cannot be provided with sufficient accuracy. The prescribed parametrization for troposphere delays impacting the signal broadcast from satellite i and received at receiver k reads:

T i k = m h (e)D hz + m w (e)D wz + m g (e) (G N cos(a) + G E sin(a)) (2.27)
where -e and a are the elevation and the azimuth of the satellite with reference to the receiver -D hz and D wz are, respectively, the zenith dry and wet delays; the latter has been parametrized as a piece-wise linear function sampled every two hours. The multiplicative terms m h and m w are the hydrostatic and wet mapping functions which allow to project zenith delays along the satellite line-of-sight. The Vienna Mapping Function (VMF1) [START_REF] Böhm | Global Mapping Function (GMF): A new empirical mapping function based on numerical weather model data[END_REF] has been used in this analysis.

-G N and G E are the N/S and E/W tropospheric gradients; the relative mapping function, m g , is provided by [START_REF] Chen | Effects of atmospheric azimuthal asymmetry on the analysis of space geodetic data[END_REF] -EOPs are set up daily as offsets plus drifts. No constraints are applied with the only exception of U T 1 -U T C which is fixed to the a priori value given its complete correlation with orbital parameters.

-satellite phase center offsets have been set up, but tightly constrained to the official conventional values provided by the IGS. It is important to note that these parameters are critical for the realization of the GNSS network scale. [START_REF] Zhu | Detecting discontinuities in GNSS coordinate time series with STARS: case study, the Bologna and Medicina GPS sites[END_REF] showed that a mean error δz in the determination of z -P COs would induce a mean scale change δs according to δs = 7.8 • δz, where δz and δs are expressed in m and ppb respectively. Accurate pre-launch calibrations, however, are not available for currently flying satellites and their estimation during the analysis is impaired by the almost perfect correlation of these parameters with satellite clock offsets, zenith wet delays and station heights [Rebischung, 2014, Fig. 4.3]. As a consequence, reliable z-P CO values are unknown and the GNSS technique cannot provide a significant contribution to the realization of the terrestrial scale. Fixing z -P COs to the conventional IGS values implies forcing the GNSS scale to rely on the ITRF2008 one.

Dataset and result validation

The dataset selected for this study encompasses the observations collected by 104 GPS/GLONASS permanent IGS stations over the period 2011-2014. The starting epoch has been chosen considering that the flying GLONASS constellation has been completely restored in 2011. According to [START_REF] Blewitt | Effect of annual signals on geodetic velocity[END_REF], four years of data should be sufficient for a reliable determination of station positions and velocities. All GNSS stations co-located with active SLR instruments are taken into account; then, a series of globally distributed sites is identified to complete the network. Besides geographical location, selection criteria privileged stations which acquired continuously over the whole time period, in order to ensure the highest possible network stability. GLONASS capable stations have been preferred, when available. The map of the selected network is reported in Fig. 2.4; SLR stations are also reported to highlight the distribution of available ground co-locations. The map clearly shows that, even in most recent years, the spatial distribution of the SLR network is far from being optimal. On the bright side, almost all sites hosting an SLR telescope are provided with a GPS/GNSS antenna. Since the large majority of SLR sites is installed in the Northern hemisphere and all co-located GNSS instruments are included in the network, also the distribution of the latters results unbalanced.

During the four years under study, 22 of the selected stations have been upgraded to GLONASS capable systems, as shown in Fig. 2.5 which displays the GLONASS tracking network at the beginning of 2011 and at the end of 2014 respectively. At this epoch, almost the 80% of the network participated in GLONASS tracking.

The different comparisons and combinations carried out for this study, see Chapter 3, take advantage of a homogeneous sampling of all the involved solutions. Technique-specific services, however, issue different recommendations for the analysis of raw data. Since GPSweek 1702, IGS is providing daily frame solutions in order to avoid any loss in GNSS sensitivity to non-tidal loadings . ILRS, on the other hand, recommends to process SLR observations to LAGEOS in weekly intervals. Therefore, in the analysis carried out for this study, GNSS solutions have been generated on a daily basis and then downsampled to SLR-comparable intervals stacking successive normal equation systems. The implemented approach is similar ments the recent findings by [START_REF] Lutz | Impact of the arc length on GNSS analysis results[END_REF], showing that long-arc estimates are "superior" to daily ones in terms of WRMS of associated EOPs and geocenter coordinates, especially when EOPs are made continuous at the boundaries of the middle day. The latter feature, however, was deliberately not implemented in the study at hand in order to minimize the discrepancies with the purely one-day approach recommended by the IGS. The exploited strategy is then comparable to the Cx case analyzed in [START_REF] Lutz | Impact of the arc length on GNSS analysis results[END_REF], which is proved to produce results of intermediate quality between 1-day and 3-day solutions with continuity conditions imposed.

The GNSS time series is then down-sampled stacking three consecutive nonoverlapping 3-days normal equation systems centered on the ILRS week. No long arc computation is performed in this step. The parametrization of the final nine-days solution, then, encompasses a single set of station coordinates and three independent orbital arcs per observed satellite; all other quantities keep the parametrization implemented in the original daily solution. For each nine-day normal equation matrix, the number of processed observations and of estimated parameters is presented in Fig. 2.7. Other statistics regarding the GNSS processing can be found in appendix A. A long term stacking of the computed time series has then been performed in order to fully characterize their information content. Figure 2.8 shows the WRMS of station position residuals between each input nine-day solution and the frame resulting from the long term stacking. Two different accumulations have been performed encompassing or not the seven transformation parameters; in the first case IC are imposed on origin and scale offsets, while MC are applied to rotations. The average values of the WRMS are reported in Tab. 2.2. Typical values of a few millimeters for horizontal components and about 5 mm for the vertical are observed when transformation parameters are estimated, while significantly higher numbers are found in the second case. When realizing long term secular frames, non-linear station displacements induced, for example, by spurious origin oscillations, geocenter motion and loadings effects leak to the transformation parameters [START_REF] Collilieux | Effect of the satellite laser ranging network distribution on geocenter motion estimation[END_REF][START_REF] Collilieux | Impact of loading effects on determination of the international terrestrial reference frame[END_REF]. Their exclusion from the stacking, then, necessarily produces an increment in station residuals. Nevertheless, if the two estimated long-term frames are compared, the observed discrepancies are negligible, see Tab. 2.3, meaning that non-linear displacements do not impact the estimation of frame parameters. Such a result was already shown by Altamimi et al. [2013a] for SLR results derived with data collected over the period 1993 -2009. It has been shown, however, that four years of data are already sufficient to mitigate the effect of periodic signals, in accordance with the findings of [START_REF] Blewitt | Effect of annual signals on geodetic velocity[END_REF].

E [mm] N[mm] U [mm]

case a 1.9 2.1 4.6 case b 4.5 10.0 10.3 Table 2.2 -Average values of the WRMS of nine-day GNSS station position residuals estimated during the long term stacking of the time series, see Fig. 2.8. In case a) transformation parameters have been set-up during the accumulation, while they were ignored in case b).

T x [mm] T y [mm] T z [mm] Scale [ppb]
0.0 ± 0.0 0.1 ± 0.0 0.0 ± 0.0 -0.03 ± 0.01 Table 2.3 -Translation and scale offsets estimated between the secular frames computed from GNSS 9-day solutions either setting up or ignoring transformation parameters during the long term stacking.

The frame information expressed by the last four years of GNSS data have also been compared to ITRF2008. For this purpose, another long term stacking has been performed setting up seven transformation parameters for each input solution and applying minimum constraints with reference to ITRF2008. The estimated offsets for origin and scale are reported in Fig. 2.9. For each 9-day solution, the origin of the network coincides with the quasiinstantaneous GNSS materialization of the Earth's Center of Mass (CoM), as the equations of motion governing the orbit evolution hold only in a frame centered on the barycenter of the physical system. For crustal-based frame, however, the actual realization of the origin is inherently provided by the coordinates of a set of points distributed over the Earth's surface which is known to move, with reference to the total Earth's Center of Mass (CM), in response to crustal deformation and fluid transportation in and between atmosphere, oceans, surface aquifers and ice-sheets [START_REF] Eanes | Low-frequency geocenter motion determined from SLR tracking[END_REF][START_REF] Wu | Geocenter motion and its geodetic and geophysical implications[END_REF]. The relative displacement of the geometric center of the solid Earth's surface (Center of Figure ,CF) with reference to CM is known as geocenter motion [START_REF] Eanes | Low-frequency geocenter motion determined from SLR tracking[END_REF][START_REF] Wu | Geocenter motion and its geodetic and geophysical implications[END_REF] and it is expected to show both linear and non-linear variations over time [START_REF] Dong | Geocenter variations caused by atmosphere, ocean and surface ground water[END_REF][START_REF] Chen | Geophysical interpretation of observed geocenter variations[END_REF].

As geocenter motion is not directly measurable, a characterization of the phenomenon can only be inferred thanks to collateral observations, a thorough review is provided in Sec. 1.3 of [START_REF] Rebischung | Can GNSS Contribute to Improving the ITRF Definition?[END_REF]. The analysis of quasi-instantaneous translation parameters estimated with reference to a secular frame, for example, can be employed to investigate the non-linear part of the signal, network shift approach. This strategy has been largely used in the literature, even if [START_REF] Collilieux | Effect of the satellite laser ranging network distribution on geocenter motion estimation[END_REF] pointed out that it is intrinsically biased by the non-uniform distribution of geodetic stations over the Earth's surface. Translation offsets are actually computed, in fact, for the Center of Network (CN), which, in consequence of local displacements, may be characterized by a relative motion with reference to CF. Temporal variations of the CN-CF vector depend on the geometry of the available tracking stations and the related artificial contribution to geocenter motion is then known as network effect. The greatest impact is expected for highly inhomogeneous configurations since nonconventional movements induced by geophysical loadings are spatially correlated. In addition, a poor distribution may generate correlations between the different transformation parameters estimated during the stacking, causing spurious signals to leak into the "geocenter motion" estimation.

Besides the intrinsic limitations of the approach, it should be considered that, depending on the technique-specific sensitivity to the geocenter, the estimated signal can reflect modeling deficiencies and observation errors rather than the actual geophysical compensation to mass displacements.

The stated shortcomings encourage some caution in the interpretation of the results shown in Fig. 2.9. Nevertheless, translational parameters show at least a qualitative agreement with theoretical expectations. A clear seasonal seasonal cycle is, in fact, identifiable for T x and T y in accordance with the experimental evidence that the dominant loadings act at annual frequency, in association to fluid exchanges between the Northern and Southern hemispheres [START_REF] Wu | Geocenter motion and its geodetic and geophysical implications[END_REF]. The signal in T z , instead, is characterized by abrupt excursions with a period of about 4 months. These oscillations likely result from technique-specific biases rather than actual geophysical signals. [START_REF] Meindl | Geocenter coordinates estimated from GNSS data as viewed by perturbation theory[END_REF], for example, found amplified analogous oscillations in a GLONASS only solution. In addition, a collinearity analysis carried out in Rebischung et al. [2014] shows that the z-component of the geocenter is extremly poorly determined by GNSS systems. In a standard analysis encompassing clock and troposphere parameters, in fact, all geocenter components correspond to quasi-singularities of the normal equation system, but the insensitivity to the z-component is further amplified by the estimation of the ECOM D 0 , B C and B S terms of the orbital model. Unresolved phase cycle ambiguities, then, may absorb part of the geocenter signal as they are set up as constant biases per pass of each satellite over each station. The higher impact is observed, again, on the determination of the z-component [Rebischung, 2014, Sec. 4.5.3]. The behavior expressed by the T z time series in Fig. 2.9, then, basically reflects modeling deficiencies and observation errors and is not significant for inferring any geophysical information.

In the analysis proposed in this study, PCOs are set up, but tightly constrained to a priori values, according to the practice adopted for the repro2 campaign. By so doing, GNSS network scale is conventionally inherited from the ITRF2008 realization [START_REF] Rebischung | IGS08: the IGS realization of ITRF[END_REF]. Non zero mean offsets in the scale panel of Fig. 2.9 are anyway expected as a consequence of modeling updates concerning Earth's radiation pressure and higher order ionospheric terms. In addition, six GPS and five GLONASS satellites have been launched since October 2012 and only preliminary PCOs were available for the processing [START_REF] Rebischung | The IGS contribution to ITRF[END_REF]. The observed annual cycle, on the other hand, can rather be interpreted as a network effect. The scale information, in fact, correlates, among other parameters, with the mean network height which is subjected to seasonal variations in consequence of non-modeled loading effects aliased in station coordinates. Loading displacements are spatially correlated and average to zero because of mass conservation, so the impact on the scale factor tends to vanish for perfectly homogeneously distributed network. As stations in Northern hemisphere clearly outnumber those distributed in the Southern, see Fig. 2.4, it is not surprising to detect a seasonal cycle in the proposed plot.

The mean values of the estimated transformation parameters shown in Fig. 2.9 are reported in Tab. 2.4. The analogous quantities computed in the framework of the repro2 campaign are also listed for validation purposes. Presented results are derived from the official IGS combined solution and from the contribution provided by the Center for Orbit Determination in Europe [START_REF] Rebischung | The IGS contribution to ITRF[END_REF]. The relevant values are obtained propagating the results of [Rebischung et al., 2016, Tab. 5 and 7] to the mid epoch of the interval analyzed in this study, 2013.0. External results are, then, computed from data collected throughout the time span exploited for the derivation of the ITRF2014, i.e. 1994 -2014. The estimates provided by CODE constitute a valuable term of reference for the analysis strategy developed in this study, as CODE's solution is computed with the Bernese Software and includes GLONASS data. The processed network, however, is about three times larger than the one exploited for this thesis; such a wealth of available data ensures that accurate estimates can be provided also relying on daily orbital arcs.

Tx Ty Tz Scale [mm] [mm] [mm] [ppb]

this study -5.0 -1.4 1.0 -0.14 IGSrepro2 -3.8 0.6 -5.6 -0.27 CODE -5.0 -0.8 -3.3 -0.14 Table 2.4 -Mean GNSS origin and scale offsets with reference to ITRF2008.

A remarkable agreement between estimates computed for this study and CODE's results is observed; differences in T x , T y and the scalef actor are, in fact, within 1 mm, and the slightly bigger discrepancy in T z is still acceptable considering the higher peak-to-peak variations and the spurious periodic signals shown by the relevant time series. The agreement with the IGS combined solution is a bit worse, but still within the variability expressed by the estimates of the different contributing ACs.

SLR

The Satellite Laser Ranging technique is based on the determination of the station-to-satellite range (i.e. radial distance) from the measurement of the roundtrip time of flight of an optical pulse shot and acquired by the same ground system after retro-reflection upon an orbiting target. The first SLR data ever acquired was obtained at the Goddard Space Flight Center in 1964, when a laser return from the Beacon Explorer 22B satellite was recorded. Five years later, the ranging capacity was extended from artificial satellites to the Moon, thanks to the panels of retroreflectors installed by Apollo 11 astronauts. In the following years, Apollo 14 and 15 and two unmanned Soviet Lunakhod missions, allowed deploying four additional arrays.

Several SLR-dedicated missions (LAGEOS, AJISAI, Starlette, Stella...) have been launched over the years in order to exploit the full potential of the technique concerning the derivation of geodetic parameters. SLR has, then, been able to contribute to the modeling of the Earth's gravity field, the measurement of global tectonic plate motion, the determination of Earth Orientation Parameters and, as previously stated, the realization of the ITRF origin and scale. In addition, SLR has been providing calibration for space-borne radar altimeters and validation for precise orbits determined for remote sensing and GNSS satellites4 .

In the context of ITRF realization, data of the two LAser GEOdynamics Satellites (LAGEOS) are the prominent contributors. Launched in 1976 and 1992 respectively, these missions rely on spacecrafts optimized for the estimation of geodetic parameters. The satellites consist of aluminum-covered brass spheres with a diameter of 60 cm and masses of 400 and 411 kg. The satellite bodies host 426 corner-cube retro-reflectors, see Fig. 2.10, homogeneously allocated over the entire surface. In order to ensure that the different sectors of the spacecraft have the same reflectivity, the optical properties of each corner-cube have been tested to support the development of an ad-hoc distribution scheme. The satellites are completely passive and have no on-board sensors or electronics. They orbit at an altitude of about 6000 km, neous both concerning spatial distribution, see Fig. 2.4, and tracking capabilities, given the great variety of technical specifications of operating systems. All the stations participating in the ILRS network are equipped with, at least, third generation lasers able to transmit well collimated pulses with high power density over large distances. Nevertheless, ranging characteristics such as single-shot and normal point precision or the signal-to-noise ratio are strongly station-dependent. As pointed out in [START_REF] Sośnica | Satellite laser ranging to GPS and GLONASS[END_REF], in fact, "almost every SLR station has a different technology". Diverse lasing media and tracking policies are currently exploited, resulting in a variety of possible ranging frequencies, laser pulse widths and repetition rates. In addition, tracking performances are impacted by the operative capabilities of the other fundamental elements of an SLR station, such as timer, photo detectors, filters and even the mounting system. A complete review of the infrastructure of the SLR ground segment is provided in [START_REF] Sośnica | Determination of precise satellite orbits and geodetic parameters using satellite laser ranging[END_REF]; a comparison of Hz and kHz station performances with reference to SLR tracking residuals is offered in [START_REF] Sośnica | Satellite laser ranging to GPS and GLONASS[END_REF].

Principle

The basic observable of the SLR technique is the two-way travel time of the optical signal shot and received by the same ground station after reflection upon targets equipped with suitable optical devices. As illustrated in Fig. 2.11, the highly collimated pulses generated by the laser are transmitted to a beam splitter which divides the signal in two parts. The fraction which is going to be transmitted is set on the optical path to the shooting device, while the rest is used to trigger the timer responsible for the measurement of the time of flight. The pointing control device orientates the telescope in order to shoot laser pulses towards the desired target and to receive the relevant return. Once the photo-detector acquires the reflected signal, the timer is stopped and the measurement is time-tagged with reference to U T C and stored.

For the sake of completeness, it should be mentioned that the primary ILRS station data product consists in normal points (NPs) which condensate the information of successive measurements. Full-rate original observations are, in fact, highly correlated over short time intervals, and forming NP allows reducing the observation noise while maintaining the substantial information inherent to the data. Each station is responsible to set minimum data criteria for computing NPs. Official ILRS guidelines require at least 6 data points (single shots) for daytime acquisitions, and 3 data points for nighttime observations in case of single photoelectrons systems with high data yield; stations with KHz return rates may select more stringent criteria.

-δρ bias represents the systematic station-and satellite-dependent difference between the theoretical and the measured range [START_REF] Thaller | Inter-technique combination based on homogeneous normal equation systems including station coordinates[END_REF]. These biases absorb both environmental effects which are not properly described by exploited models and actual systematic errors deriving from calibration issues, signal strength variations, lapses of the timer triggering mechanism or inaccurate characterization of the distance between the geometric reference point and the electric unit acquiring the signal.

-δρ CoM expresses the distance between the retro-reflectors and the Center of Mass (CoM) of the satellite. The measured time of flight, in fact, determines the ranging up to the optical devices on-board the spacecraft, but motion laws governing the orbit evolution require the satellite position X i to be referred to the Center of Mass. For SLR-specific geodetic targets, a careful pre-launch calibration of δρ CoM has been performed, but [START_REF] Otsubo | System-dependent center-of-mass correction for spherical geodetic satellites[END_REF] already showed that station-and satellite-specific corrections are needed to achieve ranging precision at the mm-level. CoM corrections are, in fact, shown to depend both on ranging systems and on the observation policy adopted by ground stations. For non spherical targets, determining δρ CoM is even more complicated because the effective reflection point becomes a function of the incidence angle [START_REF] Sośnica | Satellite laser ranging to GPS and GLONASS[END_REF]. In addition, most of these targets are active satellites with moving parts and fuel consumption throughout the mission lifetime. As a consequence, not only the location of the optical center, but also that of the CoM are more troublesome than for passive LAGEOS-like satellites.

-ǫ i r are random measurements errors While ρ rel are ρ tro are modeled in the same way for all kind of SLR observations, the target-dependent parametrization selected for δρ bias and δρ CoM are detailed in the following Sec. 2.3.2. It should be noted that, contrary to GNSS systems, no synchronization terms appear in Eq. 2.28 since just one timer is involved in the measurement of the pulse time of flight. The absence of clock and atmospheric parameters to estimate during the analysis, together with the stable orbits of the LAGEOS satellites explain why SLR is more sensitive to geocenter motion than GNSS. δρ bias and δρ CoM , on the other hand, directly correlate with station heights and accurate estimates are needed in order to provide high quality scale information.

Major shortcomings of the SLR technique are the sparse and unbalanced distribution of the ground network and the impossibility to establish the station-satellite optical link in bad weather conditions.

Analysis Strategy

SLR observations to the two LAGEOS satellites and to GPS/GLONASS vehicles carrying on-board retro-reflector arrays are analyzed in this study. For the first group of data, the ILRS has long established analysis standards that are carefully followed in the derivation of the presented results; for the second group, instead, no official guidelines have ever been issued. The strategy selected for this project, then, results from a mixture of references to the available literature, e.g. [START_REF] Thaller | Combination of GNSS and SLR observations using satellite co-locations[END_REF][START_REF] Thaller | Geocenter Coordinates from GNSS and Combined GNSS-SLR Solutions Using Satellite Co-locations[END_REF], and independent choices which will be detailed in the following.

SLR observations to LAGEOS

For SLR observations to LAGEOS satellites, the standard analysis recommended by the ILRS has been followed. Data are processed in weekly batches, provided that at least 10 N P per station are available. For each week, a single set of station coordinates is estimated imposing loose constraints with reference to the SLRF2008 positions which are used as a priori values. Contrary to the GNSS analysis, a single orbital arc is set-up over the whole interval. The parametrization selected for the orbit encompasses the 6 initial osculating elements, a constant acceleration in the along track direction, and two once-per-revolution accelerations both in along-track and cross-track directions. The relevant a priori information was kindly provided by the Italian Space Agency (ASI) as, until March 2016, no official orbital products were delivered by the ILRS. Concerning the EOPs, daily offsets are set-up for xpole, y-pole and LOD. The same a priori values exploited for GNSS analysis have been imposed. Range biases are set-up according to recommendations of the ILRS Analysis Working Group (AWG) which identified a time-dependent subset of stations for which range biases should be either estimated per pass or fixed to official values resulting from a multi-year adjustment of the observations. The AWG is responsible for maintaining the Data Handling File where all this information is detailed5 . Analogously, time-dependent δρ CoM corrections are applied according to ILRS estimates6 .

SLR observations to GPS and GLONASS

SLR observations to GPS-35, GPS-36 and all flying GLONASS satellites enable tying the GNSS and the SLR techniques via one of their possible co-location in space. As already explained in Sec. 1.4.2, the actual link is performed by the common esti-mation of orbital parameters with the contribution of both observation types. The same orbit parametrization and the same a priori values adopted for GNSS data must then be used in the derivation of the normal equation system relative to this last group of measurements. Station coordinates are consequently sampled at 9-day intervals. In agreement with the standard ILRS analysis, atmospheric corrections are modeled with the Mendes-Pavlis mapping function and daily offsets are set-up for x-pole, y-pole and LOD. Center of mass corrections are kept fixed to the ILRS official values; possible errors in the tabulated information are then mapped into the range bias parameters. Following the findings of [START_REF] Thaller | Combination of GNSS and SLR observations using satellite co-locations[END_REF], these lasts have been set up in the form of one offset per station and group of satellites. Given the limited number of SLR observations to high, complex targets such as GNSS satellites, estimating range biases at the 9-day level considerably increases the noise of the solution and it is not convenient to set-up one parameter per satellite-station pair, if empirical evidences suggest that this is not strictly necessary. Concerning the stations for which range biases have been estimated, two different solutions have been computed: in the first one, range biases have been set-up for all the stations, while in the second only the sites listed in the Data Handling File for the analysis of LAGEOS data are considered.

Dataset and result validation SLR observations to LAGEOS

The number of SLR observations to LAGEOS satellites is reported in Fig. 2.12 together with the volume of parameters estimated in each weekly solution. If GNSS observations constantly increased over the four years under study (Fig. 2.7), an average bulk of about 3000 LAGEOS acquisitions per week is maintained throughout the time span. Week-to-week oscillations in the amount of recorded data are however considerably more marked than those shown by GNSS, as a consequence of the SLR dependency on weather conditions.

Tracking statistics concerning annual yield, data quality, and operational compliance vary from station to station and only a limited subnetwork complies with ILRS criteria for the identification of high performing core sites7 , see Fig. 2.13 for the map of available fiducial stations.

The higher quality of the data provided by core systems is reflected by the WRMS of station position residuals computed during the long term stacking of weekly estimates. During the accumulation, seven transformation parameters have been set-up between each weekly solution and the derived long term frame; IC have been imposed on origin and scale parameters, while the orientation has been real-

E [mm] N[mm] U [mm]

All stat 6.24 6.71 6.12 Core stat 4.48 4.78 4.12 Table 2.5 -Average values of the time series displayed in Fig. 2.14.

As proved in Fig. 2.15, weekly estimates of core station coordinates also show a better agreement with the ILRS official solution submitted for the computation of ITRF2014. The displayed residuals result from a 7-parameter comparison. Estimated values are systematically 1.5/2-times larger when the whole network is used.

The frame information expressed by SLR observations to LAGEOS satellites has also been compared to ITRF2008. For this purpose, a long term stacking has been computed imposing MC on all the seven transformation parameters set up between each input solution and the resulting secular frame. The MC conditions have been implemented exploiting core stations only. The time series of translation and scale factor estimates are presented in Fig. 2.16, average values are reported in Tab. 2.6. For comparison purposes, the same parameters have been additionally evaluated on the basis of SLR contributions to the computation of ITRF2014. In particular, the combined ILRS solution and the independent estimates provided by a subset of the official ACs 8 have been considered. As shown in Tab. 2.6, results obtained in this study are within the variability observed for the estimates of the different official ACs. Concerning the translation, the worst agreement is found in T z which, analogously to the GNSS case, is the most poorly determined component. The observed lower quality is consistent with the slight collinearity issue between T z and orbital osculating elements detected by Rebischung [2014, App. 4.4.3]. A few mm offset is also observed in T y , in accordance with the discontinuity which can be identified in 2010 when screening the entire ILRS time series [START_REF] Altamimi | Preliminary results and ilrs contribution[END_REF] A.C.

T x T y T z scale [mm] [mm] [mm] [ppb]
This study -0.84 -2.49 4.77 -0.2 ILRS sol.

-0.15 -2.99 3.00 -0.5 BKG sol.

-0.57 -2.54 2.92 -0.53 ASI sol.

-1.80 -3.06 3.41 -0.23 ESA sol.

-1.74 -2.86 2.70 -0.9 GFZ sol. 0.80 -3.96 1.02 -0.18

Table 2.6 -Mean SLR origin and scale offsets with reference to ITRF2008.

SLR observations to GNSS satellites

The number of SLR observations to GNSS satellites is presented in Fig. 2.17. A clear increment is observed from August 2014, when the first LAser Ranging to GNSS s/c Experiment (LARGExperiment) took place [START_REF] Noll | Summary of Results from ILRS GNSS Tracking Campaigns[END_REF]. During this campaign, ILRS stations have been asked to modify their usual schedule in order to improve GNSS-tracking performances. Different strategies have then been explored in successive experiments in order to identify the best possible operational policy. Also for observations to GNSS satellites, ranging statistics vary considerably from station to station, as can be seen in Fig. 2.18 presenting the total number of data acquired at each site between January 2011 and December 2014. The three best performing stations, namely Zimmerwald (7810), Yarragadee (7090) and Graz ( 7839), collected about 45% of the data; the four best stations provided more records than the rest of the network combined. Analogous findings have been obtained by [START_REF] Sośnica | Satellite laser ranging to GPS and GLONASS[END_REF], for the period 1994 -2013. Additional plots showing the time series of daily acquisition at each site are presented in App. A. Validating the results obtained from SLR observations to GNSS satellites is not straightforward. On the one hand, ILRS does not provide any official product based on the analysis of such data; on the other, for poorly performing SLR stations the number of acquisitions is insufficient to set up a non singular normal equation system including station coordinates, range biases, EOPs and the complete orbit parametrization selected for GNSS satellites. A trial solution encompassing just station positions and range biases has then been computed and compared to the official ILRS contribution to the computation of ITRF2014. Results are presented in Fig. 2.19. A general agreement at the 2 cm level is observed.

Summary

The analysis strategies developed for this study provide frame results of comparable quality with reference to the products delivered by the official ACs of the dedicated IAG Services. IGS and ILRS solutions computed in the framework of the reprocessing campaign for the realization of the ITRF2014 have been used as reference term for the validation.

Concerning the analysis of GNSS observations, the inclusion of a high percentage of GLONASS data in a rather sparse network still limits the accuracy of the computed solution. In this study, it has been necessary to implement an orbit parametrization relying on three-day arcs in order to obtain station coordinate estimates comparing to IGS official products within a few mm in the horizontal components and 4 -5 mm in the vertical. Improvements on shorter time-spans might result, instead, from the exploitation of the newECOModel [START_REF] Arnold | CODEś new solar radiation pressure model for GNSS orbit determination[END_REF] or of a box-wing approach in the description of the forces acting on the spacecrafts.

A remarkable agreement is found between the results presented in this study and the repro2 CODE's solution (which also encompasses GLONASS data and was obtained with the Bernese Software). Average translation and scale offsets computed between these two solutions and ITRF2008 compare at the mm level, with the only exception of T z . The discrepancy observed for this component can still be considered reasonable, given the noise level of the time series and marked spurious oscillations. The behavior shown by T z is likely to reflect just modeling deficiencies and observation errors rather than actual geophysical signals. Besides being theoretically explained by collinearity issues [START_REF] Rebischung | Can GNSS Contribute to Improving the ITRF Definition?[END_REF], such insensitivity is also reflected in the direct comparison with the equivalent time series derived from LAGEOS observations, see Fig. 2.20. The other components, on the other hand, are characterized by a seasonal signal in phase with the one sensed by SLR and even the amplitude is nicely reproduced in case of T x . In the interpretation of the results, however, it should always be considered that the network shift approach allows studying the non linear part of the geocenter motion only up to the possible time variations of the CF-CN vector. Spurious signals depending on the selected observation geometry, network effect, are likely to be aliased in the estimated parameters.

Concerning the scale factor, only a small average offset (∼ 0.2 ppb) is observed with reference to ITRF2008. The adoption of IGS official calibrations for satellite P COs conventionally force the GNSS-realized scale into the ITRF2008 one, but [START_REF] Rebischung | The IGS contribution to ITRF[END_REF] pointed out that modeling updates applied in the repro2 campaign might explain the observed variation. with reference to ITRF2008. To facilitate the comparison, both time series have been detrended and smoothed with a Vondrak filter with a 7 cpy cutoff frequency [START_REF] Vondrak | A contribution to the problem of smoothing observational data[END_REF].

The analysis of SLR observations to LAGEOS satellites has shown that, even in the most recent years, tracking performances are strongly station-dependent. Coordinates computed for core stations show a higher week-to-week consistency and a better agreement with the ILRS official solution, as demonstrated by the reduction of the relevant WRMSs displayed in Fig. 2.14 and 2.15.

The long term frame computed from LAGEOS data collected during the period 2011 -2014 shows non negligible translation offsets with reference to ITRF2008 in the y-and z-direction. If it should be considered that T z is the noisiest and worst determined component, the discrepancy observed for T y is consistent with the discontinuity affecting the relevant historical time series in 2010. An analogous jump, also characterizes the historical series of the scale factor. In this case, however, successive observations tend to restore pre-offset values and, according to the ILRS combined solution, data collected in the period 2011 -2014 show the expected half a ppb offset with reference to ITRF08. Further investigations are, anyway, required to understand the nature of the detected discontinuities in order to reliably assess the realization of the frame defining quantities. As a side note, it can be pointed out that the scale parameter is the one associated with the largest variability in the estimates provided by the different ACs.

Finally, SLR observations to GNSS satellites have been taken into account. Also in this case, tracking performances are dramatically station-dependent and it has been shown that the four best performing stations acquire more than one half of the total volume of data. For poorly performing stations, it is even impossible to compute a fully parametrized solution over the three-day time span imposed by the selected orbit estimation process. A partial validation of the developed analysis strategy was then performed computing a trial solution encompassing just station positions and range biases. On average, the agreement with the SLR official solution is found at the level of a couple of cm.

Chapter 3

Results

GNSS and SLR solutions computed according to the analysis strategies detailed in Chapter 2 can be employed for the derivation of combined TRFs relying exclusively on the space ties on-board GNSS satellites. In this chapter, the results of both long-term and quasi-instantaneous computations are presented and discussed. Future perspectives are also investigated by means of simulations addressing the enhancement of frame parameter precision resulting from technological improvements and/or strengthening of the ground network.

The following questions are confronted in this chapter:

1. Can the space tie on-board GNSS satellites provide a reliable and effective alternative to the use of terrestrial ties in the derivation of Terrestrial Reference Frames?

2. Can this link transfer the origin and scale information from one technique-specific frame to the other?

3. How does the frame realization based on the space tie approach compare with the traditional one?

4. What gain in the precision of frame parameters can be expected in consequence of technological improvements?

5. Should GNSS providers be encouraged to add laser retro-reflectors on all GNSS satellites ?

General combination scheme

Results presented in this thesis are obtained following the multi-technique combination scheme detailed in Fig. 3.1. The synthesis of the different pieces of information is carried out at the normal equation level, but it is reminded that all the available observations are treated with the same analysis software, ensuring the highest possible consistency in data handling. The proposed strategy, therefore, approximates the combination at the observation level, as already discussed in Sec. 2.1.

Daily microwave and optical observations concerning GNSS satellites are separately pre-processed and the corresponding normal equation matrices are derived according to the parametrization detailed in Chapt. 2. As orbit parameters and EOPs are set up in both systems, these can be stacked according to Eq. 2.14. The weights employed in the combination are computed following Eq. 3.1, presented in the next section. Three consecutive pre-combined systems are then accumulated in a unified solution relying on a single orbital arc. Orbit parameters are then reduced and 9-day normal equations centered on the ILRS week are derived from the stacking of three consecutive non-overlapping 3-day systems. The sampling of the obtained time series is, therefore, comparable to the one selected for the analysis of LAGEOS data and the overall combination can finally be computed.

Two different approaches have been followed. In one case, the combination has been performed at the weekly level, stacking the normal equation matrices referred to the same mean epoch according to Eq. 2.14. In the other, an ITRF-like combination has been realized: distinct long term solutions have been independently computed from the two time series and the resulting frames have been successively stacked following Eq. 1.15. Given that both input solutions encompassed SLR station positions and EOPs, it was not necessary to exploit terrestrial ties. Once again, the different contributions have been weighted according to Eq. 3.1. During the stacking, transformation offsets have been set up between each GNSS+SLRtoGNSS normal equation and the combined solution, in order to transfer to the final frame the origin and scale information realized by LAGEOS observations. The presented approaches reflect two different philosophies that can be followed for the datum definition. When realizing quasi-instantaneous frames, in fact, all the available observations, suitably weighted, are taken into account for the realization of the defining parameters. In the ITRF-like process, instead, the realization of origin and scale relies exclusively on LAGEOS observations.

Weights

During the inter-technique combination, the different contributions have been weighted according to:

w i = 1 σ 0i σ 0GN SS 2 • χ 2 i DOF (3.1)
where:

-σ 0i , a priori σ of unit weight, is a user-selected parameter which should represent the standard deviation of the i-th group of input data (since it is used for this purpose in the definition of the weight matrix employed in the least square adjustment) -

χ 2 i
DOF , a posteriori unit variance factor, expresses the consistency between adjusted parameters and original observations

For each contributor, therefore, the weight to be used in the combination is computed according to the noise of the raw observations scaled for a measure of the fit quality. If σ 0i captures the variability of the match between observations and estimates,

χ 2 i DOF ∼ 1; bigger(/smaller
) values indicate that σ 0i has been under(/over) estimated. Accounting for the a posteriori unit variance factor in the definition of the weights, then, ensures that the diverse inputs are scaled according to a realistic estimate of the relevant observation noise.

The values of σ 0i assumed in this study for the different observation types are reported in Tab. 3.1; the computed time series of χ 2 i DOF are shown in Fig. 3.2. For the adjustment of SLR observations to GNSS satellites the trial solution described in Sec. 2.3.3 was exploited. Derived weights are, therefore, somewhat optimistic considering that the number of DOF is increased. However, the singularities arising in fully parametrized normal matrices make it impossible to compute a complete solution without including GNSS observations. Observation type σ 0i [m] GNSS 0.001 SLR to LAGEOS 0.01 SLR to GNSS 0.05 Table 3.1 -A priori σ of unit weight. The reported time series of χ 2 i DOF show that SLR results are more scattered than GNSS ones, reflecting the greater week-to-week variability of the network geometry and of the amount of collected observations. GNSS values, on the other hand, present a clear seasonal cycle which can be related to the unbalanced number of stations located in the two hemispheres. This can be clarified analyzing baselinewise the time series of observation residuals: for stations located in the Northern hemisphere, a seasonal variation is detected which is in phase with the cycle shown in Fig. 3.2 -panel a, while the opposite holds for Southern sites. Given that, in the selected network, stations in the Northern hemisphere outnumber those in the Southern, the effects cannot cancel out.

As an example, relevant results evaluated for the IGS stations GODZ (Greenbelt, Maryland, U.S.) and ALIC (Alice Spring, Northern Territory, Australia) are reported in Fig. 3.3. In the daily double-difference adjustment, GODZ is typically coupled with the co-located site GODE and with USNO and BRMU stations that are installed, respectively, at the US Naval Observatory in Washington and in the Bermuda Island. ALIC is generally associated to the "nearby" stations DARW, STR1 and YAR3 in Darwin (Northern Territory), Canberra (Capital Territory) and Yarragadee (Western Australia). The signals portrayed in the two panels clearly show inverted seasonality. The flat line in GODZ plot is referred to the very-short GODZ-GODE baseline; the stability of the pertinent residuals derive from the absence of differential loadings and atmospheric artifacts between these co-located sites.

RMS of observation residuals [mm]

RMS of observation residuals [mm] In any case, the observed temporal variations in

χ 2 i DOF
, either in the form of weekto-week scatter or of seasonal cycle, are negligible when compared to the differences observed for technique-specific a priori unit weights. In the computation of the weights, then, only the average values of

Results

The results of the combination are presented and discussed in the following sections. Products of the stacking at the weekly level are considered in Sec. 3.3, while Sec. 3.4 focuses on the ITRF-like combination. As anticipated in Sec. 2.3.2, two different solutions have been computed in order to investigate the effects of possible inaccuracies in the official values of space tie vectors:

-CASE A: SLR range biases to GNSS satellites estimated only for those stations indicated by the ILRS Analysis Working Group (AWG) in the framework of the adjustment of LAGEOS data. In this case, the computation basically relies on nominal space tie vectors, computed as the sum of Laser Ranging Arrays (LRA) offsets and GNSS satellite PCOs recommended by ILRS and IGS respectively.

-CASE B: SLR range biases estimated for each station and each group of satellites (i.e. two biases per station: one for the two GPS and one for all the GLONASS spacecrafts). It should be pointed out that, with such a parametrization, biases absorb both station-dependent mismodeled tracking effects and inconsistencies in the tabulated tie vectors. The choice of setting up a single parameter per station and group of satellites is justified by the theoretical assumption that range biases are primarily affected by ranging distance, tracking policies and technical characteristics of the retro-reflectors [START_REF] Otsubo | Centre-of-mass correction issues: Toward mmranging accuracy[END_REF]. Findings of [START_REF] Thaller | Combination of GNSS and SLR observations using satellite co-locations[END_REF] further support this approach. [START_REF] Sośnica | Satellite laser ranging to GPS and GLONASS[END_REF], on the other hand, indicate that satellitedependent effects might also affect the ranging, but setting up these additional parameters might be critical for stations acquiring just a limited number of returns from GNSS satellites.

Quasi-instantaneous approach

In order to evaluate the linking efficiency of the space ties on-board GNSS satellites, the results of the combination at the weekly level have been compared to technique-specific adjustments carried out at the same epoch. A 9-day GNSS solution based on 3-day orbital arcs has been explicitly computed for this purpose, while LAGEOS weekly results, Fig. 3.1 -B, have been exploited as a reference for the SLR technique. The comparison has been performed estimating the full set of transformation parameters and the full variance/covariance matrix has been employed for the weighting.

Translation and scale offsets estimated between the frames realized, on a weekly basis, by GNSS-only solutions and by the GNSS subnetwork of the fully combined adjustment are presented in Fig. 3.4; equivalent results are obtained for both CASE A and CASE B solutions. Analogous transformation parameters computed for the SLR network are reported in Fig. 3.6 and 3.7 for CASE A and CASE B respectively. The same scale is employed for the y-axis of all plots. It can be noted that the frame realized by the GNSS subnetwork in the fully combined solution does not differ at all from the one computed on the basis of GNSS-only data. This means that space ties are not conveying any LAGEOS-derived information to the GNSS subnetwork.

The results concerning translation offsets are in line with findings of [START_REF] Thaller | Geocenter Coordinates from GNSS and Combined GNSS-SLR Solutions Using Satellite Co-locations[END_REF]. In that study, geocenter coordinates were set-up as explicit parameters in the observation adjustment and different time series were derived using either microwave-only or microwave plus optical data. For orbital parametrization choices comparable to those implemented in this work, no differences were found on the geocenter estimates carried out introducing or neglecting LAGEOS observations and exploiting the space ties on board GNSS satellites as the only inter-technique link. For the sake of completeness, it is reported that the alternative orbit parametrization investigated in that study envisaged the estimation of once-per-revolution empirical orbit parameters without applying any constraints. In that case, the addition of LAGEOS data contributed to the reduction of draconitic signals in the geocenter time series, but just up to the level which can be achieved setting up the recom-mended constraints.

After the findings of Rebischung et al. [2014], it can be further commented that SLR observations to GNSS satellites do not directly impact the mechanisms which determine the GNSS insensitivity to geocenter motion. The latter is, in fact, mainly due to the necessity of simultaneously estimating satellite clock parameters and (microwave) atmospheric delays, to which SLR data provide no additional information. A limited improvement could have been expected for the z-component, since this has been shown to partially correlate with a subset of ECOModel parameters [START_REF] Rebischung | Can GNSS Contribute to Improving the ITRF Definition?[END_REF]. In operational condition, however, GNSS data are so much more abundant and precise than SLR observations to GNSS satellites that they dominate the orbit estimation process and no detectable effects are found on T z .

Concerning the scale, then, it should be reminded that by fixing P COs to official IGS08 values, the GNSS scale is conventionally inherited from ITRF2008. If zP COs are freed in the fully combined solution, a mean scale bias of ∼ 0.6 ppb is observed with reference to the GNSS-only standard solution. The result is in excellent agreement with the multi-year adjustment presented in [START_REF] Thaller | Precombined GNSS-SLR Solutions: What Could be the Benefit for the ITRF[END_REF]. When estimating zP COs in addition to the standard parametrization, the scale realization of the GNSS network becomes almost singular. If the co-location on board GNSS satellites offered an efficient link between the two techniques, the scale information as materialized by SLR observations to LAGEOS would be transferred to the combined solution and then mapped into GNSS zP COs. To test this hypothesis, a new reference GNSS solution has been computed freeing zP COs. From the numerical point of view, the derived normal equation can be inverted because it is almost, but not exactly singular, see Rebischung [2014, Fig. 4.3]. The derived results are, however, rather meaningless as they just basically reproduce modeling errors and observation uncertainties. When the fully combined solution with free zP COs is compared to this new reference solution, the previously observed bias vanishes completely, Fig. 3.5. In the fully combined case, then, no LAGEOS scale information was transferred to the GNSS frame via the space ties.

No relevant frame discrepancies are observed also when focusing on the SLR subnetwork, Fig. 3.5 and 3.6. Larger result uncertainties are observed, mainly because of the contribution of noisy SLR observations to GNSS satellites, but no significant mean offsets are observed. It is interesting to note that discrepancies in T z are less stable when range biases to GNSS satellites are not estimated. As the z-component of the geocenter motion is the worse determined by SLR, modeling errors and inconsistencies show up more easily in this parameter. In CASE B solution, range biases seem to systematically absorb the sources of frame discrepancy. This further confirms the original findings of [START_REF] Otsubo | GLONASS laser ranging accuracy with satellite signature effect[END_REF], suggesting that all SLR stations provide biased observations when tracking GNSS satellites. In that study, however, only SLR observation residuals to fixed orbits were considered. The presented plots show that, in the combined solution, each technique-specific subnetwork keeps its own independent frame realization. The selected space ties are not providing any actual link, and the two realizations coexist autonomously within the "combined" adjustment. To better understand the implications of these findings a 1D sketch representation of the obtained results is provided in Fig. 3.8. In the drawing it is assumed that an SLR and a GNSS station are installed at the same distance x from the geocenter, Fig. 3.8 -a. In each technique-specific frame, different (biased) realizations of the center of mass and of the scale are materialized, implying that the actual coordinates of the two points in each relevant system are, in general, different from one another even if the stations are physically co-located in the same point Fig. 3.8 -b. The application of effective ties during the combination leads to a unified frame in which all subnetworks rely on the same physical parameters. Due to systematic errors and modeling deficiencies, the materialized origin and scale will still be biased, but the relative positioning of points will be correctly reproduced, at the expenses of specific-network rigid transformations and/or deformations Fig. 3.8 -d. When the employed ties fail to actually link the two frames, however, frame information are not homogenized and, even within the final "combined" solution, stations keep positioning according to their own intrinsic frame realization Fig. 3 (c) When employed ties fail to actually link the two frames, stations keep positioning, in the "combined" solution, according to their own intrinsic frame realization. (d) In properly combined frame, instead, all subnetwork adopt the same origin and scale and spatial co-location is restored. The author is grateful to P. Rebischung for providing the original idea for this representation.

.8 - c. CM X GNSS X SLR X SLR X GNSS X SLR X GNSS a) b) c) d)

Long-term approach

The quality of the GNSS-SLR space ties on-board GNSS satellites have been further tested on the basis of the computation of secular reference frames. As described in Sec. 3.1, these frames have been computed stacking the long term accumulations of the LAGEOS-only solutions and of the pre-combined time series derived from GNSS data and SLR observations to GNSS satellites. The origin and the scale of the LAGEOS solution have been transferred to the combined frames. These lasts have been evaluated against the results of a traditional ITRF computation, Sec. 3.4.2 and have also been checked for internal consistency, Sec. 3.4.3.

A piece-wise linear parametrization has been assumed for modeling station coordinates; the list of discontinuities selected for the computation of ITRF20141 has been exploited for this purpose. Such list results from a visual inspection of station position residuals, in agreement with the findings of [START_REF] Gazeaux | Detecting offsets in GPS time series: First results from the detection of offsets in GPS experiment[END_REF] showing that manual methods typically perform better than automated ones in the detection of jumps in GNSS time series. With the growing number of GNSS receivers around the world, however, the development of reliable automated screening procedures shall also be encouraged. In this context, the potential of the blind STARS algorithm [START_REF] Rodionov | A sequential algorithm for testing climate regime shifts[END_REF] are investigated in App. C.

The long term stacking of the pre-combined GNSS+SLRtoGNSS time series has been performed estimating seven transformation parameters between each nine day solution and the resulting secular frame. Internal constraints have been applied to the realization of the origin and the scale, while the orientation singularities have been recovered through minimum constraints, with reference to ITRF2008, implemented for a fiducial subset of well distributed GNSS stations. Two different solutions have been computed either imposing or neglecting Equal Velocity Constraints (EVC) for co-located ground sites. Strictly speaking, only the latter implementation follows the pure space tie approach; the former, on the other hand, provides a closer analogy with the classical ITRF computation and is based on the reasonable assumption that no relevant intra-site motions should be detected. Exceptions to such a hypothesis may be related to monument instabilities and particularly unstable geological environments [Sarti et al., 2013b].

Comments on the estimation of range biases to GNSS satellites

All the results that will be presented in Sec. 3.4.2 and 3.4.3 have been computed considering CASE A parametrization for the analysis of SLR observations to GNSS satellites. In order to perform consistent computations for CASE B solutions, in fact, it would first be necessary to fix range biases to reliably estimated values. An attempt was made to retrieve such quantities from the long term stacking of the fully combined weekly frames derived from CASE B solutions, Sec. 3.3. The obtained results are presented in Fig. 3.9. As a consequence of the large evaluated uncertainties, the estimated GLONASS RBs are significantly different from zero just for about half of the SLR stations; concerning the GPS, a clear indication of the RB value is found in a few cases only. SLR observations to GNSS satellites are likely to be biased for virtually all the stations [START_REF] Otsubo | GLONASS laser ranging accuracy with satellite signature effect[END_REF][START_REF] Thaller | Combination of GNSS and SLR observations using satellite co-locations[END_REF] , but the values computed within this study have been judged of too poor quality to justify a back-substitution in the original solution.

The biggest error bars observed in Fig. 3.9 correspond to estimates based on sparse datasets, confirming that setting up one RB per satellite might be critical for poorly performing stations. The accuracy of the results might, instead, benefit from the availability of satellite-specific corrections for LRA offsets and from modeling the elevation-dependency of the satellite signature effect [START_REF] Sośnica | Satellite laser ranging to GPS and GLONASS[END_REF].

Comparison with traditional ITRF computation

In order to investigate how the space tie approach compares to the traditional one, a reference GNSS-SLR frame was computed exploiting the set of terrestrial ties selected for the realization of ITRF2014. The combination was realized on the basis of a GNSS-only and a LAGEOS-only secular frames; SLR observations to GNSS satellites were not considered. Technique-specific stackings were performed estimating seven transformation parameters between each input adjustment and the resulting cumulative solution. Internal constraints were applied to the realization of the physical frame parameters while minimum constraints with reference to ITRF2008 were imposed on rotations. Following [START_REF] Altamimi | an improved solution of the international terrestrial reference frame[END_REF][START_REF] Altamimi | A new release of the international terrestrial reference frame for earth science applications[END_REF], an uncertainty of at least 3 mm was assigned to each tie vector component and the different ties were weighted so that no normalized residuals larger than 3 were detected. The two technique-specific networks were additionally linked through equal velocity constraints imposed at co-located ground sites. During the inter-technique combination, fourteen transformation parameters were set up for the GNSS solution, so that SLR observations to LAGEOS satellites were responsible for the reference frame definition. Analogous transformation offsets were set up also for the pre-combined GNSS + SLRtoGNSS long term solution during the combination with the secular LAGEOSonly frame. The resulting offsets and the relative rates are reported in Tab. 3.2 and Tab. 3.3 respectively. In both tables, the upper entry is referred to the classical computation, while middle and lower ones report the results of the combination based on the space tie approach either imposing or neglecting equal velocity constraints for co-located ground stations. The consistency between the fully combined frame and the input ones can be evaluated from Tab. 3.4 where the WRMS of station position and velocity residuals are listed for each group of observations.

T 1 T 2 T 3 D R 1 R 2 R 3 [mm] [mm] [mm] [ppb] [mas] [mas] [mas]
Terrestrial ties -5.7

1.3 0.1 -0.24 -0.05 -0.04 0.07 (1.8) (1.7) (1.7) (0.26) (0.07) (0.07) (0.07) Space ties -3.9 1.3 -9.2 0.52 -0.05 -0.23 -0.01 EVC (2.2) (2.2) (1.9) (0.26) (0.13) (0.10) (0.13) Space ties -4.2 1.6 -8.7 0.54 -0.08 -0.15 -0.04 NO EVC

(2.0) (1.9) (1.6) (0.22) (0.10) (0.08) (0.10) In order to determine whether the space ties on-board GNSS satellites offer a link of comparable strength with reference to the traditional approach, it is necessary to investigate the formal errors of the transformation offsets and of their relative rates. These values, in fact, provide an estimate of the uncertainties with which the frame information realized by SLR observations to LAGEOS satellites is transferred to the other input solution. Table 3.3 shows that, when equality constraints are imposed on the velocities of co-located stations, the rate errors are comparable for both tying approaches. If EVC are ignored, however, the uncertainties increase of a factor of 3, indicating that the combination which rigorously relies on the space co-location is considerably weaker than the one taking advantage of terrestrial links. Together with their uncertainties, rate values also increase.

The weakness of the investigated space link is also reflected in the WRMSs of station position and velocity residuals computed for the LAGEOS-only frame, Tab. 3.4. As shown in Fig. 3.8, the reduction of independent solutions into a combined adjustment requires, in general, a certain (limited) deformation of the original solution. The SLR network realized on the basis of LAGEOS observations is systematically associated to lower residuals when space ties are exploited in place of terrestrial ones. Relevant values even drop to zero when EVC are not applied. The small inconsistencies observed for the GNSS network in case of a traditional frame computation also disappear when space ties are exploited. The high WRMS values observed for the frame realized by the SLR subnetwork in the GNSS+SLRtoGNSS frame reflect, instead, the high noise of such solution and its poor comparison with the standard adjustment of LAGEOS results, see Fig. 2.19.

The uncertainties associated to transformation offsets are found to be comparable for both tying approaches. The interpretation of this result is however delicate. In the original time series, in fact, the origin and scale information of each input solution was only loosely determined. Before running the long term stacking, then, suitable constraints have been applied taking the precision of the datum definition at the 0, 1 mm level. GNSS stations have been used as a reference in this phase. Such pre-processing measure is functional to the reliable application of internal constraints during the accumulation, but perturbs the original statistical information concerning the precision of the datum realization. As these constraints have been applied to both GNSS-only and GNSS+SLRtoGNSS time series, the observed comparable uncertainties might, to some extent, be artificial. Further investigations are demanded in order to unambiguously characterize the effects of the preprocessing.

Concerning the values of the estimated transformation offsets, the higher discrepancies with the traditional approach are found for T z and for the scale parameter. In principle, differences are acceptable given the different nature of the ties linking the exploited techniques. It should be reminded, however, that, since the CASE A parameterization is exploited, possible errors and inconsistencies in the space tie vectors cannot be absorbed by suitable parameters and necessarily map into transformation offsets. In addition, it should be considered that the SLR solution based just on observations to GNSS satellites compares to the LAGEOS solution only at a few cm level, see Fig. 2.19.

Internal consistency evaluation

Besides the comparison with the traditional ITRF computation, the internal consistency of the stacking can be checked. If the selected space ties performed a successful link, SLR and GNSS stations should be expressed in the same reference frame within the GNSS+SLRtoGNSS solution. Seven transformation parameters have then been estimated at epoch 2013.0 between the overall combined frame and the technique-specific subnetworks within the GNSS+SLRtoGNSS long term adjustment. Transformation parameters comparable to those listed in Tab. 3.2 were expected for both groups. As it can be seen in Tab. 3.5, however, only the GNSS subnetwork comply with such expectations while results computed for SLR stations present evident inconsistencies.

T 1 T 2 T 3 D R 1 R 2 R 3 [mm] [mm] [mm] [ppb]
[mas]

[mas]

[mas] EVC GNSS subnet -3.9 1.3 -9.2 0.52 -0.049 -0.233 -0.010 (0.0) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0) SLR subnet -6.6 -0.7 10.1 -0.93 -0.001 0.055 -0.005 (4.2) (4.1) (4.0) (0.63) (0.159) (0.175) (0.161) NO EVC GNSS subnet -4.2 1.6 -8.8 0.55 -0.082 -0.15 -0.036 (0.0) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0) SLR subnet 3.0 -11.2 14.9 -1.30 -0.340 -0.149 0.280 (5.3 ) (5.2) (5.1) (0.80) (0.202) (0.222) (0.203) Table 3.5 -Translation T, scale D and orientation R offsets estimated between the fully combined secular frame and the technique-specific subnetworks within the GNSS+SLRtoGNSS long term stacking. Relevant formal errors are indicated within brackets.

Simulations

The co-location of GNSS and SLR techniques on-board GNSS satellites is, in principle, particularly appealing for the computation of the ITRF as it allows to directly link two of the four contributing systems relying exclusively on satellites already encompassed in the current analysis. The approach might even support the investigation of technique-specific biases improving the accuracy of the derived results.

In previous sections it has been shown that, at present, the quality of the link provided by such space co-location is insufficient to replace terrestrial ties. It is however worth investigating whether improvements in tracking performances might strengthen the provided connection and eventually realize a reliable inter-technique bond. In order to investigate this possibility, a series of simulations has been performed with the software originally developed in the framework of [START_REF] Rebischung | Can GNSS Contribute to Improving the ITRF Definition?[END_REF].

According to [START_REF] Kuang | Analysis of orbital configurations for geocenter determination with GPS and low-Earth orbiters[END_REF], "the observability of a quantity in an observation system is reflected by its corresponding variance, or formal error, in the full covariance matrix obtained through the associated estimation process". The synthetic study presented in this thesis, then, has been carried out envisaging a series of different tracking scenarios, see Sec. 3.5.1, and deriving the associated normal equa-tion matrices with reference to the parametrization exploited in standard analysis. Once the design and the minimal constraint matrices are, respectively, computed and set up, the variance/covariance matrix of the system can be retrieved according to Eq. 1.22. Concerning implicit parameters 2 , such as frame origin and scale, the relevant uncertainties can be evaluated according to [START_REF] Sillard | A review of algebraic constraints in terrestrial reference frame datum definition[END_REF], see also Rebischung [2014, Sec. 3.2.4].

Before presenting the implemented simulations and the relative results it is necessary to add a relevant consideration. Rephrasing the statement of [START_REF] Kuang | Analysis of orbital configurations for geocenter determination with GPS and low-Earth orbiters[END_REF], variations in the information content brought to the computation of each parameter in distinct observing scenarios imply changes of the associated variances. The reasoning, however, is based on the assumption of normally distributed errors and, therefore, it applies exclusively to the precision attained for each quantity. The accuracy of the estimates, which is directly impacted by systematic biases, cannot be addressed with this approach.

Investigated scenarios

Thirty different simulations have been carried out. Three of them encompass exclusively GNSS observations, see Tab. 3.9, and serve as a reference for assessing the precision improvements deriving from the additional consideration of SLR data to LAGEOS and GNSS satellites. A further reference solution, presented in App. B, was also computed on the basis of laser returns from LAGEOS spacecrafts. The general options implemented in the simulations are presented in Tab. 3.6, while the characteristics adopted for the 26 multi-technique cases are reported in Tab. 3.7.

For both GNSS and SLR ground networks, a homogeneous distribution was implemented, notwithstanding the actual partitioning between oceans and inland areas. This choice ensures that the estimated uncertainties are independent from any network effect. The adopted daily time-span appears in conflict with the recommended weekly-basis for the derivation of SLR products. ILRS guidelines, however, are intended to ensure that enough data are available for the analysis in order to provide reliable, robust results. In these simulations, however, a wealth of information is guaranteed even on a daily basis. The same observation rate, in fact, is adopted for all data types meaning that, as long as the satellite is in the station field of view, a NP is formed every 5 min. This assumption is far from being realistic for SLR, especially concerning day-time acquisitions and/or high targets such as GNSS satellites, but it guarantees a wealth of information even on a daily basis. In 2 Following Rebischung [2014, Sec. 3.2.1] implicit parameters are those quantities y with respect to which the partial derivatives of the explicit model parameters x are a known vector λ = ∂x ∂y of R p addition, this choice simplifies the interpretation of the results in term of changes of the variable under study, namely the SLR network dimension and the achieved precision of ranging data. 3.6 -General options implemented in the simulations The parametrization selected for the simulations is presented in Tab. 3.8. For the sake of simplicity, no radial biases are, in general, estimated for SLR observations; the only exception is case B25pl, for which LRA offsets for GNSS satellites are set up. Simulated systems are additionally provided with the constraints necessary for coping with known singularities. As introduced in Sec. 2.2.2, the list encompasses no-net-rotation conditions and constraints on the DUT1 parameter. In addition, since clock offsets are explicitly estimated, also epoch-wise zero-mean requirements must be implemented, see Rebischung [2014, Sec. 2.1.2 The selected case identifiers consist of a capital letter specifying the block followed by the number of stations included in the SLR network and, if pertinent, by a 1/2-digit code referring to some peculiar characteristic of the solution. In particular m stands for multi-GNSS cases (GPS+GLONASS), while p refers to the estimation of zPCOs.

Aspect

Results

In this section, the formal errors of the physical frame parameters computed according to the different simulated scenarios are presented. Table 3.9 lists the reference values computed for three simulations encompassing exclusively GNSS observations. In the first case, only GPS data are considered; in the second, GLONASS observations are added; in the third, the estimation of zPCOs is introduced. For each of these options, the subsequent plots, Fig. 3.10, 3.11 and 3.12, show the relative precision improvement which can be obtained adding SLR measurements of increasing precision and abundance. Presented values refer to the frame parameters as realized on the basis of the GPS/GNSS subnetwork only. 3.9 -Formal errors of the physical frame parameters computed for simulations based exclusively on GNSS data. These results are taken as a reference for the evaluation of the effects resulting from the inclusion of additional SLR data to the computation.

GPS

According to the results reported in Tab. 3.9, adding GLONASS observations considerably reduce the formal errors on all the estimated parameters of about 30 -35%. For the most part, this trivially follows from the increment in the number of observations, but [START_REF] Rebischung | Can GNSS Contribute to Improving the ITRF Definition?[END_REF] also pointed out that exploiting different constellations contributes in reducing collinearity issues. In particular, stationdependent parameters are common to both systems and are, therefore, less able to absorb the second order signatures of the geocenter (see his Sec. 4.5.4).

The estimation of zP COs, has negligible effects on the precision of the origin realization, while it degrades the formal errors of the scale factor of about one order of magnitude. As no-net-scale constraints have not been set up in the simulation, a corruption was expected given the almost perfect correlation of the scale parameter with the simultaneous estimate of tropospheric zenith delays, clock offsets and zP COs.

handled for all the contributing systems.

When zP COs are estimated, it can be observed that the quasi-singularity of the scale of the GNSS network is not recovered by the addition of SLR observations. For cases presented in Fig. 3.12, in fact, the precision of the SLR realization of the scale parameter3 systematically exceeds that of the GNSS realization. If the space tie was able to transfer the scale information from one network to the other, more relevant precision gains would be expected.

The modest observed improvements are consistent with the evidence that the GNSS weak ability to sense the geocenter (and the scale) depends primarily from the collinearity issues related to the simultaneous estimation of tropospheric and clock parameters (and zPCOs). None of these quantities, however, is included in the SLR observation equation, which means that SLR measurements do not provide any direct information to their determination. On the other hand, the partial correlation of T z with orbit parameters explains why this component is the one for which the greatest effects are observed.

Apart from frame parameters, however, it is interesting to study the contribution of SLR measurements to the precision of orbital elements, given that SLR observations have been long exploited for the validation of GNSS orbital products. Pertinent results are shown in App. B.

Summary and conclusive remarks

In this chapter, several combination approaches have been implemented in order to evaluate the performance of the space ties on-board GNSS satellites in the derivation of terrestrial reference frames. Such a combination approach automatically overcomes some of the most problematic aspects of the terrestrial ties currently used for the computation of the ITRF, namely their spatial distribution and the frequency of their updates. Nevertheless, the actual strength of the provided link had to be thoroughly investigated in order to characterize the precision and the accuracy of the combined frame parameters.

Despite the conceptual appeal of such a combination approach, results presented in this chapter show that the technique co-locations on-board GNSS satellites cannot substitute terrestrial ties in the derivation of reference frames. The provided link is, in fact, unable to transfer any scale or origin information from one subnetwork to the other. This has been verified at all possible time scales.

At the weekly level, it has been shown that the independent subnetworks maintain their technique-specific frame realization within the "combined" solution. In particular, the GNSS origin and scale do not benefit at all from the combination with LAGEOS observations which provide a more accurate materialization of the frame parameters. Even when zP COs are freed, bringing the estimation of the GNSS scale on the verge of instability, the two technique frames remains distinct and just coexist within the combined adjustment. In this context, it should be considered that the GNSS limited sensitivity to geocenter motion and scale depends on the simultaneous estimation of parameters which are not directly impacted by SLR observations realizing the space ties.

From the perspective of long term frame computation, it has been shown that space ties on-board GNSS satellites are insufficient to synthesize a unique frame from GNSS data and SLR observations to GNSS satellites. When comparisons are attempted with reference to their combined long term solution, in fact, inconsistent results are found for the distinct technique-specific subnetworks.

When compared to the traditional ITRF adjustment, the rigorous space tie approach leads to increased transformation rate uncertainties and to an improved consistency of the independent input solutions with the combined frame. Both these features indicate that the link provided by the investigated space ties is considerably weaker than the one realized by terrestrial ties.

Since currently available SLR observations to GNSS satellites are limited in num-3.6. SUMMARY AND CONCLUSIVE REMARKS 123 ber and considerably noisier than GNSS data and SLR observations to LAGEOS satellites, it has additionally been studied whether significant improvements in the availability and quality of these records would have an impact on the realization of the frame. This study has been performed with simulations providing the formal errors associated to all estimated parameters under different observing scenarios. According to [START_REF] Kuang | Analysis of orbital configurations for geocenter determination with GPS and low-Earth orbiters[END_REF], variations in these uncertainties reflect changes in the observability of the associated quantities. The simulations assumed a perfectly distributed network of high performing stations acquiring a valid return every five minutes during each satellite flyover. Under the best investigated conditions, i.e. when the same tracking precision was implemented for all the involved groups of observations, a precision improvement of 10 -15% was found for T z , while the other frame parameters remained basically unimpacted. The observed effects are rather modest especially considering the technological and infrastructural gap between the state of the art tracking conditions and the implemented scenarios.

In the light of the results presented in this thesis, the co-location of the GNSS and SLR techniques on-board GNSS satellites do not provide a viable tying approach for the realization of terrestrial reference frame. In this context, even improvements in the quality and abundance of the SLR tracking of GNSS satellites are likely to have just a limited impact in the near future. The most relevant contribution of SLR measurements to GNSS vehicles remains the microwave-independent orbit validation process. Any recommendation to GNSS providers concerning retro-reflectors payloads should be issued on the basis of studies aiming at optimizing such a function of the relevant SLR data.

Appendix A Technique-specific adjustmentssupplementary material

In this appendix ancillary information concerning the technique-specific adjustments presented in Chapt. 2 are provided. The following plots illustrate statistics about parameters estimated in the analysis of GNSS data and SLR observations to LAGEOS satellites; the number of station-dependent SLR observations to GNSS satellites is also displayed.

Resolving GNSS initial phase ambiguities is crucial to obtain high quality results. Besides increasing the noise of the final solution, in fact, unresolved ambiguities which need to be set up in the final adjustment may partly absorb geophysical signals and/or impact the GNSS sensitivity to other parameters (e.g. geocenter motion). The percentage of resolved ambiguities for GPS and GLONASS satellites is presented in panel a of Fig. A.1 and A.2 respectively. Results obtained with each different strategy are presented separately. Panel b of the same figures shows the number of parameters initially set up. The increment of GLONASS capable stations within the network is reflected in the increment of the relative ambiguities screened over time. These plots reveal that about 80% of total ambiguities are resolved for GPS satellites, while the equivalent statistics for GLONASS is limited to 40%.

The number of parameters estimated in GNSS analysis is displayed in Fig. A.3 (a). Different groups, i.e. station coordinates, tropospheric delays, orbital elements, EOPs and unresolved initial phase ambiguities, are presented separately. The latters are largely the major contributors to the total number of unknowns in the final adjustment. The relative abundance of the other estimated quantities is shown in So far, simulation results have been presented focusing on the possible benefits that SLR measurements can provide to GNSS-only results. The complementary effects, however, can also be studied. Precision improvements obtained for the SLR realization of frame parameters are presented in Fig. B.7 and B.8. Displayed gains refer to a LAGEOS-only solution computed from the observations collected by a network of 25 stations homogeneously distributed over the globe; an a priori variance factor of 1 cm is assumed. The relevant formal errors are presented in Tab. B.1. The reader interested in a comparison with the corresponding values computed for GNSS-only solution is addressed to Tab. 3.9. It is important to stress, however, that the accuracy and stability of the frame parameters cannot be investigated on the basis of these results. GNSS estimates are characterized by smaller formal errors as a consequence of the greater number of observations, but this does not conflict with the evidence that the quality of the frame origin and scale realization provided by SLR data is actually superior, as discussed in Sec. 1.3.3 and references therein.

C.1 Introduction

The global navigation satellite systems (GNSS) technique provides nowadays a key contribution for monitoring the shape of the Earth's surface and its kinematics and variations on very different spatial and temporal scales. The analysis of the GNSS data allows estimating station positions with 2-to 3-mm repeatability for the horizontal components, and about 4-to 6-mm for the vertical [START_REF] Dach | International gnss service technical report[END_REF]. Station displacements from continuous time series are typically estimated with precisions higher than 1 mm/year [START_REF] Blewitt | Geodetic Observations and Global Reference Frame Contributions to Understanding Sea-Level Rise and Variability[END_REF]. A recent study by [START_REF] Santamaría-Gómez | Correlated errors in GPS position time series: Implications for velocity estimates[END_REF], using a global network of stations and applying upto-date IERS standards, finds an average uncertainty in the order of 0.2 mm/year for the vertical component in the case of time series spanning over a decade or longer. Given these precisions and the dense spatial coverage, GNSS time series contribute to the realization of accurate reference frames [START_REF] Altamimi | an improved solution of the international terrestrial reference frame[END_REF], which are required to support current and future developments in Earth system studies. However, to meet scientific challenges such as monitoring of the longterm variability of global sea level, an improved terrestrial reference frame accurate at a level of 1 mm and stable at 0.1 mm/year is needed [START_REF] Blewitt | Geodetic Observations and Global Reference Frame Contributions to Understanding Sea-Level Rise and Variability[END_REF]. Notwithstanding other requirements, to achieve this objective, the potential of GNSS data series shall be fully exploited.

Almost two decades of continuous GNSS observations have shown that coordinate time series are characterized by linear and non-linear variations as well as by sudden jumps. According to [START_REF] Gazeaux | Detecting offsets in GPS time series: First results from the detection of offsets in GPS experiment[END_REF], these last features can be defined as "sharp change of the mean resulting in a long-lasting effect on parameters, such as velocity estimation". If the discontinuities are not properly accounted for and removed, they might corrupt the reliable modeling of the long-period trends and of non-linear variations likely originated by geophysical and climate changerelated phenomena. There are a number of different causes for such discontinuities. Examples are: earthquakes, changes in the station equipment, antenna mounting problems, multipath, vandalism and data analysis procedure (e.g., change of the reference system).

Several strategies have already been proposed to address the detection of discontinuities in GNSS time series. Different authors have exploited various approaches such as likelihood maximization [START_REF] Picard | A statistical approach for array CGH data analysis[END_REF][START_REF] Bos | Fast error analysis of continuous GPS observations[END_REF], Detection Identification and Adaption (DIA) procedures [START_REF] Perfetti | Detection of station coordinate discontinuities within the Italian GPS Fiducial Network[END_REF][START_REF] Ostini | Fodits: A new tool of the bernese gps software to analyze time series[END_REF][START_REF] Roggero | Discontinuity detection and removal from data time series[END_REF], Bayesian methods [START_REF] Vitti | Sigseg: a tool for the detection of position and velocity discontinuities in geodetic time-series[END_REF], Hidden Markov Models (HMM) [START_REF] Kehagias | Time series segmentation with shifting means hidden markov models[END_REF] and Neyman-Pearson test [START_REF] Teunissen | Testing theory[END_REF]. The Detection of Offsets in GPS Experiment (DOGEx) [START_REF] Gazeaux | Detecting offsets in GPS time series: First results from the detection of offsets in GPS experiment[END_REF] has recently benchmarked these and other methods showing "the statistically significant high performance of hand-picked solutions compared to automated solutions in terms of epoch detection and a posteriori velocity estimation". In fact, of the 25 submitted solutions, only 5 were based on a manual search of the offsets; nevertheless, two of the three best solutions (true positive "TP" greater than 20 %, false positive "FP" and false negative "FN" less than 40% each) were hand-picked.

In the present study, the potentiality of a different method, namely the Sequential t test Analysis of Regime Shifts (STARS) based on the mean, originally proposed by [START_REF] Rodionov | A sequential algorithm for testing climate regime shifts[END_REF]; [START_REF] Rodionov | Application of a sequential regime shift detection method to the Bering Sea ecosystem[END_REF]; [START_REF] Rodionov | Use of prewhitening in climate regime shift detection[END_REF] for climatic studies, has been examined and tested. STARS is a rather simple and fast procedure which relies on only three user-defined parameters of straightforward interpretation. STARS does not require any a priori information on the time of occurrence and on the size of the jumps, as highly desirable given the general unavailability of such data. These characteristics justify the effort of setting up a suitable working procedure allowing the use of the STARS method in the analysis of GNSS data sets.

This study describes the procedure we developed to efficiently apply the STARS algorithm to identify discontinuities in GPS time series. A synthetic data set has been generated to test the method performance. In addition, the procedure has been applied to daily coordinate time series of three GPS EUREF/IGS stations (MSEL, MEDI and BOLG) located in the southeastern Po Plain, in Italy, where results from other observing techniques can provide a sound assessment of the methodology.

C.2 The methodology

For the readers' convenience, the first part of this section outlines the basic features of the STARS method, the complete description can be found in [START_REF] Rodionov | A sequential algorithm for testing climate regime shifts[END_REF]; [START_REF] Rodionov | Application of a sequential regime shift detection method to the Bering Sea ecosystem[END_REF]; [START_REF] Rodionov | Use of prewhitening in climate regime shift detection[END_REF]. The remaining part of the section provides a detailed explanation of the procedure we have developed to properly adapt STARS to GNSS data series.

The algorithm we implemented for this study is based on the version presented by [START_REF] Rodionov | Application of a sequential regime shift detection method to the Bering Sea ecosystem[END_REF]; however, we choose to maintain the minimum length of a regime as inRodionov [2004]. A regime is meant to be the data span delimited by two consecutive shifts identified by the algorithm. STARS makes use of three user-defined parameters: the cutoff length L, representing the minimum time interval between two consecutive discontinuities, the significance level p, of the exploited two-tailed Student t test, and the Huber parameter H [START_REF] Huber | Robust Estimation of a Location Parameter[END_REF] used to compute a weighted mean over the L-day intervals according to [START_REF] Rodionov | Use of prewhitening in climate regime shift detection[END_REF]. The application of the algorithm requires the preliminary evaluation of σ L , which is the mean standard deviation computed over all possible L-day intervals in the time series.

The actual discontinuities identification process is developed in two steps. First, for each point x j of the time series, it is checked whether its distance from the local mean m 1 , evaluated over the previous window [x j-L , xj -1], is significant according to a t test performed at the user-defined significance level p. In other words, x j is found to be a jump candidate if either 140APPENDIX C. IDENTIFYING DISCONTINUITIES IN GNSS TIME SERIES

x j < m 1 -t 2σ 2 L L (case 1) (C.1) or x j > m 1 + t 2σ 2 L L (case 2) (C.2)
where t is the value of the t distribution with 2L -2 degrees of freedom at the given significance level p. The local mean m 1 , is estimated as follows. A mean m 0 , and a standard deviation σ 0 , are computed at first giving equal weight to all data x n (n = j-L, . . . , j-1). Then, a weighted mean of the same data is calculated using weights defined on the basis of the Huber parameter H. Let us define d n = x n -m 0 , namely the data departure from the mean; if d n < Hσ 0 , the corresponding data are given unit weight, otherwise the weight is d -1

n . The estimate of the weighted mean is repeated by updating m 0 and σ 0 to derive the m 1 final value. This procedure is intended to minimize the influence of the noisiest data when computing the mean.

In the second step, jump candidates are accepted if the following inequality holds for all the progressive means x jk , within the L days after day j

x jk < m 1 -t 2σ 2 L L (case 1) (C.3) x jk > m 1 + t 2σ 2 L L (case 2) (C.4)
for k = 1 . . . L, where

x jk = 1 k k-1 l=0 x j+l (C.5)
This allows discriminating between significant discontinuities and local fluctuations. Once a jump is detected, its size is evaluated as the difference between the representative means of consecutive regimes.

STARS was originally developed to analyze annual time series of climatic data. Our work is concerned with long-period GPS series where geophysical information but also several perturbing and undesired signals and/or outliers are present. Therefore, we choose to introduce a pre-processing phase in which outliers are identified and eliminated from the data sets. The adopted cleaning algorithm works by sequentially applying a 2σ ′ and a 3σ ′′ rejection criterion to 6-month portions of the time series, σ ′ being the original variance of the series over the analyzed period, and σ ′′ the correspondent quantity after the first cleaning cycle. Once outliers have been removed, the series are detrended and a mean seasonal cycle is computed by stacking the data. The discontinuities are then detected by applying the STARS algorithm to the detrended and deseasoned series (residuals). The Up, East and North coordinates have been treated independently because jumps might not affect all three position components in the same way.

To tailor the STARS method to GNSS time series, we devised a strategy to properly inspect the algorithm outcomes. This implies a multiple screening of the data sets by adopting four different values of the L parameter, namely L = 10, 15, 20 and 25 days. The choice of an appropriate value for L is of prominent importance. A value too small would lead, in fact, to an overfitting of the series behavior; this is why we disregarded regimes shorter than L, differently from what proposed in [START_REF] Rodionov | Application of a sequential regime shift detection method to the Bering Sea ecosystem[END_REF]. On the other hand, a value too big might produce an artificial misplacement of the detected discontinuities. Moreover, a further constraint on L is given by the estimation process of the jump size. In fact, the difference between local means is only meaningful for those intervals in which the effects of colored noise are non-significant [START_REF] Williams | Offsets in Global Positioning System time series[END_REF]. Provided that abrupt and long lasting changes of the mean should be detected independently of the exploited screening length, it has been decided to retain, among all the detected jumps, only those common to the four cutoff lengths. Common jumps are defined as those discontinuities that are found in the run where L = 25 days, and for which also the tests with shorter time windows identified a jump within ±L run /2, L run being 10, 15 and 20 days. The window L = 25 days has been chosen as reference in the attempt to minimize possible overfitting. Concerning the other parameters, in all cases, we retained p = 0.1 and H = 3.

We complete the discontinuity selection process by imposing a σ L -dependent threshold for the minimum size of an identifiable jump. This is done to account for the noise level of the time series at the reference screening length. Indeed, small jumps (at the millimeter level) can be detected if and only if the extent of data scattering allows it. We choose to make the threshold depend on σ L through an exponentially decaying function. The parameters of this exponential decay are obtained empirically by fitting the curve to a set of points (σ L -threshold pairs, see Fig. C.1) chosen arbitrarily on the basis of visual inspection of several GPS series. It will be shown in the following section that this parameterization fully applies to the selected synthetic case study.

The resulting equation for the threshold, given in percentages of σ L , is

y = A exp σ L b + b 0 (C.6)
where A = 360 ± 2, b = 0.75 ± 0.03 mm and y 0 = 82.2 ± 0.8. For example, if σ L = 3 mm, then the minimum size of the offset is 2.7 mm. Finally, if consecutive jumps, in opposite directions, are detected within a period of three months, the jump acceptance threshold is compared to the difference of the respective magnitudes. The choice of a 3-month interval is motivated by the possible presence of irregular seasonal fluctuations which cannot be completely removed by subtracting the mean seasonal cycle from the data. Contrary to other authors [START_REF] Picard | A statistical approach for array CGH data analysis[END_REF], we did not introduce a rejection criterion based on the jump occurrence frequency. In fact, given the process we adopted to estimate the jump size, limiting the number of detectable offsets per year might lead to an erroneous correction of the time series. This is evident, for example, when abrupt fluctuations, lasting a few months, are present in the time series. In such cases, the algorithm estimates more than one jump to properly characterize the event. As already stated, these jumps are corrected only if the difference of their magnitudes exceeds the jump acceptance threshold. The number of detected jumps, besides depending on L and p, can be somehow controlled by increasing or decreasing the magnitude of the minimum identifiable offset. However, once the parameters have been fixed, the proposed method identifies all those jumps which are statistically significant according to the performed test.

Once the discontinuities have been identified, the coordinate time series are corrected accordingly, and the mean seasonal cycle is then re-estimated. The long-term trend of the corrected and deseasoned time series is finally computed with the CATS software [START_REF] Williams | CATS: GPS coordinate time series analysis software[END_REF] using the maximum likelihood estimation algorithm. The adopted noise model is a combination of variable white and power law with the spectral index estimated.

C.3 Study of synthetic time series

The method performance has been tested against a database of 50 synthetic daily time series over a period ranging from 14 to 18 years. As shown in Eq. C.7, such time series, y(t), have been generated by adding the different components typically recognized in real GNSS data sets, i.e., an intercept a, a linear trend b, a seasonal cycle c(t), a series of jumps j(t), and a combination of white and flicker noise n(t). Missing epochs were also introduced to properly represent realistic conditions. The synthetic series were thus obtained from:

y(t) = a + b + c(t) + j(t) + n(t) (C.7)
The seasonal cycle has been modeled by summing annual and semiannual signals, each of which has been represented as the sum of a sine and a cosine wave of different amplitudes. The values for these amplitudes and for the linear trend parameters have been randomly selected from the database of the global network of stations processed by the JPL (http://sideshow.jpl.nasa.gov/post/series.html). A subset of these stations was also exploited to distribute the missing epochs and to adjust the length of the synthetic time series which were originally created spanning 8,000 days. The missing epochs were introduced exactly where they were found in the real records of long running stations. The resulting time series last between 14 and 18 years. The white component of the noise has been obtained with the MAT-LAB random generator (randn function), while the flicker noise component has been simulated according to [START_REF] Bos | Fast error analysis of continuous GPS observations[END_REF], equations C.3 and C.4. Noise amplitudes were chosen so that the σ L of the deseasoned residuals of the synthetic series ranged between 0.8 and 3.4 mm. Considering that the vertical component is typically two to three times noisier than the horizontals, the selected range of values includes both the low level noise which can be observed for the horizontal components of a well-behaving instrument located in a stable environment, as well as the larger noise associated with the vertical component. In an attempt of properly mimicking a realistic behavior, the highest noise levels have been associated with the seasonal cycle parameters exhibited by the vertical components of the selected real stations. As in [START_REF] Gazeaux | Detecting offsets in GPS time series: First results from the detection of offsets in GPS experiment[END_REF], the time occurrence of jumps follows a binomial distribution; on average, six jumps per series have been simulated, the actual numbers ranging between 1 and 12. Jumps magnitudes have been assigned by sampling a Pareto distribution and by multiplying the obtained value by +1 or -1 with equal probability. The cutoff and the shape parameter of the Pareto distribution have been set to 2 mm and to 1.5, respectively. This means that the magnitude of the smallest simulated offset is 2 mm. The accuracy of the IGS weekly GNSS coordinates is in the order of 3 and 6 mm for the horizontal and vertical components, respectively (http://igs.org/components/prods.html). Given these bandwidths, the reliable identification of 2 mm offsets can appear as a rather optimistic goal to be routinely achieved. However, [START_REF] Gazeaux | Detecting offsets in GPS time series: First results from the detection of offsets in GPS experiment[END_REF] pointed out that, especially for long time series, even the effect of repeated small offsets can significantly impair the estimate of station velocities.

As done for the DOGEx experiment [START_REF] Gazeaux | Detecting offsets in GPS time series: First results from the detection of offsets in GPS experiment[END_REF], three types of events have been considered to assess the quality of the method, namely the TP, the FN, and the FP events. Among the detected jumps, only those which were identified The choice of ±12 days, significantly higher than the ±2 days adopted in the DOGEx experiment, is motivated by the fact that L = 25 days has been selected as the reference run. The screening length L, in fact, sets the minimum distance between two consecutive jumps. As a consequence, if a simulated offset is too close in time to the occurrence of the previously detected jump, the algorithm would be forced to delay detecting such an offset until the adopted time window has elapsed. The outcomes of this test are summarized in If 100% represents the sum of all the possible events, i.e., simulated offsets plus FPs, Fig . C.3a shows that the TPs amount to 48% of the total, the FNs to 24% and the FPs to 28%. If the focus is restricted to the 2/3 of the simulated data set which are meant to correspond to North and East coordinates (i.e., those series with smaller noise level and seasonal cycle amplitudes exhibited by the horizontal components of the exploited real stations), these percentages change to 58, 20 and 22%, respectively. These results are not directly comparable with those obtained in the framework of DOGEx for two main reasons. One point is that in DOGEx the noise is not necessarily time constant at each site and one other aspect is the definition of the TP. Concerning the first issue, the generating function of the timedependent noise is not known. Therefore, even if this behavior was taken into account in our simulated series, this would have likely generated inconsistencies in a one-to-one comparison. One of the key elements of the procedure we devised is the multi-wavelength screening of the time series. This implies that the identification of a TP can occur only within ± half of the longest window used in the analysis which, in our work, is L = 25 days. Therefore, the metric introduced by DOGEx can hardly be exploited in the evaluation of our approach.

Figure C.3b presents the size difference between the estimated and simulated jumps in the case of TP events. A code color has been adopted in the histogram to distinguish between vertical-like (blue) and horizontal-like (gray) components on the basis of their respective noise level and seasonal cycle. The same code color was used in Fig. C.4 which shows the comparison between the estimated and simulated trends. For about half of the horizontal-like components, the agreement is better than ±0.2 mm/year, which is the average uncertainty found by [START_REF] Santamaría-Gómez | Correlated errors in GPS position time series: Implications for velocity estimates[END_REF] for time series of lengths comparable to those analyzed in this work. If the whole data set is taken into account, the percentage of series showing a comparable agreement between simulated and estimated trends reduces to 40% (62% if a bandwidth of ±0.3 mm/year is considered). This is mainly due to the detrimental effect of FP events which are more likely to be detected in the noisiest time series where the data scatter might influence the local mean estimation. Moreover, for particularly noisy series, the stacking of the residuals might be inadequate to properly remove the seasonal cycle. This could leave in the data set fluctuations likely to be interpreted as jumps by the algorithm. Therefore, when processing noisy time series, a manual check of the software outcomes would be a valuable option. This practice has already been suggested by other authors [START_REF] Perfetti | Detection of station coordinate discontinuities within the Italian GPS Fiducial Network[END_REF] and gets further evidence from the findings of the DOGEx experiment. In fact, it has been pointed out by [START_REF] Gazeaux | Detecting offsets in GPS time series: First results from the detection of offsets in GPS experiment[END_REF] that, in general, the manual approach, though not applicable when processing large networks of stations, performs better than fully automated methods tested to A further test, exploiting the prewhitening procedure proposed by [START_REF] Rodionov | Use of prewhitening in climate regime shift detection[END_REF], has been carried out on the noisiest time series. The prewhitening routine was added by Rodionov to the original algorithm to remove from the time series the fingerprint of colored noise. In fact, the STARS method is based on a t test working rigorously for Gaussian processes only. The results show that there is no improvement using the prewhitening procedure for the data set under study. This could be explained by the negligible effect of the colored noise over the selected time windows [START_REF] Williams | Offsets in Global Positioning System time series[END_REF].

The outcomes of the analysis of the synthetic database were also useful to identify possible causes of FN events. Among the undetected jumps, some were those originally simulated with magnitude smaller than the detection threshold estimated for the specific time series. Other FN occurrences were due to underestimation of the magnitudes, thus causing an erroneous rejection of the offsets by the method. Moreover, due to the minimum length, L, imposed to a regime, the runs at different screening wavelength might detect the same jump at slightly distinct epochs, following the definition outlined in Sect. C.2. In this regard, it has been tested that disregarding the minimum regime length [START_REF] Rodionov | Application of a sequential regime shift detection method to the Bering Sea ecosystem[END_REF] leads to an increased number of FPs, as many more discontinuities are detected by the algorithm. Therefore, the chosen set up for the proposed strategy seems to offer the best possible trade-off between FN and FP events. The Bernese software version 5.0 [START_REF] Dach | Manual of Bernese GPS Software Version 5.0[END_REF] was used to process the MSEL and BOLG GPS data by adopting high-accuracy International GNSS Service (IGS) products (i.e., ionospheric files, satellite orbits, and Earth rotation parameters). Five IGS stations in Europe were chosen as fiducial sites in the network adjustment procedure. The results are given in the IGS08; the coordinates derived prior to the reference system change (from IGS05 to IGS08) have been converted by means of the antenna calibration updates and the recommended Helmert transformation parameters [START_REF] Rebischung | IGS08: the IGS realization of ITRF[END_REF]. For MEDI, the coordinate time series were downloaded from the JPL archive (http://sideshow.jpl.nasa.gov/post/ series.html).

C.4 The area and the stations

C.5 Results

The methodology, as described previously, has been applied to a number of GPS coordinate time series affected by discontinuities of various nature. For each of the analyzed stations, the residual North, East and Up coordinates are presented in C.1 lists the estimated long-period trends before and after the corrections have been applied. For each trend, the estimated spectral index of the power law noise is listed. A sound explanation has been found for many of the different offsets marked on the plots. For example, as described in the previous section, the BOLG Up component shows repeated and sudden upward winter jumps, clearly visible in the years 2005, 2006, 2007, 2009 and 2011. Resettings of the antenna pole took place in May 2009 and May 2012 and are clearly recognizable in the time series. A good clue for the capability of the method to properly correct the data set is given in Fig. C.8b, where the difference between the original Up time series and the corrected one drops to zero after the first repositioning of the antenna. By contrast, a slight misalignment between the two series can be detected following the second resetting. This is due to the occurrence, at the time of the resetting, of a Mw 5.9 earthquake (May 20th, 2012), with epicenter about 50 km from the station. The earthquake signature is observed in the Up and also in the North component, in agreement with the north-south compression mechanism of this seismic event [START_REF] Piccinini | Source complexity of the May 20, 2012, Mw 5.9, Ferrara (Italy) event[END_REF][START_REF] Cesca | Source modelling of the m5-6 emilia-romagna, italy, earthquakes (2012 may 20-29)[END_REF]. Tests and replacement of the antenna cable, performed in March 2005, have also been identified in the North coordinate. For the East, the discontinuity due to the change from IGS05 to IGS08 is properly detected on April 17th, 2011.

The comparison between the MSEL and MEDI time series (Figs. C.9 C.10) shows that MEDI is characterized by a larger noise level in all three components; the standard deviation of the MEDI Up residuals is about 1.5 times larger than that of MSEL, while a factor of 2 is found for the North and East components. The MSEL Up presents an anomalous bump lasting a few years, from about 2000 to 2004, in conjunction with cycle 23 of the solar flux, which was characterized by two very strong maxima at the beginning of 2000 and towards the end of 2001. The receiver firmware was upgraded only in February 2004 and jumps identified in this period might be associated with this forcing (Fig. C.9a). Other examples of detected discontinuities are those associated with the May 20th, 2012, earthquake mentioned above. The jump was recognized in both MSEL and MEDI North components (Figs. C.9a C.10a), while it was detected only in the MSEL Up series because of the higher noise level of MEDI Up. The antenna of the MEDI station was cleaned at the end of March 2009 since, during the years, guano had accumulated within the concentric cylinders of the choke ring design. The guano-originated noise seems to affect both the Up and East coordinates during the period 2006-2009. It shall be pointed out that in those periods where the MEDI noise level becomes significantly higher than the σ L value (from January 2006 to about March 2009 for East and Up; the whole year 1999 for North), we have rescaled the jump acceptance threshold according to Eq. C.6. In this case, σ L represents the average standard deviation of the L-day intervals evaluated over the period of increased noise.

The JPL makes available (http://sideshow.jpl.nasa.gov/post/tables/table3. html) the discontinuities identified in the analyzed series of station coordinates by using analyze tseri [START_REF] Dong | Estimating regional deformation from a combination of space and terrestrial geodetic data[END_REF], included in the QOCA package (http: //qoca.jpl.nasa.gov). In the JPL strategy, a jump always affects all three coordinates, while we treat the three components independently. This poses a problem in a direct comparison of the jumps identified by the two approaches. Concerning the MEDI station, over the analyzed time period, JPL finds three jumps (October 1996, April 2006, March 2009). The only common jump is the one we detected in 1996 in the North component. The JPL strategy does not identify the May 20th, 2012, earthquake signature. Additionally, the two jumps of 2006 and 2009 have a major impact on the estimate of the Up linear trend, which turns out to be significantly larger (higher subsidence) than that estimated either with our procedure or using VLBI data (see Table C.1).

The trends of the MSEL and MEDI Up series prior and after corrections compare within the errors. There are significant differences in the original horizontal velocity components (∆ Ṅ = 0.59 ± 0.22 mm/year, ∆ Ė = -1.07 ± 0.16 mm/year), despite the fact that the two stations are very close to each other. Such differences are reduced, in particular in the East component by about a factor of 2 (∆ Ṅ = 0.45±0.13 mm/year, ∆ Ė = -0.48 ± 0.14 mm/year), after removing the identified discontinuities and the re-estimated mean seasonal cycle (Table C.1). The remaining discrepancies suggest the presence of intrasite motions in the Medicina observatory area, which were also detected by a few ground surveys carried out with a total station [Sarti et al., 2013a]. For a quantitative comparison of the GPS measurements with the terrestrial observations, we estimated the difference between the MEDI and MSEL 3D velocities over the common period 2001-2007. This velocity vector difference has also been estimated without removing the mean seasonal cycle from the GPS data. The results (Table C.2) agree within the respective uncertainties.

The trends have also been compared with those of the Medicina VLBI data. The VLBI solutions were kindly provided by D. Macmillan (2013, private communication). The only correction we applied to these coordinate series is that for the displacements caused by the May 20th, 2012, earthquake. The jump magnitudes were also provided by D. Macmillan (2013, private communication). The estimated linear trends and the associated errors are displayed in Table C.1. The seasonal signal was not estimated and removed from the VLBI time series because of the relatively uneven data sampling. The corrected VLBI and GPS coordinate trends are in satisfactory agreement; in particular, for the Up component, the estimated velocities agree within the relevant errors, while for the horizontal trends the applied corrections lead to an improved comparison between the GPS and VLBI estimates. The small differences between the results of the GPS and VLBI techniques might arise from the foundation depth of the monuments. The VLBI pilings are deeply founded on layers of sand and gravelly sand, while both GPS antennas are mounted on shallower pillars resting on over-consolidated plastic clays and silty clays [START_REF] Romagnoli | Influence of soil consolidation and thermal expansion effects on height and gravity variations[END_REF]Sarti et al., 2013a]. Differences in the technical realization of the two GPS pillars might also contribute to the observed signals.

C.6 Conclusions

Discontinuities in GPS coordinate time series are a common problem and arise from a number of different causes. Given the detrimental effect they have on the reliable estimation of long-term behaviors, it is of prominent importance to properly clean the data sets. This is also a key factor when attempting a geophysical interpretation of the time series, where the sought signals might be faint and hard to separate from the data noise.

We have implemented an analysis procedure, based on the STARS algorithm, which has proved to be well suited for the identification and size estimation of discontinuities in the GPS coordinate time series. Advantages of this approach are the rather simple conceptual background, the independence from a priori information on the epoch and magnitude of jumps, the easy interpretation of the user-defined parameters and the fast execution time. Differently from other approaches [START_REF] Ostini | Fodits: A new tool of the bernese gps software to analyze time series[END_REF][START_REF] Santamaría-Gómez | Correlated errors in GPS position time series: Implications for velocity estimates[END_REF] which remove seasonal oscillations by means of mathematical functions such as sinusoids, our analysis uses the mean seasonal cycle estimated from the data. We have developed an automated procedure, which helps when processing a large number of stations. However, when the dimension of the data set allows it, a manual check of the software outcomes would be a valuable option. Actually, it has been pointed out by [START_REF] Gazeaux | Detecting offsets in GPS time series: First results from the detection of offsets in GPS experiment[END_REF] that the manual approach, in general, performs better than existing automated or semi-automated methods.

To assess the capability of the STARS approach and of the implemented procedure, we have created a synthetic database reproducing realistic GNSS time series featuring a linear trend, seasonal cycle, discontinuities, data gaps and white plus flicker noise. This analysis has shown that the methodology allowed identifying 58% of TP events when the σ L noise level of the residual series, spanning almost two decades, ranges between 0.8 and 2.2 mm. If the whole data set is considered, this means including also the higher noise series, the percentage of the TPs decreases to 48%. The differences between the estimated and simulated linear trends agree within ±0.2 mm/year in the 40% of the total cases, this statistics improves to 62% if a bandwidth of ± 0.3 mm/year is considered.

It should be pointed out that these encouraging results can hardly be assessed within the metrics introduced by the DOGEx, where the best solutions were required to obtain a TP percentage larger than 20% and FP and FN statistics lower than 40% each. A direct comparison between the results of our case study and the DOGEx findings is, in fact, problematic due to the intrinsic nature of the procedure that we have developed and to the differences in the generation of the series noise.

The procedure has been applied to real data sets, namely three GPS stations, BOLG, MSEL and MEDI, located in the southeastern Po Plain, in Italy. The BOLG GPS data series were severely affected by a number of discontinuities mainly due to the antenna installation, cable replacement and earthquake occurrence. Since the epoch of these events was known, testing our procedure over this series was a reliable proof of the method quality.

At Medicina, the analysis has identified only one discontinuity common to both GPS stations, namely the jump in the North component in conjunction with the May 20th, 2012, earthquake. Most of the other jumps, instead, are instrument dependent, which means that they resulted from different antenna and/or receiver behaviors as well as from monument stability and its environment. The removal of all the detected discontinuities has enhanced the agreement between the coordinate trends of the two stations, in particular, those of the horizontal components (Table C.1). A further comparison has been performed between the difference of the MEDI and MSEL 3D velocities, as derived from the corrected GPS coordinates, and the baseline change from terrestrial geodetic measurements [Sarti et al., 2013a]. The results agree within the respective uncertainties. The availability, at the observatory site, of VLBI observations made it possible to perform an additional check for assessing the capability of the methodology to properly correct the time series. In particular, for the trends of the horizontal components, after removing the discontinuities, the agreement with the VLBI data improves up to 0.4 mm/year in the case of MEDI East. The remaining small differences between the results of the GPS and VLBI techniques might arise from the foundation depth of the monuments.

χ 2 i

 2 DOF values computed from the weekly/9-day adjustments of a) GNSS, b) SLR-to-GNSS and c) SLR-to-LAGEOS observations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 3.3 RMS of the observation residuals computed for each baseline including either GODZ or ALIC stations. . . . . . . . . . . . . . . . . . . . 100 3.4 Translation and scale offsets estimated between the frames realized, on a weekly basis, by a GNSS-only solution and by the GNSS subnetwork of the fully combined adjustment. . . . . . . . . . . . . . . . 102 3.5 Scale difference between fully-combined and GNSS-only weekly frames both estimated with free zP COs. . . . . . . . . . . . . . . . . . . . . 104 3.6 Translation and scale offsets estimated between the frames realized, on a weekly basis, by a LAGEOS-only solution and by the SLR subnetwork of the fully combined adjustment -CASE A. . . . . . . . . . 104 3.7 Translation and scale offsets estimated between the frames realized, on a weekly basis, by a LAGEOS-only solution and by the SLR subnetwork of the fully combined adjustment -CASE B. . . . . . . . . . 105 3.8 Sketch representation of inter-technique combinations with/without applying effective tying constraints. . . . . . . . . . . . . . . . . . . . 106 3.9 GPS and GLONASS range biases. . . . . . . . . . . . . . . . . . . . . 109 3.10 Relative improvement in the formal errors of reference frame parameters with respect to a GPS-only solution. . . . . . . . . . . . . . . . 119 3.11 Relative improvement in the formal errors of reference frame parameters with respect to a GPS+GLONASS solution. . . . . . . . . . . . 119 3.12 Relative improvement in the formal errors of reference frame parameters with respect to a GPS-only solution encompassing the estimation of zP COs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 A.1 Assessment of ambiguity resolution strategies for GPS observations. . 126 A.2 Assessment of ambiguity resolution strategies for GLONASS observations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 A.3 Number of parameters set up in GNSS analysis. . . . . . . . . . . . . 127 A.4 Number of parameters set up in SLR analysis of LAGEOS observations.127 A.5 Number of SLR observations to GNSS satellites collected daily at each site. Stations are sorted according to their ID code. The color code indicates the number of different targeted satellites. . . . . . . . 129 B.1 Relative improvement in the formal errors of orbital parameters (initial state vector) with respect to a GPS-only solution. . . . . . . . . . 132 B.2 Relative improvement in the formal errors of empirical orbit parameters with respect to a GPS-only solution . . . . . . . . . . . . . . . . B.3 Relative improvement in the formal errors of orbital parameters (initial state vector) with respect to a GPS+GLONASS solution. . . . . B.4 Relative improvement in the formal errors of empirical orbit parameters with respect to a GPS+GLONASS solution . . . . . . . . . . . B.5 Relative improvement in the formal errors of orbital parameters (initial state vector) with respect to a GPS-only solution. . . . . . . . . . B.6 Relative improvement in the formal errors of empirical orbit parameters with respect to a GPS-only solution encompassing the estimation of zP COs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.7 Relative improvement in the formal errors of reference frame parameters with respect to a LAGEOS-only solution with 25 stations and sigma=1 cm. The characteristics of displayed cases are detailed in Tab. 3.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.8 Relative improvement in the formal errors of reference frame parameters with respect to a LAGEOS-only solution with 25 stations and sigma=1 cm. The characteristics of displayed cases are detailed in Tab. 3.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.1 Empirical relation between GNSS time series noise level and minimum size of identifiable discontinuities. . . . . . . . . . . . . . . . . . . . . C.2 Example of synthetic time series . . . . . . . . . . . . . . . . . . . . . C.3 Assessment of the method performance . . . . . . . . . . . . . . . . . C.4 Comparison between the originally simulated trends and the velocities of the corrected time series. . . . . . . . . . . . . . . . . . . . . . . . C.5 Case study area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.6 Medicina site map . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.7 BOLG station -Technical drawing of the antenna mount . . . . . . . C.8 Results -BOLG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.9 Results -MSEL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.10 Results -MEDI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  of the WRMS of nine-day GNSS station position residuals estimated during the long term stacking of the time series. . . . . 2.3 Translation and scale offsets estimated between the secular frames computed from GNSS 9-day solutions either setting up or ignoring transformation parameters during the long term stacking. . . . . . . . 2.4 Mean GNSS origin and scale offsets with reference to ITRF2008. . . 2.5 Average values of the time series displayed in Fig. 2.14. . . . . . . . . 2.6 Mean SLR origin and scale offsets with reference to ITRF2008. . . . . 3.1 A priori σ of unit weight. . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Translation T, scale D and orientation R offsets estimated from GNSS/GNSS + SLRtoGNSS frame with respect to the SLRtoLA-GEOS frame. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Translation Ṫ, scale Ḋ and orientation Ṙ rates estimated from GNSS/GNSS + SLRtoGNSS frame with respect to the SLRtoLAGEOS frame. . . . 3.4 WRMS of station position and velocity residuals estimated between the fully combined frame and each input solution. . . . . . . . . . . . 3.5 Translation T, scale D and orientation R offsets estimated between the fully combined secular frame and the technique-specific subnetworks within the GNSS+SLRtoGNSS long term stacking. Relevant formal errors are indicated within brackets. . . . . . . . . . . . . . . . 3.6 General options implemented in the simulations . . . . . . . . . . . . 3.7 List of simulated cases investigated in this study. . . . . . . . . . . . 3.8 Parametrization adopted for the quantities set up in the simulations. 3.9 Formal errors of the physical frame parameters computed for simulations based exclusively on GNSS data. . . . . . . . . . . . . . . . . . B.1 Formal errors of the physical frame parameters resulting from a simulated LAGEOS-only solution based the SLR observations collected from a network of 25 stations. . . . . . . . . . . . . . . . . . . . . . . C.1 Coordinate trends computed before and after removing the identified jumps. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.2 Difference between the MEDI and MSEL 3D velocities derived from GPS observations and estimate of the baseline change from terrestrial measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156

Figure 1

 1 Figure 1.1 -Visual representation of the seven transformation parameters needed to transform one reference frame into another. From Pollet [2011, chapter 2]

Figure 1

 1 Figure 1.2 -Map of ITRF co-located ground sites active in 2014, courtesy of Z. Altamimi. Yellow stars, red dots, blue circles and green square represent VLBI, GPS/GNSS, DORIS and SLR stations respectively.

Figure 1

 1 Figure 1.3 -Visual representation of space and terrestrial tie vectors -adapted from Thaller et al. [2015]. Ties connect technique-specific reference points either (a) onboard the satellite or (b) at ground co-located stations. As the link in space is realized by the joint estimation of common orbit parameters, space ties relate the reference points of the relevant payloads via the satellite center of mass.

Figure 2

 2 Figure 2.1 -Representation of GNSS carrier phase modulation on the basis of a superimposed code: (a) original carrier, (b) modulating code, (c) resulting signal.In the case of GPS two codes are employed. The C/A (clear acquisition) code has a frequency of 1/10 of the fundamental one maintained by the clocks and is modulated on L1 only; the P (precise) code, instead, provides a ten times higher accuracy and is superimposed to both L1 and L2. The navigation message, identifying the satellite that broadcast the signal, its health status and other operative information, is also modulated on both carriers.

Figure 2

 2 Figure2.2 -Scheme of GNSS receivers working principle. The green dot indicates the phase of the signal at time t -τ , which is maintained by the wave broadcast at that epoch. The signal copy locally generated within the receiver, instead, continues to evolve during the flight time τ , so that, at the acquisition epoch t the cumulative phase difference consists of a certain number of complete oscillations (a) plus a subinteger discrepancy (b). Modulations are ignored in this sketch for the sake of clarity.

Figure 2

 2 Figure2.3 -Scheme of the satellite-Sun oriented frame: the D -axis is the satellite-Sun direction (i.e. it is normal to the plane in which the solar panels lie), the Y -axis is the rotation axis of the solar panels and the B -axis completes the right-handed trihedron.

Figure 2

 2 Figure 2.6 -WRMS of station position residuals resulting from a comparison between daily (blue line) and three-day (red-line) solutions with the IGS repro2 combined results. Horizontal lines indicate the upper limit of the general agreement shown by the solutions provided by the official IGS ACs. Long-arc estimates clearly improve the comparison; the coarse results obtained in the second half of year 2014 correspond to a period of data shortage of MAUI station. The network hole left in an area which was already suffering from a rather sparse coverage resulted in a degradation of the estimates for sites in the Pacific region.

Figure 2

 2 Figure 2.7 -Number of observations (blue line -left axis) and of estimated parameters (red line -right axis) for each nine-day solution. The increment of GLONASS capable systems is reflected in the trends shown by both curves.

Figure 2

 2 Figure 2.8 -WRMS of GNSS station position residuals estimated between each nine-day solution and the derived long term frame. During the stacking, the seven transformation parameters are set-up in case a) and ignored in case b).

Figure 2

 2 Figure2.9 -Translation and scale offsets computed during the long term stacking of GNSS 9-day solutions. MC have been imposed with reference to ITRF2008.

Figure 2 Figure 2 Figure 2

 222 Figure 2.12 -Number of observations (blue line -left axis) and of estimated parameters (red line -right axis) for each weekly SLRtoLAGEOS solution.

Figure 2 .

 2 Figure 2.15 -WRMS of station position residuals resulting from a comparison between weekly SLRtoLAGEOS solutions and ILRS combined results. Computations have been carried out considering all stations (blue line), or core sites only (red line).

Figure 2 . 16 -

 216 Figure 2.16 -Translation and scale offsets computed during the long term stacking of SLRtoLAGEOS weekly solutions. MC have been imposed with reference to ITRF2008.

Figure 2

 2 Figure 2.17 -Number of SLR observations to GNSS satellites for each 9-day solution. A clear acquisition increment is observed starting from August 2014 coinciding with the beginning of the first LARGE experiment.

Figure 2 .

 2 Figure 2.18 -Total number of SLR observations to GNSS satellites collected by each ILRS site over the period 2011 -2014. Stations are sorted in descending order of acquired data.

Figure 2

 2 Figure 2.19 -WRMS of station position residuals resulting from a comparison between 9-day estimates based on SLR observations to GNSS satellites and weekly ILRS combined solutions derived from LAGEOS data.

Figure 2

 2 Figure2.20 -Comparison of translation offsets computed during the long term stacking of GNSS (blue) and SLRtoLAGEOS (red) solutions with MC imposed with reference to ITRF2008. To facilitate the comparison, both time series have been detrended and smoothed with a Vondrak filter with a 7 cpy cutoff frequency[START_REF] Vondrak | A contribution to the problem of smoothing observational data[END_REF].

  Figure 3.2 -Time series of χ 2 i DOF values computed from the weekly/9-day adjustments of a) GNSS, b) SLR-to-GNSS (trial solution) and c) SLR-to-LAGEOS observations. Average values are reported in the displayed text boxes.

Figure 3

 3 Figure 3.3 -RMS of the observation residuals computed for each baseline including either GODZ or ALIC stations.

χ 2 i

 2 DOF have been used.

Figure 3 . 4 -

 34 Figure 3.4 -Translation and scale offsets estimated between the frames realized, on a weekly basis, by a GNSS-only solution and by the GNSS subnetwork of the fully combined adjustment.

Figure 3

 3 Figure 3.5 -Scale difference between fully-combined and GNSS-only weekly frames both estimated with free zP COs.

Figure 3 . 6 -

 36 Figure 3.6 -Translation and scale offsets estimated between the frames realized, on a weekly basis, by a LAGEOS-only solution and by the SLR subnetwork of the fully combined adjustment -CASE A.

Figure 3

 3 Figure 3.7 -Translation and scale offsets estimated between the frames realized, on a weekly basis, by a LAGEOS-only solution and by the SLR subnetwork of the fully combined adjustment -CASE B.

Figure 3

 3 Figure 3.8 -1D-sketch representation of inter-technique combinations with/without applying effective tying constraints. (a) Physical situation: two stations are located at a certain distance from the center of mass. (b) In each system-specific frame, the actual station coordinates depend on the inherent realization of the frame origin and scale.(c) When employed ties fail to actually link the two frames, stations keep positioning, in the "combined" solution, according to their own intrinsic frame realization. (d) In properly combined frame, instead, all subnetwork adopt the same origin and scale and spatial co-location is restored. The author is grateful to P. Rebischung for providing the original idea for this representation.

Figure 3

 3 Figure 3.9 -Range biases estimated from the long term stacking of fully combined weekly frames derived from CASE B solutions, see Sec. 3.3. RBs are estimated per station and group of satellites, results concerning the two GPS spacecrafts equipped with retro-reflector arrays are presented in panel a, while the analogous results for all GLONASS vehicles are shown in panel b. Error bars are consistently estimated during the stacking.

  Fig. A.3 (b) and reflects the parametrization choices detailed in Sec. 2.2.2. The equivalent plots relative to SLR analysis are displayed in Fig.A.4.

Figure

  Figure A.3 -Panlel a): number of parameters estimated in GNSS analysis; panel b): pie chart displaying the relative abundance of set-up parameters excluding initialphase ambiguities.

Figure A. 4 -

 4 Figure A.4 -Panlel a): number of parameters estimated in the SLR analysis of LAGEOS observations; panel b): pie chart displaying the relative abundance of set-up parameters.

Figure A. 5 -

 5 Figure A.5 -Number of SLR observations to GNSS satellites collected daily at each site. Stations are sorted according to their ID code. The color code indicates the number of different targeted satellites.

  The very high gains shown in the plots are a direct consequence of the large number of observations added with respect to the LAGEOS-only solution. Formal errors of the physical frame parameters resulting from a simulated LAGEOS-only solution based the SLR observations collected from a network of 25 stations.

Figure C. 1 -

 1 Figure C.1 -The exponential decay function describing the empirical relation between the noise level, over the 25-day interval, and the minimum size of identifiable discontinuities. The function has been derived empirically by fitting a set of arbitrarily chosen σ L -threshold pairs

  Figure C.2 -Example of a generated synthetic time series (upper panel). The lower panel shows the residual series and the occurrence of FP (cyan line), FN (light magenta line) and TP (green line) events

  Fig. C.3.

Figure

  Figure C.3 -Assessment of the method performance. a Frequency of occurrence of true-positive (TP), false-positive (FP) and false-negative (FN) events. b Size difference between the estimated and simulated jumps in the case of TP events. Offsets found in vertical-like time series are shown in blue, while those found in horizontal-like time series are shown in gray

  146APPENDIX C. IDENTIFYING DISCONTINUITIES IN GNSS TIME SERIES date.

Figure C. 4 -

 4 Figure C.4 -Comparison between the originally simulated trends and the velocities of the corrected time series. Vertical-like time series are shown in blue, while horizontal-like data sets are shown in gray

  The case study area is in northeastern Italy and it encompasses the city of Bologna and the nearby smaller town of Medicina (Fig. C.5). Since many years, the Department of Physics and Astronomy (DIFA) of the University of Bologna runs several permanent GPS stations, two of which are located at Bologna (BOLG) and Medicina (MSEL) in the Po River Valley. In addition, at the Medicina Radioastronomy Observatory, the Italian Space Agency (ASI) operates another GPS system (MEDI). These stations are characterized by long and continuous data records, MEDI and MSEL started acquiring data in March and July 1996, respectively, while BOLG initiated in May 1999. The BOLG and MSEL sites belong to the EUREF EPN (European Permanent Network), while MEDI is part of the IGS (International GNSS Service) network and is also a reference station for the International Terrestrial Reference Frame (ITRF). The MSEL and MEDI stations were chosen for this study because of their vicinity, about 30 m (Fig. C.6), which should entail common environmental and tectonic forcings. Moreover, there is also a VLBI antenna in close proximity of the two GPS sites, which provides the opportunity for a valuable inter-technique comparison.

Figure

  Figure C.5 -Case study area: southeastern Po Plain, Italy. Geographical location of the city of Bologna and of the nearby town of Medicina where the GPS stations are installed

Figure

  Figure C.7 -BOLG station. Technical drawing of the antenna mount

  Figs. C.8a, C.9a and C.10a. In these plots, vertical red lines indicate the epochs at which discontinuities have been detected. In addition, green lines indicate the epoch of hardware replacement, firmware upgrade and change of the reference system. Figures C.8b, C.9b and C.10b show the difference between the original and the corrected series. Table

Figure C. 8 -

 8 Figure C.8 -BOLG station: a Up, East and North residuals of the original GPS coordinate time series. Red lines indicate the detected discontinuities; green lines indicate cable replacement and reference system change. b Difference between the original and the corrected Up, East and North coordinate time series; the scales are different from those of panel a to clearly distinguish the magnitude of the discontinuities

  1.1 Visual representation of the seven transformation parameters needed to transform one reference frame into another. From Pollet [2011, chapter 2] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Map of ITRF co-located ground sites active in 2014. . . . . . . . . . . 1.3 Visual representation of space tie and terrestrial tie vectors -adapted from [Thaller et al., 2015] . . . . . . . . . . . . . . . . . . . . . . . . Number of available observations and of estimated parameters for each nine-day GNSS solution. . . . . . . . . . . . . . . . . . . . . . . 2.8 WRMS of GNSS station position residuals estimated between each nine-day solution and the derived long term frame . . . . . . . . . . . 2.9 Translation and scale offsets computed during the long term stacking of GNSS 9-day solutions. MC have been imposed with reference to ITRF2008. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Scheme of the optical path within corner-cube retro-reflectors. . . . . 2.11 Simplified scheme of SLR station working principle. . . . . . . . . . . 2.12 Number of available observations and of estimated parameters for each weekly SLRtoLAGEOS solution. . . . . . . . . . . . . . . . . . . 2.13 Map of SLR core stations available to this study. . . . . . . . . . . . 2.14 WRMSs of SLR station position residuals estimated between each weekly solution and the derived long term frame . . . . . . . . . . . . 2.15 Comparison of weekly SLRtoLAGEOS solutions with ILRS combined results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.16 Translation and scale offsets computed during the long term stacking of SLRtoLAGEOS weekly solutions. MC have been imposed with LIST OF FIGURES 2.19 WRMS of station position residuals resulting from a comparison between 9-day estimates based on SLR observations to GNSS satellites and weekly ILRS combined solutions derived from LAGEOS data. . . 90 2.20 Comparison of translation offsets computed for GNSS (blue) and SLRtoLAGEOS (red) solutions with MC imposed with reference to ITRF2008. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 3.1 General combination scheme. . . . . . . . . . . . . . . . . . . . . . . 97 3.2 Time series of

2.1 Scheme of GNSS carrier phase modulation. . . . . . . . . . . . . . . . 2.2 Scheme of GNSS receivers working principle. . . . . . . . . . . . . . . 2.3 Scheme of the satellite-Sun oriented frame . . . . . . . . . . . . . . . 2.4 Network map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 GLONASS network map . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Comparison of daily and three-day solutions with IGS repro2 combined results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 reference to ITRF2008. . . . . . . . . . . . . . . . . . . . . . . . . . . 2.17 Number of SLR observations to GNSS satellites for each 9-day solution. 2.18 Total number of SLR observations to GNSS satellites collected by each ILRS station. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  ].

	8 BKG Bundesamt fűr Kartographie und Geodasie
	ASI Agenzia Spaziale Italiana
	ESA European Space Agency
	GFZ GeoForschungs Zentrum

Table 3 .

 3 2 -Translation T, scale D and orientation R offsets estimated from GNSS/GNSS+SLRtoGNSS frame with respect to the SLRtoLAGEOS frame. Relevant formal errors are indicated within brackets.

		Ṫ1	Ṫ2	Ṫ3	Ḋ	Ṙ1	Ṙ2	Ṙ3
		[mm/y] [mm/y] [mm/y] [ppb/y] [mas/y] [mas/y] [mas/y]
	Terrestrial ties	0.2	0.0	-0.6	-0.04	0.005	-0.005	0.015
		(0.2)	(0.2)	(0.2)	(0.03) (0.005) (0.005) (0.006)
	Space ties	0.2	0.2	-0.6	0.0	-0.04	0.006	0.027
	EVC	(0.3)	(0.3)	(0.3)	(0.04) (0.044) (0.032) (0.043)
	Space ties	0.5	0.5	-2.2	-0.15	-0.006	0.028	-0.008
	NO EVC	(1.0)	(0.9)	(0.8)	(0.09) (0.069) (0.039) (0.071)

Table 3 .

 3 3 -Translation Ṫ, scale Ḋ and orientation Ṙ rates estimated from GNSS/GNSS+SLRtoGNSS frame with respect to the SLRtoLAGEOS frame. Relevant formal errors are indicated within brackets. Table3.4 -WRMS of station position and velocity residuals estimated between the fully combined frame and each input solution.

		Pos Res WRMS		Vel Res WRMS
		E	N	U	E	N	U
		[mm] [mm] [mm]	[mm/y] [mm/y] [mm/y]
	Terrestrial ties
	GNSS	0.0	0.0	0.1	0.0	0.1	0.2
	SLR	2.1	1.8	0.3	1.7	1.8	0.4
	TIES	2.9	3.2	6.5
	Space ties -EVC
	GNSS	0.0	0.0	0.0	0.0	0.0	0.0
	SLRtoGNSS	6.6	9.1	10.8	0.4	0.5	0.6
	SLRtoLAG	0.8	0.7	0.1	0.7	0.5	0.1
	Space ties -NO EVC
	GNSS	0.0	0.0	0.0	0.0	0.0	0.0
	SLRtoGNSS	6.7	8.8	11.1	0.3	0.2	0.9
	SLRtoLAG	0.0	0.1	0.0	0.0	0.0	0.0

Table 3 .

 3 7 -List of simulated cases investigated in this study. The table is divided in four blocks which group simulations assuming the same values of σ 0 for SLR observations to LAGEOS and GNSS satellites.

	].

  TableC.1 -Coordinate trends computed before and after removing the identified jumps. For each trend, the estimated spectral index (SI) of the power law noise is listed.The VLBI results have been corrected only for the May 20th, 2012, earthquake offsets as provided by D. MacmillanMEDI-MSEL baseline change (mm/year)Original Corrected GPS without seasonal cycle 1.1 ± 0.7 0.6 ± 0.4 GPS withseasonal cycle 1.2 ± 0.7 0.8 ± 0.5 Ground surveys[Sarti et al., 2013a] 0.7 ± 0.2 TableC.2 -Difference between the MEDI and MSEL 3D velocities derived from GPS observations and estimate of the baseline change from terrestrial measurements over the 2001-2007 time period. Both the original and corrected baseline change estimates are listed for the MSEL and MEDI stations. The first row displays the result obtained by removing the mean seasonal cycle from the GPS time series, while the second line presents the result including the mean seasonal cycle

	HMM	Hidden Markov Model
	IAG	International Association of Geodesy
	IAU	International Astronomical Union
	IC	Internal Constraints
	IDS	Trend International DORIS Service BOLG	SI	MEDI	SI	MSEL	SI VLBI Medicina SI
	Original data (mm/year) IERS International Earth Rotation and Reference Systems Service
	IGS	U International GNSS Service -7.51 ± 1.16 1.3 -1.85 ± 0.17 0.8 -1.87 ± 0.11 0.9	-1.85 ± 0.18	0.7
	ILRS	E International Laser Ranging Service 21.55 ± 0.16 1.2 22.77 ± 0.14 1.0 21.70± 0.08 1.1	22.03 ± 0.07	1.1
	ITRF	N International Terrestrial Reference Frame 19.18 ± 0.24 1.3 17.31 ± 0.22 1.3 17.90 ± 0.05 1.0	17.67 ± 0.03	0.2
	Corrected data (mm/year) ITRS International Terrestrial Reference System
	IUGG	U International Union of Geodesy and Geophysics -7.89 ± 0.47 1.1 -2.12 ± 0.17 0.8 -2.00 ± 0.08 0.7	-1.92 ± 0.16	0.7
	IVS	E International VLBI Service 21.52 ± 0.11 1.1 22.34 ± 0.13 0.9 21.86 ± 0.06 1.0	22.04 ± 0.07	1.1
	LAGEOS	N Laser Geodynamics Satellites 19.18 ± 0.24 1.1 17.47 ± 0.13 1.1 17.92 ± 0.03 0.8	17.65 ± 0.03	0.1
	LARGE	LAser Ranging to GNSS Experiment
	LOD	Length of Day
	LRA	Laser Retro-reflector Array
	MC	Minimum Constraints
	NAVSTAR NAVigation Satellite Timing and Ranging Global Positioning System
	NP	Normal Point
	PCV	Phase Center Variations
	PCO	Phase Center Offsets
	QIF	Quasi-Ionosphere-Free
	QUASAR QUASi stellAR objects
	RB	Range Bias
	RMS	Root Mean Square
	RP	Reference Points
	SAO	Satellite Antenalla Offset
	SLR	Satellite Laser Ranging
	STARS	Sequential T-test Analysis of Regime Shifts
	TGP	Tide Generating Potential
	TP	True Positive
	TRF	Terrestrial Reference Frame
	TRS	Terrestrial Reference System
	UN	United Nations
	UT	Universal Time
	UTC	Coordinate Universal Time
	VLBI	Very Long Baseline Interferometry
	VMF1	Vienna Mapping Function
	WRMS	Weighted Root Mean Square

Full text available at: http : //ggim.un.org/docs/A 69 L53 E.pdf . Relevant information can be found also on the official website of the United Nations Global Geospatial Information Management (UN-GGIM) Working Group on the Global Geodetic Reference Frame (GGRF), www.unggrf.org .

available at ftp://igs-rf.ensg.eu/pub/

The complete list of satellites currently tracked by the ILRS network is available at http: //ilrs.gsfc.nasa.gov/missions/satellite_missions/current_missions/

The Data Handling File is available at http://ilrs.dgfi.tum.de/fileadmin/data_ handling/ILRS_Data_Handling_File.snx

Time-dependent CoM corections recommended by the ILRS for the analysis of LAGEOS data are available at http://ilrs.dgfi.tum.de/fileadmin/data_handling/com_lageos.txt

Core sites must meet the requirements detailed during the VI General Assembly of the ILRS, held in Nice March 28th, 2001. The meeting minutes are available through the ILRS website.

The lists of discontinuities employed in the computation of ITRF2014 are available at http:// itrf.ensg.ign.fr/ITRF_solutions/2014/doc/ITRF2014-soln-gnss.snx, for GNSS sites, and at http://itrf.ensg.ign.fr/ITRF_solutions/2014/doc/ITRF2014-soln-slr.snx, for SLR ones.

The precision of SLR frame parameters for the different cases analyzed in this study are shown in App. B. Relevant information can be retrieved combining the values presented in Tab. B.1, referred to a LAGEOS-only solution, with the gains reported in Fig. B.7 and B.8.

Acknowledgments

This work presents the results of three years of studies and efforts. I am deeply grateful to all the people who contributed to this achievement. I want to warmly thank my supervisors, Susanna Zerbini and Zuheir Altamimi for their guidance throughout the project, for their support and for the confidence they placed in me.

I am unspeakably thankful to Paul Rebischung, for the fruitful discussions, the brilliant ideas and his endless patience (not to mention the quiches!).

Thanks to my colleagues in Bologna, Maddalena and Efisio, and to those who welcomed me in Paris, Myriam, Ana, Leslie, Julien, Samuel, it was a real pleasure to work side-by-side and enjoy more than mere office hours.

Thanks to Cinzia and Pippo, for their help when I was absolutely lost with SLR data, and for the pleasant chats at conferences and meetings.

SLRtoLAGEOS:

LAGEOS initial state vector 6 per satellite LAGEOS empirical accelerations 5 per satellite Table 3.8 -Parametrization adopted for the quantities set up in the simulations. Concerning SLR observations to GNSS satellites one of the most significant aspects to consider is the fact that the number of acquired data varies tremendously from station to station, as can be seen in Fig. A.5. On a daily basis, for most of the stations, the number of recorded observations is actually insufficient to derive a non singular normal equation system including station coordinates, range biases, EOPs and the complete orbit parametrization selected for GNSS satellites.

Appendix B

Simulations -Supplementary material

In this appendix, results concerning the impact of SLR observations to the definition of GNSS orbit parameters are presented. Each plot illustrates the relative precision gain with reference to the corresponding GNSS-only solution. The investigated parameters include the six components of the initial state vector, i.e. satellite position and velocity at the beginning of the orbital arc, and the five empirical parameters absorbing unmodeled effects, see Sec. 2.2.2. The illustrated results complement the findings presented in Sec. 3.5.2 concerning physical frame parameters.

Within the simulations, the formal errors under study are computed satellitewise. Outcomes are here presented in aggregated form with box-plots extending from the lower to the upper quartile values of the data. A horizontal red line indicates the median, whiskers reach the minimum and the maximum points.

The presented plots show that the precision improvement obtained for orbital elements is considerably higher than the one observed for frame parameters. For the best SLR tracking scenario investigated with these simulations, median gains superior to 60% are systematically observed. Assuming a ranging precision at the cm level, instead, improvements in the order of 10% could be achieved if a network of about 50 stations perfectly distributed would be available. Comparable results hold for all the estimated parameters. It should be pointed out, however, that investments and technical improvements required to achieve these tracking capabilities are extremely demanding considering the current operational conditions.

Appendix C Identifying discontinuities in GNSS time series

When a piece-wise linear parametrization is selected for station coordinates in the realization of a terrestrial reference frame, it is fundamental to reliably locate the discontinuities affecting the relevant time series. In [START_REF] Bruni | Detecting discontinuities in GNSS coordinate time series with STARS: case study, the Bologna and Medicina GPS sites[END_REF], a screening strategy has been designed in order to apply the Sequential t test Analysis of Regime Shifts (STARS) [START_REF] Rodionov | A sequential algorithm for testing climate regime shifts[END_REF] to geodetic data-sets. Method and findings are reported in the following. The final publication is available at Springer via http://dx.doi.org/10.1007/s00190-014-0754-4

Global navigation satellite systems (GNSS) data are a fundamental source of information for achieving a better understanding of geophysical and climate-related phenomena. However, discontinuities in the coordinate time series might be a severe limiting factor for the reliable estimate of long-term trends. A methodological approach has been adapted from [START_REF] Rodionov | A sequential algorithm for testing climate regime shifts[END_REF]; [START_REF] Rodionov | Application of a sequential regime shift detection method to the Bering Sea ecosystem[END_REF]; [START_REF] Rodionov | Use of prewhitening in climate regime shift detection[END_REF] to identify both the epoch of occurrence and the magnitude of jumps corrupting GNSS data sets without any apriori information on these quantities. The procedure is based on the Sequential t test Analysis of Regime Shifts (STARS) [START_REF] Rodionov | A sequential algorithm for testing climate regime shifts[END_REF]. The method has been tested against a synthetic data set characterized by typical features exhibited by real GNSS time series, such as linear trend, seasonal cycle, jumps, missing epochs and a combination of white and flicker noise. The results show that the offsets identified by the algorithm are split into 48% of true-positive, 28% of false-positive and 24% of false negative events. The procedure has then been applied to GPS coordinate time series of stations located in the southeastern Po Plain, in Italy. The series span more than 15 years and are affected by offsets of different nature. The methodology proves to be effective, as confirmed by the comparison between the corrected GPS time series and those obtained by other observation techniques. I am honestly grateful to Dr Richard Biancale and Prof Carla Braitenberg for accepting to review the manuscript. I am looking forward to your comments and suggestions.

List of abbreviations

Voglio ringraziare la mia famiglia per il supporto continuo, l'affetto e il buonumore. Grazie di cuore.

Thanks to my friends. Thanks to the new ones who met me in this strange pe-riod of a student's life; and thanks to the historical ones...who still love me despite this last strange period of a student's life! Your support and the joy you bring into my life are the sparkle I look for every day.

Thanks to Davide, who was sure he would have gone through a third thesis unmentioned, but apparently was wrong! Thanks for taking good care of me and for multiplying the number of hours that my days have.