Analysis of Data From the GPS Reference Station at AAU Using GAMIT

GPS Technology - 2007

Group - 07gr1049

Isaac Nii Noi Tetteyo

Spring 2007 A ALBORG U NIVERSITY

Faculty of Engineering and Science

Aalborg University

Institute of electronic systems TITLE: Analysis of data from the GPS reference station at AAU Using GAMIT PROJECT PERIOD: 1. September 2006 - 16. August 2007 GROUP: 07gr1049 GROUP MEMBER: Isaac Nii Noi Tetteyo SUPERVISOR: Kai Borre

NUMBER OF DUPLICATES: 3 NUMBER OF PAGES IN APPENDIX: 5 TOTAL NUMBER OF PAGES: 70

Abstract This project concerns determination of the velocity of the GPS reference station at Aalborg University as a result of some geophysical phenomenon taking place. This is done by analyzing GPS data Obtained at the site over certain period of time using GAMIT/GLOBK software package. In the project, we describe the data archiving, processing procedures, and provide main results in terms of position time-series and velocities. We conclude that the reference station is in motion but to obtain accurately the magnitude and direction of motion demands certain time span of data to be processed.

Preface

This document reports the work of group 07gr1049 in the 9th and 10th semesters. It comprises of nine chapters. The rst gives an introduction and project objectives, the second and the third overview GPS theory and basic positioning methods and techniques. GAMIT/GLOBK overview is described in the fourth chapter followed by problem statement in the fth chapter. Data processing and analysis with some gures and results are presented in the sixth chapter. In the seventh chapter, we deal with land uplift in Fennoscandia with some baseline investigations. The nal results and discussion then comes in the eighth chapter. Conclusion and recommendations are given in the last chapter. Attached to the report is CD containing all processing results obtained in the project. The reader is welcome to use it. Some useful GAMIT input and output le formats, tables, etc. are described in the appendix. References are given as abbreviations of literature titles in square brackets that correspond to entries in the bibliography section at the end of the document. List of gures and list of tables are also provided. I want to thank my supervisors Professor Kai Borre and Laust Olsen with sincere heart for guiding me through the project, as well as assistant professor Henrik Have Lindberg for providing valuable information and advice on difculties with GAMIT, which arose at various stages of my project work. I am really grateful for all your support and assistance. Aalborg, August 13, 2007

Isaac Nii Noi Tetteyo

.....................
(signature)

Contents
1 Introduction 1.1 The Danish GPS Centre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Project Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 GPS Overview and Theory 2.1 General GPS System Description . . . . . . . . . 2.2 Geodetic Reference Systems . . . . . . . . . . . 2.2.1 World Geodetic System 1984 . . . . . . 2.2.2 International Terrestrial Reference Frame 2.3 GPS Observables . . . . . . . . . . . . . . . . . 2.4 Problems Associated with GPS Measurements . . 2.4.1 Ephemeris Errors . . . . . . . . . . . . . 2.4.2 Satellite Clock Errors . . . . . . . . . . . 2.4.3 Receiver Clock Errors . . . . . . . . . . 2.4.4 Ionospheric Effects . . . . . . . . . . . . 2.4.5 Tropospheric Effects . . . . . . . . . . . 2.4.6 Multipath . . . . . . . . . . . . . . . . . 2.4.7 Overall Error Budget . . . . . . . . . . . 2.4.8 Dilution of Precision . . . . . . . . . . . 2.5 GPS Applications . . . . . . . . . . . . . . . . . 3 Positioning Techniques 3.1 Point Positioning . . . . . . . . . . . 3.2 Precise Point Positioning . . . . . . . 3.3 Relative Positioning . . . . . . . . . . 3.3.1 Single Difference . . . . . . . 3.3.2 Double Difference . . . . . . 3.4 Comparison Between DGPS and PPP 5 5 6 8 8 8 9 9 10 10 11 11 11 11 12 12 12 13 13 15 15 15 16 16 17 17 18 18 19 19 21 22

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

4 Overview of GAMIT/GLOBK 4.1 GAMIT Processing Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 Parameter Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Overview of GLOBK Processing . . . . . . . . . . . . . . . . . . . . . . . . . 5 Problem Statement 5.1 Scope and limitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

CONTENTS 6 Data Processing and Analysis 6.1 Data Acquisition . . . . . . . . . . . . . . . . . . . . 6.2 Data Preparation for GAMIT Processing . . . . . . . . 6.2.1 RINEX observation and navigation les . . . . 6.2.2 Preparing the L-le . . . . . . . . . . . . . . . 6.2.3 Creating the station information le . . . . . . 6.2.4 Creating a scenario le . . . . . . . . . . . . . 6.2.5 Control les for the analysis (sittbl. and sestbl.) 6.3 Generating G- and T-les from external ephemerides . 6.3.1 Global Files . . . . . . . . . . . . . . . . . . . 6.3.2 Preparing GAMIT to Run . . . . . . . . . . . 6.3.3 Running GAMIT . . . . . . . . . . . . . . . . 6.3.4 Output les after running GAMIT . . . . . . . 6.3.5 Evaluating the solutions . . . . . . . . . . . . 6.4 Automatic Batch Processing . . . . . . . . . . . . . . 6.4.1 Processing Data from AAUC Site . . . . . . . 6.5 Data Preparation for GLOBK . . . . . . . . . . . . . . 6.5.1 Running Glred . . . . . . . . . . . . . . . . . 6.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 23 23 24 25 25 25 26 29 30 31 35 36 36 37 39 40 41 43 52 53 58 61 61 66 67 67 69 69 69

Land Uplift in Fennoscandia 7.1 Baseline Investigation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Final Results and Discussion Conclusion and Recommendations 9.1 Recommendations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8 9

A GAMIT File Formats A.1 Summary of GAMIT Processing . . . . . . . . . . . . . . . . . . . . . . . . . A.2 GAMIT Input and Output Files . . . . . . . . . . . . . . . . . . . . . . . . . . B GLOBK Processing Summary, File Formats, and Examples B.1 Summary of globk analysis of GPS data . . . . . . . . . . . . . . . . . . . . B.2 Example of globk/glorg output . . . . . . . . . . . . . . . . . . . . . . .

Chapter

1
Introduction

The technological advancement in Global Positioning System (GPS) with regard to data handling and processing has provided the ability and versatility in provision of precise GPS observations in a fraction of time. Many GPS applications such as modeling post-seismic velocity and displacement demand high accuracy positioning for which extensive data acquisition and processing is required. This is due to the fact that observations obtained with a GPS receiver is associated with some biases and errors and their impact affect the accuracy of the nal results. There is also some kind of geophysical effects which cause gradual displacements of point positions from their actual positions over a period of time. Most of these errors and occurrence can be eliminated or mitigated through the process of Differential GPS (DGPS) or Precise Point Positioning (PPP). DGPS is only applicable when the stations are closely spaced (e.g. 35 Km). But for GPS stations which are located further apart (e.g. 500 Km) they do not observe the same satellites and as such this technique is no longer applicable. These techniques will be discussed in detailed later in the report.

1.1 The Danish GPS Centre

The Danish GPS Center (DGC) was established 1996 as a research centre at Aalborg University. The basic motivation for DGC is to follow, support, and inuence the technological progress of satellite based navigation. DGC has its roots in surveying, geodesy, mathematical modelling, and digital signal processing. DGC continuously operates a reference station at Aalborg University. This station provides GPS carrier phase and code range measurements. Surveyors, engineers, scientists, and students can apply these data to position points at which user data have been collected. The objectives of the DGC are: To follow recent trends and latest evolutions in international satellite based positioning systems at the system level. To acquire knowledge about novel receiver technology; the antennas, the analog radio frequency parts, and the digital hardware for data processing and interfacing. To build up an understanding of the appropriate digital algorithms and the computational aspects related to the execution of these and their relation to specied accuracy. to develop and establish a useful collection of optimized software for post processing of GPS data with reference to still improving position accuracy, to build up networks for communication and discussion of GPS related information. 5

CHAPTER 1. INTRODUCTION To professionally assist small and medium sized industrial companies reaching for technical solutions involving GPS technology. To establish and maintain international contacts at universities and companies To participate in international research and development projects To maintain acquaintance with and encourage the development of new commercially available GPS products and to carry out systematic tests for evaluation of these. To generally promote the interest and understanding of the global navigation satellite systems (GNSS). [1] . From the above objectives of DGC, it is very clear that the accuracy of the position obtained from the GPS receiver is very crucial for the various users.

1.2

Project Objectives

There are some geodynamic phenomena such as earth crust movements and other deformations of the earth which changes position of points over time. Plate tectonic motion causes stations on earths crust to move horizontally. such movements can be 5-10 cm per year or even larger. [? ]. Therefore, for precise positioning, a time tag has to be associated with the coordinates of stations. An alternative is to estimate velocity in addition to the position coordinates In order to estimate accurately the positions and velocities of these points, GPS observations over several months has to be processed. This kind of investigation also requires a high precision geodetic measurements which cannot be achieved with a stand-alone GPS receiver. This requires that the reference station is tied to a network of stations with distances of hundreds of kilometers or more away from the vicinity in order for the extent of movements and velocities to be determined. The DGC reference station is part of a network with inter station distances of 500 km or more. There are various geophysical GPS research analysis software packages such as GAMIT/GLOBK, GYPSY, and Bernese used to investigate this kind phenomenon but in this case it is only GAMIT/GLOBK which will be used. There are several reasons why GAMIT/GLOBK is used, rstly it is the only software among the three which is available on the university servers/workstations. Secondly, it thoroughly estimate almost all the parameters affecting positioning by GPS observations. These parameters are estimated by precise a priori information about; satellite and their orbits, antenna phase centres, polar motion of the Earth, tectonic plate movements, position of the Sun and Moon. GAMIT/GLOBK will be described in details later in the report. The main aim of this project is to analyze the data from the GPS reference station at DGC using GAMIT/GLOBK software. From the analysis we expect the following: 1. Detection of earth crust movements. 6

1.2. PROJECT OBJECTIVES 2. Estimate the rate of movement or velocity. 3. Determine the direction of movements.

Chapter

GPS Overview and Theory

In order to understand and appreciate the GPS and how it works, it is very important to gain basic knowledge of the principles behind GPS and that is the main objective of this chapter. It specically deals with brief overview of GPS System, Geodetic Reference Systems, GPS Observables, GPS Measurement Errors, and GPS Applications.

2.1

General GPS System Description

The Global Positioning System (GPS) is a satellite-based navigation system for accurate and instantaneous position determination and timing. The primary objective of GPS was to meet the needs of the US military and national security, in regards to positioning and timing, on a 24hour per day basis all around the world and under all weather conditions. The satellites transmit at frequencies L1 at 1575.42 MHz and L2 at 1227.60 MHz modulated with two types of codes and the navigation message. The codes are Precise Code (P-Code) and the Coarse Acquisition Code (C/A Code). The P-Code generated at the fundamental frequency (i.e. 10.23 MHz) is available on both both L1 and L2. The C/A Code modulated on only the L1 frequency.The navigation message is generated at a low frequency of 50 Hz and is modulated on both L1 and L2 carrier frequencies. It contains information on the ephemerides of the satellites, GPS time, clock behaviour, and system status messages. The satellite constellation consists of 29 to 30 active satellites and three spare evenly spaced in six orbital planes. There are control/monitoring station networks around the world for monitoring the status and health of the GPS satellites.

2.2

Geodetic Reference Systems

In order to estimate accurately the position of a receiver on the earths surface based on satellite positions and ranges, it is very necessary that a common geodetic coordinate system is dened for all satellites and the receiver. A geodetic coordinate system is dened as a set of rules for specifying how coordinates values to be assigned to positions on the surface of the earth dened in the X, Y , and Z axes [Fre96]. This requires that the point must be co-rotating with the Earth in its diurnal motion in space. This is referred to as Terrestrial Reference System (TRS). The earth rotational axis will not hold still, it wanders slightly with respect to the solid earth in a very slow oscillation called polar motion. The actual displacement caused by the wandering does not exceed 12 metres [Sic01]. Nevertheless, TRS would be useless if the earths rotational axis wobbles. This can be computed through Very Long Baseline Interferometry (VLBI), Satellite Laser Ranging (SLR) and GPS observations. From these measurements a mean position 8

2.2. GEODETIC REFERENCE SYSTEMS and origin of the earths rotational pole has been estimated. This is called Conventional Terrestrial Pole (CTP). Accurate positioning within TRS is also affected by number of phenomena such as plate tectonic movements, solid earth tides, and ocean loading displacements. These corrections can be computed through complex mathematical expressions which are included in the GAMIT/GLOBK software package. The Conventional Terrestrial Reference System (CTRS) is dened with its origin at the centre of mass of the earth, the geocentre. The Z-axis passes through the CTP, the X-axis is line from the geocentre through the intersection of the zero meridian with the equator, and the Y -axis is extended from the geocentre along a line perpendicular from the X-axis in the same mean equatorial plane. They both rotate with the earth as part of the right-handed orthogonal system. The three dimensional cartesian coordinates (x, y, z) derived from this system are known as Earth-Centred-Earth-Fixed (ECEF) Coordinates. Theoretically, CTRS is an ideal TRS, but in practice, this is not the case since the geocentre of the earth is an imaginary quantity. One way of estimating CTRS is by using globally distributed points on the earths surface. The more accurate the coordinates, the denser the concentration of the points, and therefore the more easier and accurate is the realization of CTRS. Its realization is known as terrestrial reference frames. There are two terrestrial reference frames used in relation to GPS. These are The World Geodetic System 1984 (WGS84) and the International Terrestrial Reference Frame (ITRF) [Lei04, Sic01, Eng01].

2.2.1 World Geodetic System 1984

The WGS84 has evolved signicantly since its creation in mid-1980s. The original WGS84 reference frame established in 1987 was realized through a set of Navy Navigation Satellite System or TRANSIT (Doppler) Station coordinates. These set of estimated coordinates had uncertainty of 1-2 m [pbtNWUC06]. However, there has been signicant improvement in realization of the current WGS84 reference frame through the use of advanced GPS techniques. The accuracy derived in 1996 has been shown to be in the order of 5 cm [pbtNWUC06]. This type of terrestrial reference frame is the most commonly used by GPS since the GPS satellites broadcast their positions in WGS84.

2.2.2 International Terrestrial Reference Frame

The best geocentric reference frame currently available is the ITRF. The ITRF is maintained by the International Earth Rotation Service (IERS) which monitors the Earth Orientation Parameters (EOP), for the scientic community through a global network of observing stations. This is done through Space Geodesy techniques such as; Very Long Baseline Interferometry (VLBI), Lunar Laser Ranging (LLR), Satellite Laser Ranging (SLR), Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS) and GPS observations. Station velocities are also estimated to account for some geophysical effects such as plate tectonic movements. The positions of the observing stations are now considered to be accurate to the centimeter level 9

CHAPTER 2. GPS OVERVIEW AND THEORY [Sic01].

2.3

GPS Observables

The two most important observables used in GPS positioning are the pseudoranges and the carrier phases. The pseudorange is related to the distance between the satellite and the receivers antenna at some epochs of emission and reception of the codes. The pseudorange is the time offset multiplied by the speed of light and it is biased by the lack of time synchronization between the clock in the GPS satellite and the clock in the GPS receiver. The foundation of a pseudorange measurement is the correlation of the codes received from GPS satellite with replicas of those codes generated within the receiver. The codes generated by the receiver are based on the receivers own clock, and the codes of the satellite transmissions are generated by the satellite clock. The equation for the pseudorange observable between receiver i and satellite k is modeled as follows Pik = k + Iik + Tik + ctk + k (2.1) i i i where Pik , k , Iik , Tik , ctk are pseudorange, true range, ionospheric delay, tropospheric delay, i i and clock bias respectively between satellite k and receiver i scaled to range units. k represent i errors which come from pseudorange measurements and multipath. The carrier phase is the phase of the received carrier with respect to the phase generated by an oscillator in the GPS receiver. The measured phase consists of fractional component i (t0 ) and integral part at an instant of reception t0 . The integral component is a number of whole wave cycles, which is measured relatively to some initial lock-on value. The receiver measures the fractional phase, and keeps the track of the changes to the phase. It is more precise than the pseudorange but the problem is that the receiver cannot distinguish one cycle of a carrier from another. The initial phase is undetermined, or ambiguous, by an integer number of cycles N . This unknown quantity is the cycle ambiguity. Carrier phase can be converted into equivalent distance by similar expression as follows k (t) = k Iik + Tik + ctk + Nik + (i (t0 ) k (t0 )) + k . i i i i (2.2)

Here, cycles phase k is not delayed but advanced by Iik . Furthermore, is the wavelength, i k (t0 ) is the unknown fractional part of the phase near the satellite at reception, and k is noise. i Equation (2.2) is very similar to the pseudorange equation (2.1), the major difference being the presence of ambiguity term N . This uncertainty can be resolved by the method of differential positioning.

2.4

Problems Associated with GPS Measurements

All GPS measurements, be they pseudorange or carrier phase, are affected by biases and errors which affect the accuracy of the position determination. There are several sources of bias with 10

2.4. PROBLEMS ASSOCIATED WITH GPS MEASUREMENTS varying characteristics of magnitude, periodicity, satellite-receiver dependency, etc. Among these biases are ephemeris errors, satellite clock errors, receiver clock errors, ionospheric effects, tropospheric effects and satellite-receiver geometry. Most of the information in this section was taken from [Lei04, Sic01, Fre96, Riz99, SB97].

2.4.1 Ephemeris Errors

The ephemeris information used to calculate the satellite position is based on the observations from the monitor stations of the space segment. The data is processed at the Master Control Station and the satellite navigation message information is uploaded to every satellite. The ephemeris error is therefore the discrepancy between the true position (and velocity) of the satellite and its broadcast ephemeris. This grows from the time of upload by a monitor station until the next upload. Parkinson & Spilker Junior (1996) estimate the root-mean-square (rms) error as 2.1 m but it might be smaller in recent years to as low as 1 m.

2.4.2 Satellite Clock Errors

Although GPS satellites use high quality cesium or rubidium atomic clocks for time-keeping and signal synchronization, there are unavoidable clock errors which change with time. Rubidium oscillator is correct to about 1012 and that of cesium is correct to about 1013 . The offset could reach 107 seconds in a day, multiplied by the velocity of light gives 26 m. As all observation made at an instant to a particular satellite, all GPS receivers are affected by equal magnitude of the same satellite clock error, the principle of differential positioning can be used to eliminate this error.

2.4.3 Receiver Clock Errors

GPS receivers are equipped with inexpensive quartz crystal oscillators with low accuracy. The offset between receiver clock time and GPS Time is the receiver clock error which affect all satellite receiver ranges measured at a particular epoch. Double differencing can be used to eliminate this error.

2.4.4 Ionospheric Effects

The ionosphere extends between approximately 50 and 1500 km above the earth and it is characterized by the presence of free electrons and positively charged atoms and molecules called ions. This is as a result of the gas molecules being excited by solar radiation. The total electron content (TEC) equals the number of free electrons in the column of unit area along which the signal travels between the satellite and the receiver. TEC varies as a function of latitude of the receiver, the season, the time of the day the observation of the satellite signal is being made and the level of solar activity at the time of observation. As the electromagnetic GPS signal 11

CHAPTER 2. GPS OVERVIEW AND THEORY Table 2.1: Overall error budget Source Impact, [m] Ephemeris data 2 Satellite clock 2 Ionosphere 4 Troposphere 0.5 - 1 Multipath 0-2 propagate through the medium, dispersion occurs and the free electrons delay the pseudorange and advance the carrier phase by equal magnitude. The amount is directly proportional to the TEC and inversely proportional to the carrier frequency. GPS frequency delays or advances can be up to 50 m for signals at the zenith to as much as 150 m for observations made at the receivers horizon [Klo91]. An effective procedure to deal with this error is to take advantage of the frequency dependence of the ionospheric effect by using a dual-frequency receiver. Measurements are made on both L1 and L2 frequency signals and combining them in a linear form, the delay is eliminated since the impact on L1 and L2 is different.

2.4.5

Tropospheric Effects

The troposphere extends from the surface of the earth to about 50 km above the earth and it is neither ionized nor dispersive. GPS signals travelling through this medium will experience delay that is a function of elevation and altitude of the receiver. Tropospheric effect is dependent on the atmospheric pressure, temperature, and humidity. The bias ranges from approximately 2 m for satellites at the zenith to about 20 m for satellites at an elevation angle of 10o [Bru93]. The propagation of GPS signal in this medium is frequency independent, therefore this effect cannot be removed by combining observations made on L1 and L2 frequencies. There are many models available for this correction.

2.4.6

Multipath

Multipath effects are propagation errors arising from interference of the direct signal by reected signals from water, metallic surfaces, and nearby buildings. The combined direct and reected signals will give rise to incorrect pseudorange observation. Errors which arise as a result of multipath cannot be reduced by the technique of DGPS, since they depend on the local reection geometry near each receiver antenna. The remedies for multipath lies in site selection and effective antenna design to lter out multipath effects using advanced signal processing.

2.4.7

Overall Error Budget

Table 2.1 presents standard errors of a single frequency receiver as given in [SB97] expressed in position uncertainty in meters. As it can be seen the most signicant errors come from 12

2.5. GPS APPLICATIONS ionosphere and then multipath, satellite clock, and others.

2.4.8 Dilution of Precision

The distribution of satellites above the observers horizon has a direct effect on the accuracy of the position determination. The accuracy of GPS position is subject to a phenomenon called dilution of precision (DOP). It is the measure of geometry of the visible satellites with respect to one another and the GPS receivers. A low DOP factor is good, a high DOP factor is bad. This is sufce to say that, when satellites are in optimal conguration for a reliable GPS position, the DOP value is low; when they are not, the DOP value is high. If all the four satellites required for 3-dimensional positioning are all crowded together in one part of the sky, then there is likely to be a less accurate position and the DOP will be high. DOP is like a warning that the actual errors in a GPS position are liable to be larger than you might expect. It must be claried that it not the errors themselves that are directly increased by the DOP factor, it is the uncertainty of the GPS position that is increased by the DOP factor. There are a number of DOP components; there is horizontal dilution of precision (HDOP) and vertical dilution of precision (VDOP) when the uncertainty of a solution for positioning has been isolated into its horizontal and vertical components, respectively. When both horizontal and vertical components are combined, the uncertainty of 3-dimensional position is called position dilution of precision (PDOP). In time domain, we have time dilution of precision (TDOP), which indicates the uncertainty of the clock. The combination of all the above components form geometric dilution of precision (GDOP). There is also relative dilution of precision (RDOP), that includes the number of receivers, the number of satellites they can handle, the length of the observing session as well as the geometry of the satellites conguration. Visible satellites which are widely and evenly spaced in the sky indicates good geometry and therefore low.

2.5 GPS Applications

GPS has become an important tool for various applications, the general overall uses are numerous and this report cannot address all the GPS applications. Therefore we will narrow the classication into navigation and surveying applications. Navigation was one of the planned primary uses of GPS. Both military and civilian uses of the system in this mode are similar in that users wish to know their spatial locations as precisely as possible. In aircraft navigation in which accuracy requirements are fairly high during approach and landing, xed receivers must broadcast range corrections (and carrier phases) to incoming aircraft so that they can compute more accurate positions as they approach the runway. GPS-aided approach and landing will be an economical answer to airport control. GPS can also be used for attitude determination of the aircraft. In surveying, GPS provides a powerful geodetic tool for monitoring global changes over a period of time which is the key for understanding long-term geodynamic phenomena. This applications includes measuring crustal deformations, volcanic uplift, plate tectonics, and earth 13

CHAPTER 2. GPS OVERVIEW AND THEORY rotation. In the past VLBI and SLR techniques have been used for this purpose. GPS is not presently capable of replacing these techniques but will be used to augment them and provide more cost-effective solutions for geodetic problems. The objective of this project falls in this category of the GPS application [HLC01].

14

Chapter

Positioning Techniques
In the previous chapters, errors associated with GPS positioning were described. In the this chapter, we described how this error are handled and mitigated using various positioning techniques. GPS positioning tecniques be it kinematic or static can be subdivided into Point Positioning, Precise Point Positioning, and Relative Positioning. These can either be real-time or postprocessed depending on the application. This chapter will focus on the basic positioning techniques in only static mode.

3.1 Point Positioning

Point positioning refers to the estimation receiver antenna coordinates xi and the receiver clock error ti using pseudorange observables. The role of the carrier phase is limited to smoothing the pseudoranges, if used at all. There are several simplifying assumptions made in point positioning, the satellite position xk at transmission time are assumed known and available from the broadcast ephemeris. While the receiver clock error is estimated at every epoch, neglecting the residual satellite clock error tk . But the satellite clock broadcast corrections must be applied. Ionospheric and tropospheric delays are are also computed from models. Hardware delays and multipath are neglected. Using the simplifying assumptions made above, we can write equations of the type i k = cti (3.1) i k is the geometric range given by i k = xk xi where

xk xi =

(xk xi )2 + (y k yi )2 + (z k zi )2

(3.2)

The four unknowns xi , yi , zi , and ti can be computed using four pseudoranges measured simultaneously. The effect of the earth rotation during signal travel time must be incorporated in equation (3.1). The basic requirement, however is that there are four satellites visible at a given epoch. Point positioning depends on the accuracy of the navigation message and the constellation of the satellites used. In practice, not just four satellites are observed but all satellites in view in order to achieve redundancy and better geometry.

3.2 Precise Point Positioning

The vast majority of GPS application utilizes the principle of relative positioning, however in the late 1990s, the Jet Propulsion Laboratory (NASA) powered a new technique that did not re15

CHAPTER 3. POSITIONING TECHNIQUES quire differencing to obtain precise position. This is labeled as Precise Point Positioning (PPP). PPP refers to a centimeter accuracy of a single static receiver using a long observation series and to subdecimeter accuracy of a roving receiver using free ionospheric and carrier phase functions [Lei04]. The receiver clock error and and the zenith tropospheric delay are estimated for each epoch. When using PPP, simplifying assumptions are avoided; all known corrections must be applied to the observations and the corrections must be consistent. The satellite position at transmissions are computed from the postprocessed precise ephemeris available from International Global Navigation Satellite System Service (IGS) or its associated processing centres. A crucial element in achieving centimeter position accuracy with PPP is accurate satellite clock corrections which are part of the precise ephemeris. The ionospheric effect is eliminated by using a dual frequency GPS receiver.

3.3

Relative Positioning

One of the popular ways of tackling the errors in GPS positioning is by the method of Relative Positioning also known as Differential GPS (DGPS). This is done by computing differences between simultaneous observations from two receivers. The vector between the two receivers is determined which is often called baseline vector or simply baseline. Basically, three main types of differences which are frequently used are Single Difference, Double Difference and Triple Difference. But only the rst two will be discussed since triple difference is not often used due to the fact that it looses its geometric strength over time and it is less stable. But triple difference is very important in estimating the integer ambiguities which become constant over time.

3.3.1

Single Difference

This is also known as between-receivers difference, and it refers to the difference in the simultaneous code and carrier phase measurements. If a signal transmitted at a given instant by a satellite k is observed by two receivers i and j closed to each other, then both receivers are observing the same satellite at the same time. They will therefore experience similar atmospheric effects, have identical orbital errors and also with the same satellite clock error. The satellite k antenna offset also cancels out. Single difference code Pij and phase k (t) observations is ij given as follows: k Pij = Pik Pjk = k + ctij + k (3.3) ij ij
k k (t) = k (t) k (t) = k + ctk + Nij + ii (t0 ) + k ij i j ij ij ij

(3.4)

From equations (3.3) and (3.4), it can be seen that factors which are not eliminated by single difference is the receiver clock errors tij , and the integer cycle ambiguities Nij in the carrier phase observable. 16

3.4. COMPARISON BETWEEN DGPS AND PPP

3.3.2 Double Difference

This is also known as between-satellite difference. It is the difference in measurement from two GPS satellites k and l as measured simultaneously at two receivers i and j. Subtracting two single difference yields the double difference and eliminates the receiver clock errors. The kl double difference code Pij and phase kl (t) is given as follows: ij
l k kl Pij = Pij Pij = kl + ij kl ij

(3.5) (3.6)

kl kl = k l = kl + Nij + kl ij ij ij ij ij

It can be seen from equations (3.5) and (3.6) that the only unknown quantity in the double kl difference observation to be determined is the difference ambiguity Nij . It plays an important role in accurate DGPS using double differences.

3.4 Comparison Between DGPS and PPP

The accuracies of DGPS degrade with increase in baseline but ofcourse, this is not the case with PPP since it does not rely on a base station. The main drawback with PPP is the time delay in the availability of the precise satellite orbits and clocks consequently, the most precise positions are not available until sometime. In DGPS, any error incorporated in the base station is propagated to the roving station but this may not be the case in PPP since all the points are processed independently of each other. From the above discussions, it can be seen that the strength of PPP lies in the weakness of DGPS and vice versa. This shows that their combination will be a very powerful technique in GPS positioning.

17

Chapter

Overview of GAMIT/GLOBK
GAMIT/GLOBK is a comprehensive GPS analysis package developed at MIT, the HarvardSmithsonian Center for Astrophysics (CfA), and the Scripps Institution of Oceanography (SIO) for estimating station coordinates and velocities, stochastic or functional representations of post-seismic deformation, atmospheric delays, satellite orbits, and Earth orientation parameters. The software is designed to run under any UNIX operating system supporting X-Windows. The maximum number of stations and atmospheric parameters allowed is determined by dimensions set at compile time and can be tailored to t the requirements and capabilities of the analysts computational environment. There are also C-shell Scripts with name beginning with sh which come with the package to control processing. The main aim of this chapter is to give a brief introduction to GAMIT/GLOBK but the rest will be discussed in detailed later in the report.

4.1

GAMIT Processing Algorithms

GAMIT incorporates difference-operator algorithms that map the carrier beat phases into singly and doubly differenced phases. These algorithms extract the maximum relative positioning information from the phase data regardless of the number of data outages, and take into account the correlations that are introduced in the differencing process. In the presence of cycle slips, initial processing of phase data is often performed using triple difference or doppler observations in order to obtain a preliminary estimates of station or orbital parameters. GAMIT software uses triple differences in editing but not parameter estimation. Rather it allows estimation of extra free bias parameters whenever automatic editor has agged an epoch as a possible cycle slip. GAMIT is composed of distinct programs which perform the functions of preparing the data for processing (makexp and makex), generating reference orbits for the satellites (arc), computing residual observations (o-cs) and partial derivatives from a geometrical model (model), detecting outliers or breaks in the data (autcln), and performing a least-squares analysis (solve). Although the modules can be run individually, they are tied together through the data ow, particularly le-naming conventions, in such a way that most processing is best done with shell scripts and a sequence of batch les set up a driver module (fixdrv) for modeling, editing, and estimation. Though the data editing is almost always performed automatically, the solution residuals can be displayed or plotted so that problematic data can be identied (cview). It must be stated emphatically that blind reliance on the solution without any thorough analysis and the implication of the results can lead to disaster, as always. 18

4.2. OVERVIEW OF GLOBK PROCESSING

4.1.1 Parameter Estimation

GAMIT incorporates a weighted least-squares algorithm to estimate the relative positions of a set of stations, orbital and Earth-rotation parameters, zenith delays, and phase ambiguities by tting to doubly differenced phase observations. Since the functional (mathematical) model relating the observations and parameters is non-linear, GAMIT produces two solutions, the rst to obtain coordinates within a few decimeters, and the second to obtain the nal estimates. The GAMIT solution is not usually used directly to obtain the nal estimates of station positions from a survey. Rather, GAMIT is used to produce estimates and an associated covariance matrix of station positions and (optionally) orbital and Earth-rotation parameters which are then input to GLOBK to estimate positions and velocities. In order not to bias the combination, GAMIT generates the solution used by GLOBK with only loose constraints on the parameters, dening the reference frame only at the GLOBK stage by imposing constraints on station coordinates. Since phase ambiguities must be resolved (if possible) in the phase processing, however, GAMIT generates several intermediate solutions with user-dened constraints before loosening the constraints for its nal solution. In parameter estimation based on least-squares, the conventional measure of goodness-oft is the 2 (chi-square) statistic, dened for uncorrelated data as the sum of the squares of each observation residual (post-t observed minus computed observation, "o-c") divided by its assigned uncertainty. The value is usually normalized by dividing by the degrees of freedom (df ), which is the number of observations minus the number of parameters estimated, so that the ideal value for properly weighted observations is 1.0. In a multi-parameter solution, correlations arise so the computation of 2 /df in GAMIT involves a complex matrix operation. Later in the report we will discuss how the 2 statistic is used to assess the quality of a GPS analysis. It is important to state that its value depends not only on the data noise and processing models, but also on how realistic are the a priori errors assigned to the phase observations and/or the quasi-observations used by GLOBK. [HKM06].

4.2 Overview of GLOBK Processing

GLOBK is a Kalman lter whose primary purpose is to combine solutions from the processing of primary data from space-geodetic or terrestrial observations. It accepts as data, the estimates and associated covariance matrices for station coordinates, earth-rotation parameters, orbital parameters, and source positions generated from analyses of the primary observations. These primary solutions are performed with loose a priori uncertainties assigned to the global parameters, so that constraints can be applied uniformly in the combined solution. There are three common modes, or applications, in which GLOBK is used: 1. Combination of individual sessions (e.g., days) of observations to obtain an estimate of station coordinates averaged over a multi-day experiment. For GPS analyses, orbital parameters can be treated as stochastic, allowing either short or long arc solutions. 19

CHAPTER 4. OVERVIEW OF GAMIT/GLOBK 2. Combination of experiment-averaged (from .1) estimates of station coordinates obtained from several years of observations to estimate station velocities. 3. Independent estimation of coordinates from individual sessions or experiments to generate time series assessment of measurement precision over days (session combination) or years (experiment combination). Some things GLOBK cannot do. 1. GLOBK assumes a linear model. Therefore any large adjustments to either station positions or orbital parameters (>10 m for stations and >100 m for satellite orbits) need to be iterated through the primary processing software to produce new quasi-observations. 2. GLOBK cannot correct deciencies in the primary (phase) analysis due to missed cycle slips, "bad" data, and atmospheric delay modeling errors. You cannot eliminate the effect of a particular satellite or station at the GLOBK stage of processing, though GLOBK can be useful in isolating a session which is not consistent with the ensemble and in some cases the effect of a station on the GLOBK solution can be reduced. 3. GLOBK cannot resolve phase ambiguities: the primary GPS solution must be strong enough on its own to accomplish this. The need to combine sessions for ambiguity resolution is the one reason one might want to perform a multi-session solution with primary observations. [HKM06] GLOBK operates through distinct programs, which can be invoked with a single command or run separately. The primary functions are to combine quasi-observations either from multiple networks and or epochs (glred or globk), and to impose on this solution a reference frame appropriate to the scientic objective (glorg). It must be emphasized that globk and glred are the same program, just called in different modes: glred to read data from one day at a time for generating time series, globk for stacking multiple epochs to obtain a mean position and /or velocity.

20

Chapter

5
Problem Statement

Geodetic measurements for high-precision GPS applications can be achieved using carrier beat phase of the GPS signal. In the previous chapters we discussed that, the dominant source of error in a phase measurement or series of measurements between a single satellite and ground station is the unpredictable behavior of the time and frequency standards (clocks) serving as reference for the transmitter and receiver. Tropospheric effects are largely removed by either applying a model which attempts to mathematically simulate the signal delay or by estimating the signal troposphere delay along with the receiver coordinates. Ionospheric effects are removed by observing both GPS frequencies (L1 and L2) and combining the two observations to derive an ionosphere-free observation. The 3-dimensional accuracy of estimated baseline, as a fraction of its length, is roughly equal to the fractional accuracy of the orbital ephemerides used in the analysis [HKM06]. The accuracy of the Broadcast Ephemerides computed regularly by the Department of Defense using pseudorange measurements from less than 10 tracking stations is typically 1 5 parts in 107 (2-10 m), well within the design specications for the GPS system but not accurate enough for the study of crustal deformation. By using phase measurements from a global network of over 100 stations, (IGS) is able to determine the satellites motion with an accuracy of 1 part in 109 (2 cm) [HKM06]. Therefore, errors in satellite positions can be reduced by using precise satellite orbits available from the IGS and any remaining error (except multipath) largely cancels over short distances. That leaves satellite and receiver clock errors as the dominant errors to be dealt with. This is also mitigated through Relative Positioning or Precise Point Positioning (PPP) depending on the distances between stations. The largest difference between relative processing and PPP is the way that the satellite and receiver clock errors are handled. PPP uses highly precise satellite clock estimates. These satellite clock estimates are derived from a solution using data from a globally distributed network of GPS receivers. Instead of between-satellite differencing to remove receiver clock error, PPP estimates these as part of the least- squares solution for the coordinates. The position of points estimated on the Earths surface are in constant motion as a results of earths crust movements and due to tectonic plate motion. The motions of these points are largely slow and smooth in nature, with the exception of regions where earthquake activity is high. In these regions, signicant surface displacements over a very short time period can and do occur. By using tracking data acquired from IGS stations as reference, the variations in site positions can be monitored through time and their velocities can be accounted. The purpose of the project is to analyze the GPS data at AAU Danish GPS Centre. This is achieved through the following: 21

CHAPTER 5. PROBLEM STATEMENT 1. Choosing 2 IGS Stations with good geometry to tie to DGC and with baselines as short as possible. In this case, IGS Stations chosen are ONSA (Onsala, Sweden) and POTS (Germany). 2. Downloading, extraction and preparation of required data using GAMIT version 10.21 3. Processing required data in 30-day batches, plotting time series and repeatability for detection of cycle slips and outliers using GLOBK. 4. Removing outliers and reprocessing to nally estimate station position and velocity using GLOBK 5. Plotting station position, velocity and its direction on a map. 6. Comparison of results obtained in each batch with each other and a conclusion drawn

5.1

Scope and limitation

Due to time constraints, limitations also exist for the project. Therefore in order to achieve the desired objectives, the following bounds are dened; 1. Only the rst 30 days of each year of data from DGC of the years 2000, 2001, 2004, 2005, and 2006 will be processed. Data is chosen in the same period of the year to minimized the unknown impact of seasonality on the GPS results, thus reducing aliasing problems in the incipient time series. 2. Only GAMIT/GLOBK is used in the data analysis and also some graphical representation will be displayed. 3. Ionospheric and Tropospheric effects, error due to multipath is not of interest to this project. 4. Displacements directly attributable to earthquakes may not be evident in these results

22

Chapter

Data Processing and Analysis

In order to achieve the desired objectives mentioned in the previous chapters, the necessary data must be acquired and in the correct format for processing. This chapter describes the process of data acquisition and preparation for processing.

6.1 Data Acquisition

First and foremost, we have the navigation and observation data les for each day of the year (doy) at the receiver for the AAU Reference Station stored on the DGC web site in compressed RINEX (The Receiver Independent EXchange) format. The required period of data is downloaded from the web site and uncompressed. The navigation data is discarded and precise ephemerides data is downloaded from Scripps Orbit and Permanent Array Center (SOPAC) data archives. Two IGS stations of POTS and ONSA to tie to ITRF are chosen. Observation data in compressed Hatanaka format for each doy of these IGS points are also downloaded from SOPAC. SOPAC converts all RINEX observation les to the Hatanaka ("d le") format. RINEX les more than 60 days old are available in Hatanaka format only, so all the observation les have to be converted to RINEX format since GAMIT can only read observation and navigation les in RINEX format. This is done by using c-shell script with the name crx2rnx which comes with GAMIT software package. In the subsequent sections, we will describe in detailed all input les required for GAMIT processing.

6.2 Data Preparation for GAMIT Processing

Before processing GPS data using GAMIT, it is necessary to understand some of its le naming conventions. This assures a unique denition for each experiment, facilitates data le management, and allows for ease of interactive processing and troubleshooting. First the information should be organized into sessions, dened as the spans during which all the stations tracks simultaneously the phase of two or more satellites. When a single session is processed in GAMIT, processing the remaining sessions is easy. We therefore describe how a single session is processed in this section. The rst step is to create a working directory and copy or link all the necessary les. The main les needed for a single session processing are: 1. RINEX observation and navigation les; 2. Station coordinates in the form of an L-le; 23

CHAPTER 6. DATA PROCESSING AND ANALYSIS 3. Receiver and antenna information for each site ( le station.info); 4. Satellite list and scenario (le session.info); 5. Control les for the analysis ( sittbl. and sestbl.); 6. G- and T-les from external ephemerides; 7. Links to the following global les: nutations (nutabl.); lunar and solar ephemerides (soltab. and luntab.); geodetic datums (gdetic.dat); leap seconds (leap.sec); spacecraft, receiver, and antenna characteristics (svnav.dat, antmod.dat, rcvant.dat) Earth rotation (pole., ut1.); ocean tides (stations.oct and grid.oct)

6.2.1

RINEX observation and navigation les

The acquisition of the RINEX observation and navigation les has been described earlier in the section. Data rate is 30 seconds, and all les contain observations for a 24 hour period, from 00 : 00 : 00 till 23 : 59 : 30 GPS time. File naming follows the convention "NNNNddds.yy" where NNNN is the site name ( e.g AAU C), ddd refers to the day of the year (e.g 024 means January 24th), s is the session number (e.g 0), and yy is the two last digit of the year. Examples of RINEX header is given below:
2 GPS Data Logger OBSERVATION DATA GRDL G (GPS) 01-Jan-2004 00:00 RINEX VERSION / TYPE PGM / RUN BY / DATE COMMENT COMMENT COMMENT MARKER NAME MARKER NUMBER OBSERVER / AGENCY REC # / TYPE / VERS ANT # / TYPE APPROX POSITION XYZ ANTENNA: DELTA H/E/N WAVELENGTH FACT L1/2 # / TYPES OF OBSERV INTERVAL TIME OF FIRST OBS TIME OF LAST OBS END OF HEADER

AAUC 1 BORRE 1 3427985.87 0 1 1 5 C1 5 2004 1 2004 1

GPS-CENTER ASHTECH Z-12 6J00 GEODETIC 603660.70 5326788.90 0 0 L1 1 1 L2 0 23 P1 0 59 P2 0 29

The RINEX observation data le contains the L1 and L2 carrier beat phases and pseudoranges, signal amplitudes, initial station coordinates and antenna offsets, start and stop times, and the identication of the satellites tracked in each receiver channel. 24

6.2. DATA PREPARATION FOR GAMIT PROCESSING

6.2.2 Preparing the L-le

The L-le contains the coordinates of all the stations to be used in the data processing. Two forms of coordinate formats are supported: geocentric (spherical) coordinates at the epoch of the observation, or Cartesian coordinates and velocities at a specied epoch. GAMIT software package contains the precise coordinates of all IGS stations in a folder called tables.The L-le used in this experiment is shown below:
Epoch 2006.3450: ONSA ONSA GPS POTS POTS GPS AAUC AAUC GPS From file itrf00.apr N57 56 19.54622 E011 41 12.45628 6428987.4781 N52 07 47.07811 E013 08 10.22143 6362329.8658 N57 03 00.69097 E009 46 33.79156 6345408.7828

Updated from laauc4.001 Updated from laauc4.001 Updated from laauc4.001

This is generated for a particular epoch from GLOBK apr le using the program gapr_to_l. In this case, geocentric coordinates is used, The station AAUC is not in the GLOBK apr le so its approximate cartesian coordinates from the RINEX le header is used to transform it to geocentric by using the script tf orm.

6.2.3 Creating the station information le

All of the receiver and antenna information specic to a particular site occupation are recorded in the station.info le. The values entered correspond to a single occupation, of either one day or a series of days. In this le format, the number of entries (columns) is variable and determined by a list beginning with the keyword SIT E. Although the number and order is arbitrary, the width (format) of each entry is rigid and hard-wired into the code. An example is shown below:
# New-style station.info written from old using conv_stnfo by 2003-11-06 09:28 * Ant Ht *SITE Station Name Session Start Session Stop ONSA ONSA Onsala 2004 1 0 0 0 9999 999 0 0 0 0.9950 POTS POTS Potsdam 2004 1 0 0 0 9999 999 0 0 0 0.0460 AAUC AAUC Aalborg 2004 1 0 0 0 9999 999 0 0 0 0.0050 rwk on

HtCod DHBCR DHPAB DHPAB

Ant N 0.0000 0.0000 0.0000

Ant E 0.0000 0.0000 0.0000

RcvCod TR8000 TR8000 ATUZ12

SwVer 2.80 3.20 7.10

AntCod ROGAOA TRBROG ASHL12

The most important among the entries in station.info are the antenna type (AntCod) and specication of how the height-of-instrument was measured (HtCod) since this directly affects the estimated heights from the analysis. This information in entered into the le in the form of keywords and later converted by GAMIT to L1 and L2 phase-center offsets. Entries for horizonal offsets (Ant N, Ant E) of the antenna from the monument are can also be entered. The station.info values are added to the coordinates of the monument in computing the antenna phase-center position.

6.2.4 Creating a scenario le

The scenario le le is also known as the session.info and it contains the start time, sampling interval, number of observations, and satellites (PRN s) to be used in generating X-les for each day. It does not correspond to the time-dependent scenarios used to program some receiver software, but rather includes all satellites to be used in the analysis. An example of session.info le is shown below: 25

CHAPTER 6. DATA PROCESSING AND ANALYSIS

# Session.info : free format, non-blank first column is comment #Year Day Sess# Interval #Epochs Start hr/min Satellites 2004 1 1 30 2880 0 0 1 2 3 4 5 6 7 8 9 10 11 13 14 15 16 17 18 20 21 22 23 24 25 26 27 28 29 30 31

The session.info le can be specic to a given experiment or contain all of the scenarios used for all the experiments processed. It can also be generated automatically to a specic experiment by the program makexp using the input start/stop time and the satellites available on the navigation le. To nd out what satellites are available on the raw data les, use the shell scripts rxscan and ficascan which scan through the RINEX or FICA les and generate the list of available satellites with their PRN s.

6.2.5

Control les for the analysis (sittbl. and sestbl.)

The site control table used in the experiment (sittbl.) has a form as shown below:
STATION POTS POTSDAM ONSA ONSALA AAUC AALBORG FIX NNN NNN NNN --COORD.CONSTR.-0.005 0.005 0.010 0.005 0.005 0.010 99.99 99.99 99.99 --EPOCH-001- * 001- * 001- * CUTOFF 15.0 15.0 15.0 APHS NONE NONE NONE CLK NNN NNN NNN KLOCK 3 3 3

The table contain the number of stations to be used in the experiment. The third column (FIX) in the le shown above indicates which, if any, of the four station coordinates are to be xed/free (Y/N) in the solution. Any a priori constraints are given under COORD.CONSTR. in units of meters for latitude, longitude, and radius. Large numbers have small weights and small numbers have large weights due to the nature of the covariance matrix (1/variance). For instance, since the local site of AAUC is unknown a value of 99.99/99.99/99.99 is used. Values commonly used on IGS site are 0.005/0.005/0.005. By specifying values for the number of observations to be used from one or more stations is specied as a range of epochs under the column (EPOCH). The elevation cutoff can be controlled under the column (CUTOFF). The column APHS is for invoking a model for variations in the phase centre of the receiving antenna. The next two entries control the way the receiver clock is handled. The CLK entry indicates whether or not an offset, rate or acceleration term is to be estimated when running GAMIT. This is most of the time set to "No" (N). KLOCK selects the way the clocks are modeled and can be very important to some receivers. For all receivers except MiniMacs, option 3 is chosen to show that receiver clock offset is to be estimated epoch-by-epoch using the pseudorange. A session control table (sestbl.) with only required and commonly used entries is shown below:
Session Table Processing Agency = MIT Station Constraint = Y Satellite Constraint = Y all ; Y/N ; Y/N e

(next two lines are free-format but all must be present) i n w M rad1 rad2 rad3 rad4...rad9;

26

6.2. DATA PREPARATION FOR GAMIT PROCESSING

10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 0.001 Type of Analysis = 1-ITER Data Status = RAW Final ARC = YES Initial ARC = NO Choice of Observable = LC_HELP 0.001 0.001 0.001 0.001 0.001...0.001

; PREFIT / QUICK / 0-ITER/1-ITER/2-ITER ; RAW/CLN YES/NO default = NO for BASELINE/KIINEMATIC, YES for RELAX/ORBIT L1_SINGLE/L1&L2/L1_ONLY/L2_ONLY/LC_ONLY/ L1,L2_INDEPEND./LC_HELP BASELINE/RELAX./ORBIT

; ; ; Choice of Experiment = BASELINE ;

Wobble Constraint = 0.003 0.0001 ; Tighten to 3 mas (10 cm) and .1 mas/day for small network UT1 Constraint = 0.00002 0.0001 ; Tighten to .1 s/day (5 cm/day) for small network Station Error = UNIFORM 10. ; 1-way L1, a**2 + (b**2)(L**2) in mm, ppm, default = 10. 0. Ionospheric Constraints = 0.0 mm + 1.00 ppm Met obs source = GPT 50 ; hierarchical list: RNX ufile GPT/STP [humid value] to match 10.2, use STP 50; new default is GTP 50 Output met = Y ; write the a priori met values to a z-file (Y/N) Use met.list = N Use met.grid = N Use map.list = N Use map.grid = N Zenith Delay Estimation = YES ; YES/NO (default no) Atmospheric gradients = YES ; YES/NO (default no) Gradient Constraints = 0.01 ; gradient at 10 deg elevation in meters; default 0.03 m Interval Zen = 2 Zenith Constraints = 0.50 ; zenith-delay a priori constraint in meters (default 0.5) Zenith Variation = 0.01 100. ; zenith-delay variation, tau in meters/sqrt(hr), hrs (default .02 100.) Zenith Model = PWL ; PWL (piecewise linear)/CON (step) Ambiguity resolution WL = 0.15 0.15 1000. 10. 500. ; FIXDRV, SOLVE > 9.26 defaults Ambiguity resolution NL = 0.15 0.15 1000. 10. 500. ; FIXDRV, SOLVE > 9.26 defaults Yaw Model = Yes ; YES/NO default = YES Antenna Model = ELEV ; NONE/ELEV/AZEL default = NONE Elevation Cutoff = 15 ; Set this in autcln.cmd Tides applied = 31 ; Binary coded: 1 earth 2 freq-dep 4 pole 8 ocean ; 16 remove mean for pole tide ; 32 atmosphere ; default = 31 Use otl.list = N Use otl.grid = Y Etide model = IERS03 ; IERS96/IERS03 Apply atm loading = N ; Y/N for atmospheric loading Use atml.list = N Use atml.grid = N Estimate EOP = 15 SV antenna offsets = Y ; estimate SV antenna offsets (include a priori with sat ICs) Update tolerance = .05 ; minimum adjustment for updating L-file coordinates, default .3 m Use N-file = Y ; Y/N (default no): automatic procedure to reweight by station Delete AUTCLN input C-files = NO ; YES/NO ; default NO Quick-pre observable = LC ; For 1st iter or autcln pre, default same as Choice of observable Decimation factor = 10 ; FOR SOLVE, default = 1 Quick-pre decimation factor = 10 ; 1st iter or autcln pre, default same as Decimation Factor Scratch directory = /tmp ; Directory for scratch files (default /tmp)

Each command is recognized by the keywords at the beginning of the line. They must begin in column one and be spelled out completely and correctly but are not case sensitive. The order of the commands in the les is not important except for the satellite constraints, which must follow the Satellite Constraint keyword. The sestbl. entries can be put into the following categories: 1. Analysis controls. 27

CHAPTER 6. DATA PROCESSING AND ANALYSIS Type of Analysis Data Status Choice of Observable Choice of Experiment 2. Data weighting Station Error Use N-le 3. Ambiguity Resolution Ionospheric Constraints Ambiguity resolution WL Ambiguity resolution NL 4. Atmospheric Parameters Zenith Delay Estimation Interval Zen Zenith Model Zenith Constraints Zenith Variation Atmospheric gradients Gradient Constraints Use met.list Use met.grid Use map.list Use map.grid Met obs source 5. Orbit parameters Initial ARC Final ARC 6. MODEL parameters Tides applied 28

6.3. GENERATING G- AND T-FILES FROM EXTERNAL EPHEMERIDES Use otl.list Use otl.grid Etide model Yaw Model 7. SOLVE parameters Estimate EOP Wobble Constraint UT1 Constraint Decimation factor 8. Cleaning parameters Quick-pre decimation factor Quick-pre observable The most important entries in the sestbl. is to check the number of zenith delay parameters, iterations and the satellite constraints.

6.3 Generating G- and T-les from external ephemerides

Orbital information is input to GAMIT in the form of a tabular ephemeris (T-) le, which contains the positions of all the satellites at 15-minute intervals throughout the observation span, or a G-le of initial conditions which are integrated by arc (in batch run) to create T-le. The Gle can be downloaded from SOPAC archives, however, integration of instantaneous position and velocity values obtained by evaluating the broadcast message parameters or extracted from external Earth-xed ephemeris (in SP 3 format) produced by another software package or an incompatible version of arc will not produce the original accuracy since model differences will cause the orbit to deviate from original orbit as the integration proceeds away from the initial epoch. The most accurate and reliable method of obtaining a starting T-le is to rst download SP 3 le(s) from IGS analysis centre and using the shell script sh_sp3fit as shown below:
sh_sp3fit -f <sp3 files> -d <yr doy1 doy2> -o <orbit name> i <equator prec> -r <radmod>

Where <sp3 files> is the list of SP3 les and <yr doy1 doy2> the year and day-of-year over which the t is to be performed. <orbit name> is the 4-character name for G- and T-les. In this case, we specify IGSF, which means IGS nal orbit. <equator prec> species the inertial reference frame and precession constant for the G- and T-les; J2000 was chosen in this case. <radmod> species the radiation-pressure (non-gravitational force) model to be used with the orbit; the default, BERNE is specied. Part of the G-le for day 001 and year 2004 is shown below: 29

CHAPTER 6. DATA PROCESSING AND ANALYSIS

4 1 12 0 0 GPST J2000 IAU76 BERNE 15 X Y Z XDOT YDOT ZDOT DRAD YRAD BRAD DCOS DSIN YCOS YSIN BCOS BSIN G-file generated by ORBFIT 22- 1-2007 16: 6:49 END PRN 1 -0.1111468080478D+05 0.1032215471537D+05 0.2197619147047D+05 -0.2165070676233D+01 -0.3163242788601D+01 0.3903048496650D+00 0.1030759198283D+01 -0.8027253642367D-02 0.2434678071936D-02 0.4476693138730D-02 0.3526824504122D-02 0.6420264243851D-02 0.5080757147817D-02 0.4238418662426D-02 -0.2578881696451D-02 PRN 2 -0.2616989999127D+05 -0.9043175826147D+03 0.5730355790335D+04 0.6506936986956D+00 -0.2320943372361D+01 0.2987428008352D+01 0.1010574943596D+01 -0.1417044458164D-02 -0.3305724382056D-02 0.4289324919345D-02 -0.2280376317698D-02 -0.2011367355511D-03 0.2421973010893D-02 0.8883030826874D-02 -0.9274838855426D-02

The rst line gives the epoch of the ICs in GPST or UTC (year, day-of-year, hours, minutes, seconds), followed by the time type (GPST or UTC), the inertial frame for the ICs (J2000 or B1950), the precession constant used (IAU76 or IAU68), and the model for non-gravitational ("radiation") accelerations (BERNE). The second line gives the number of ICs plus force-model parameter to be read from the G-le. This is followed by one or more comment lines, terminated by END. The initial conditions for each of the satellites are given as Cartesian components of the position and velocity vectors in units of km and km/s.

6.3.1

Global Files

These les are called global because they can be used for many experiments over the time interval for which they are valid (usually for at least a year). They are usually found in the directory gg/tables in the GAMIT software. These tables are linked to each working directory in order to minimized the amount storage space that will be needed if we store them in each working directory. The global les needed for a single session processing are explained below: Nutations (nutabl.) This is the nutation table and contains nutation parameters in tabular form for transforming between an inertial and Earth-xed system. 30

6.3. GENERATING G- AND T-FILES FROM EXTERNAL EPHEMERIDES Lunar and solar ephemerides (soltab. and luntab.) soltab. is the solar tabular ephemeris, which is the a tabulation of the position of the Earth with respect to the sun. luntab. is lunar tabular ephemeris which contains the Moons position x, y, and z for exact Julian date. Geodetic datums (gdetic.dat) This is the table of paramters of geodetic datums which are specied by the standard ellipsoid parameters, semi-major axis (in meters) and inverse attening, and cartesian offsets (in meters) from the geocenter. leap seconds (leap.sec) This a table of jumps (leap seconds) in TAI-UTC since 1 January 1982. spacecraft, receiver, and antenna characteristics (svnav.dat, antmod.dat, rcvant.dat) svnav.dat gives the correspondence between spacecraft numbers and PRN numbers for each GPS satellite, its mass, and its yaw parameters. The table is updated after each launch or change in yaw status. antmod.dat shows the table of antenna phase center offsets and variations as a function of elevation and azimuth. rcvant.dat is a table of correspondences between 6-character codes and the full names of receivers and antennas used in the RINEX les. Earth rotation (pole., ut1.) pole. is pole table and contains polar motion values in tabular form for interpolation in different GAMIT modules. ut1 is a table which contains TAI-UT1 values. ocean tides (stations.oct and grid.oct) These tables contain ocean tidal loading components of stations and global grid respectively.

6.3.2 Preparing GAMIT to Run

To begin pre-processing for a single day, we execute the script links.com within the working directory to create links directly to the GAMIT global les on the path /gg/tables . All other les necessary for processing is then copied into the working directory. We then execute the script makexp in the working directory to generate most of the additional les needed to complete preprocessing. To run makexp we execute the script as follows:
makexp -expt aauc -orbt IGSF -yr 2004 -doy 001 -nav brdc0010.04n -sinfo 30 00 00 2880

Where aauc is the 4-character experiment code, IGSF is the 4-character orbit type (i.e IGSF for IGS nal), 2004 being the year, 001 is the day-of-year, brdc0010.04n is the navigation le, and 30 00 00 2880 is the session span information sampling interval, hour, minute, and number of epochs respectively. We then run in the following order:
sh_check_sess -sess 001 -type gfile -file gigs4.001 makej brdc0010.04n jbrdc4.001 sh_check_sess -sess 001 -type jfile -file jbrdc4.001 makex aauc.makex.batch

31

CHAPTER 6. DATA PROCESSING AND ANALYSIS

fixdrv daauc4.001

Where the gigs4.001 is the G-le generated by using the script sh_sp3fit as described before, the script sh_check_sess is optional but assures that the satellites requested in session.info are available in the orbital (G- and T-) les. In order to account properly for clock effects in the phase observation, we need to supply extra information regarding the behaviour of the satellite and station clocks. This is done by running the program makej to create the (J-) le of satellite clock values from the navigation message. The output le name is jbrdc4.001 . We again use sh_check_sess to assure that the satellites requested in session.info are available on the J-le. Part of the J-le (jbrdc4.001) generated from the experiment above is given below:
SV clock terms from brdc0010.04n isaacnt MAKEJ 9.96 2006/3/7 09:00:00 YEAR DOY HR MN SEC(UTC) WKNO SOW(GPST) PRN XEAF0 (i4,1x,i4,2i3,1x,f10.7,2x,i4,1x,f14.7,2x,i2.2,2x,3d16.8) 2004 1 1 59 47.0000000 1251 352800.0000000 01 0.32232516D-03 2004 1 3 59 47.0000000 1251 360000.0000000 01 0.32233726D-03 2004 1 5 59 31.0000000 1251 367184.0000000 01 0.32236241D-03 2004 1 5 59 47.0000000 1251 367200.0000000 01 0.32234937D-03 2004 1 7 59 31.0000000 1251 374384.0000000 01 0.32237498D-03 2004 1 9 59 47.0000000 1251 381600.0000000 01 0.32238802D-03 2004 1 11 59 47.0000000 1251 388800.0000000 01 0.32240059D-03 2004 1 13 59 47.0000000 1251 396000.0000000 01 0.32241363D-03 2004 1 15 59 47.0000000 1251 403200.0000000 01 0.32242620D-03 2004 1 17 59 47.0000000 1251 410400.0000000 01 0.32243924D-03 2004 1 19 59 47.0000000 1251 417600.0000000 01 0.32245182D-03 2004 1 21 59 47.0000000 1251 424800.0000000 01 0.32246439D-03 2004 1 23 59 31.0000000 1251 431984.0000000 01 0.32248441D-03 2004 1 23 59 47.0000000 1251 432000.0000000 01 0.32247743D-03 2004 1 1 59 47.0000000 1251 352800.0000000 02 -0.24903379D-03 2004 1 3 59 47.0000000 1251 360000.0000000 02 -0.24908269D-03 2004 1 5 59 47.0000000 1251 367200.0000000 02 -0.24913158D-03 2004 1 7 59 47.0000000 1251 374400.0000000 02 -0.24918001D-03 2004 1 9 59 47.0000000 1251 381600.0000000 02 -0.24922891D-03 2004 1 11 59 47.0000000 1251 388800.0000000 02 -0.24927780D-03 (SunOS) Library ver. 10.59 XEAF1 XEAF2 0.17053026D-11 0.17053026D-11 0.18189894D-11 0.17053026D-11 0.18189894D-11 0.18189894D-11 0.18189894D-11 0.18189894D-11 0.18189894D-11 0.18189894D-11 0.18189894D-11 0.18189894D-11 0.18189894D-11 0.18189894D-11 -0.68212103D-11 -0.68212103D-11 -0.68212103D-11 -0.68212103D-11 -0.68212103D-11 -0.68212103D-11 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00

The rst line is a header constructed by makej , and the second line contains titles for guidance. The third line is a format statement used to read the entries that follow. Each data line contains the coefcients transmitted by the satellite for its own clock. The formula to be used in computing the SV clock offset is given as follows: ts = ts ts = a(0) + a(1) (t t0 ) + a(2) (t t0 )2
(0) (c) (c)

(6.1)

where ts is the time read by the satellites clock and ts is "true" GPS time. The coefcients a , a(1) , a(2) are given in the last three columns (XEAF0, XEAF1, XEAF2) and refer to the reference (c) epoch t0 given in columns 6 and 7 as GPS week number and seconds of week. The numbers in the rst 5 columns give the reference time in GPST. Program makex takes as input the scenario le (session.info), station information le (station.info), satellite clock (J-) le, Broadcast ephemeris (RINEX navigation) le, station coordinates (L-) le, and RINEX observation le, and creates X and K-les for input to fixdrv and model. To run makex , we type only makex aauc.makex.batch where aauc.makex.batch is the name of the control le generated after running the program makexp . The control le contains pointers to the input les and a list of station-days to be processed. The le aauc.makex.batch is given below: 32

6.3. GENERATING G- AND T-FILES FROM EXTERNAL EPHEMERIDES

infor 1 sceno 1 session.info rinex 1 ./ ficaf 0 coord 1 laauc4.001 stnfo 1 station.info xfile 1 ./x svclk 1 ./jbrdc4.001 clock 1 ./k wrorb 0 rdorb 1 ./brdc0010.04n extra 0 site year doy sn sw ver (a4,1x,a4,1x,a3,1x,a1,2x,a3,1x,f4.2) AAUC 2004 001 1 ASH 7.10 onsa 2004 001 1 TRB 2.80 pots 2004 001 1 TRB 3.20

The X-les are the key organizational structure because all X-les for a given session are written with the same start and stop times, selection of satellites, and sampling interval. The process of creating the X-les acts as a lter, catching most of the problems with missing or invalid data, mismatched time tags, and poorly behaved receiver clocks that would cause greater loss of time if discovered later. The K-le generated after running the program makex contains the values of the station clock offset during observation span from pseudorange. The offset from GPS time (or UTC) of each receivers clock must be accounted for in modeling the theoretical value of the phase observations at each epoch. If the positions of the receiver and a satellite are known, along with the offset of the satellites clock, then the pseudorange observation provides a direct measure of the receiver clock offset given by: p1 tr = + ts (6.2) c where is the calculated range to the satellite, p1 is the observed pseudorange and c is the speed of light. An accuracy of about one microsecond in receiver-clock offset is necessary to achieve an accuracy of one millimeter in the estimated baseline vector. In order to achieve this accuracy, the computation is performed using the station and satellite positions (from the L- and T-les) calculated for the theoretical phase observable. In this case, one microsecond (300 m in pseudorange) in the theoretical values is easily achieved when the P-code range is available. In GAMIT, the receiver clock-offset is computed using an average of values calculated from all of the satellites visible at each epoch, detecting and removing anomalous values caused by pseudorange outliers or bad SV clock values. In fact it is not necessary to provide GAMIT with any more information about a receivers clock than that incorporated in the pseudoranges at each epoch. One of the reasons why it is necessary to generate a more explicit model of the receiver-clock behavior at an earlier stage in the processing, is to provide a way of detecting poor receiver performance. Part of the K-le kaauc4.001 obtained is shown below:
AAUC AAUC AAUC AAUC AAUC AAUC AAUC AAUC 25 4 24 1 20 13 27 25 2003 2003 2003 2003 2003 2003 2003 2004 365 365 365 365 365 365 365 1 23 23 23 23 23 23 23 0 59 59 59 59 59 59 59 1 47.0000 47.0000 47.0000 47.0000 47.0000 47.0000 47.0000 47.0000 0.08279178 0.07059966 0.07715355 0.06865501 0.07496044 0.06889270 0.08108126 0.08294009 0.00005798 -0.00003471 0.00008678 0.00032231 -0.00019668 -0.00003172 0.00087292 0.00005798 0.00030194 0.00024291 0.00031915 0.00020849 0.00018628 0.00020822 0.00018479 0.00022382

33

CHAPTER 6. DATA PROCESSING AND ANALYSIS

AAUC AAUC AAUC AAUC AAUC AAUC AAUC AAUC AAUC AAUC AAUC AAUC AAUC AAUC 4 24 1 20 13 27 25 4 24 1 20 13 27 25 2004 2004 2004 2004 2004 2004 2004 2004 2004 2004 2004 2004 2004 2004 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 3 3 3 3 3 3 3 5 47.0000 47.0000 47.0000 47.0000 47.0000 47.0000 47.0000 47.0000 47.0000 47.0000 47.0000 47.0000 47.0000 47.0000 0.07052486 0.07687638 0.06864850 0.07511012 0.06870177 0.08073944 0.08309097 0.07045608 0.07660291 0.06864661 0.07526239 0.06851576 0.08039855 0.08324433 -0.00003472 0.00008678 0.00032231 -0.00019668 -0.00003172 0.00087292 0.00005798 -0.00003472 0.00008678 0.00032231 -0.00019668 -0.00003172 0.00087292 0.00005798 0.00016219 0.00023844 0.00012967 0.00010677 0.00012898 0.00010483 0.00014568 0.00008149 0.00015772 0.00005085 0.00002727 0.00004974 0.00002490 0.00006753

The rst column gives the station code, followed by the satellite PRN number, year, day of year, hours, minutes, and seconds (always UTC). The eighth column is the observed pseudorange to the satellite at the time given, in units of seconds, followed by the offset of the satellite clock (from GPS time) computed from the transmitted clock corrections. The nal number is the receiver clock correction computed using (6.2). Any differences in receiver clock values computed using the data from different satellites at the same time are due to errors in the pseudorange measurements, the satellite clock models, or the geometrical models (station coordinates and satellite ephemeris). For instance, in kaauc4.001 the corrections computed at 00h 01m 47s using data from satellites 1 and 13 differ by up to 0.7 microseconds, equivalent to radial position error of 210 meters. This must be probably due to errors in the broadcast ephemeris. Finally, program fixdrv is run using the D-le daauc4.001 generated from program makexp to read the analysis controls, and creates a batch le for GAMIT processing. D-le contains the names of the satellite ephemeris (T-), clock (I- and J-), and observation (X- or C-) les to be used. fixdrv also take as inputs sestbl. and sittbl. les. The D-le, daauc4.001 is shown below:
1 1 laauc4.001 tigsf4.001 iaauc4.001 jbrdc4.001 3 xaauc4.001 xonsa4.001 xpots4.001

The integer on the rst line gives the number of independent solutions to be performed; if greater than one, it indicates that you have concatenated essentially independent runs for serial processing. The integer on the second line indicates the number of sessions in this solution. Lines from 3 to 6 contain the name of the coordinate (L-) le to be read, the ephemeris (T-) le to be read, and the name of the I-le to be read or created. The I-le is optional and can be left blank or set to NONE if we dont have K-les readily available from which to calculate clock rates. The integer in line 7 shows the number of stations and the last three lines contain the names of the X-les to be read. The output of fixdrv is primary B-le which contains a sequence of generated secondary B-les for controlling the batch mode of data processing. It invokes the GAMIT modules in the a appropriate order for the requested type of analysis. An example of a primary B-le baauc4.bat generated from daauc4.001 is as shown below: 34

6.3. GENERATING G- AND T-FILES FROM EXTERNAL EPHEMERIDES

#!/bin/csh/ # Remove any existing GAMIT.status, .warning, .fatal files /bin/rm GAMIT.status >& /dev/null /bin/rm GAMIT.warning >& /dev/null /bin/rm GAMIT.fatal >& /dev/null # # Generation of yaw file yawtab yigsf4.001 tigsf4.001 yigsft.001 30 # # Initial O-Cs and editing model < baauc4.001 model < baauc4.002 model < baauc4.003 model < baauc4.004 sh_autedit -base autcln autcln.cmd.prefit + daauc4.001 4 >! autcln.out # # Quick solutions cfmrg < baauc4.006 solve < baauc4.007 set rms = grep Pre qaauc1.001 | awk $6 > 1 {print "no";exit} if ( $rms == "no" ) echo Quick solution rms is too high if ( $rms == "no" ) scandd maauc4.001 if ( $rms == "no" ) exit model < baauc4.008 model < baauc4.009 model < baauc4.010 model < baauc4.011 sh_autedit -base autcln autcln.cmd.prefit + daauc4.001 b >! autcln.out # # Full solution cfmrg < baauc4.012 solve < baauc4.013 # # Check the quality of the final solution set rms = grep "Postfit nrms" qaauc4.001 | awk $6 > 1 {print "no";exit} if ( $rms == "no" ) echo Full solution rms is too high if ( $rms == "no" ) scandd maauca.001 if ( $rms == "no" ) exit

6.3.3 Running GAMIT

To execute the analysis run in the foreground, we type: csh baauc4.bat To run in the background, we type: gbat baauc4.bat All the required modules will then be executed automatically to perform the processing sequence as dened in the batch le. As each module runs, it writes into les GAMIT.status, GAMIT.warning and GAMIT.fatal messages recording the progress of the run, allowing for monitoring in order to determine where problems arose if occurred. 35

CHAPTER 6. DATA PROCESSING AND ANALYSIS

6.3.4

Output les after running GAMIT

A lot of les are generated running GAMIT successfully, but the main output les which are of importance so far as this project is concerned are Q-le which contain the record of the analysis after running solve, autcln.sum which is the auto clean summary le, and the H-le which contains the covariance matrix and parameter adjustments for solution generated with loose constraints. H-le is the main le used as input to GLOBK. A section of the Q-le is shown below:
Total parameters: 168 live parameters: 84 Prefit nrms: 0.13375E+01 Postfit nrms: 0.25184E+00 Total parameters: 168 live parameters: 84 Prefit nrms: 0.13375E+01 Postfit nrms: 0.25184E+00 -- Uncertainties not scaled by nrms Label (units) 1*AAUC GEOC LAT 2*AAUC GEOC LONG 3*AAUC RADIUS 4*ONSA GEOC LAT 5*ONSA GEOC LONG 6*ONSA RADIUS 7*POTS GEOC LAT 8*POTS GEOC LONG a priori N56:50:16.19641 E009:59:14.14762 6363.1889293661 N57:13:13.30635 E011:55:31.86251 6363.0439733452 N52:11:35.04261 E013:03:57.92879 Adjust (m) 0.1081 0.1009 0.0113 -0.0018 0.0014 -0.0038 0.0020 -0.0014 Formal 0.0050 0.0051 0.0175 0.0042 0.0043 0.0092 0.0042 0.0043 Fract Postfit 21.6 N56:50:16.19990 19.7 E009:59:14.15359 0.6 6363.18894065 -0.4 N57:13:13.30629 0.3 E011:55:31.86259 -0.4 6363.04396957 0.5 N52:11:35.04267 -0.3 E013:03:57.92871

dms dms km dms dms km dms dms

Baseline vector (m ): AAUC (Site 1) to ONSA (Site 2) X -57327.21471 Y(E) 108216.35809 Z 22997.98041 L 124603.75920 +- 0.00917 +- 0.00351 +- 0.01396 +- 0.00308 Correlations (X-Y,X-Z,Y-Z) = 0.47287 0.90012 0.47634 N 44135.90239 E 116518.70021 U -1229.40179 L 124603.75920 +- 0.00342 +- 0.00304 +- 0.01644 +- 0.00308 Correlations (N-E,N-U,E-U) = -0.00757 0.02457 0.07218 Baseline vector (m ): AAUC (Site 1) to POTS (Site 3) X 372703.86995 Y(E) 278416.60188 Z -297997.65736 L 552473.15108 +- 0.01078 +- 0.00582 +- 0.01469 +- 0.00513 Correlations (X-Y,X-Z,Y-Z) = 0.32748 0.84608 0.31046 N -510630.22057 E 209559.80331 U -23833.78585 L 552473.15108 +- 0.00492 +- 0.00542 +- 0.01767 +- 0.00513 Correlations (N-E,N-U,E-U) = -0.10587 -0.13980 0.05534 Baseline vector (m ): ONSA (Site 2) to POTS (Site 3) X 430031.08466 Y(E) 170200.24380 Z -320995.63777 L 562968.07745 +- 0.00782 +- 0.00507 +- 0.00989 +- 0.00463 Correlations (X-Y,X-Z,Y-Z) = 0.22913 0.71635 0.24717 N -557035.93081 E 77665.21644 U -24741.50698 L 562968.07745 +- 0.00467 +- 0.00486 +- 0.01180 +- 0.00463 Correlations (N-E,N-U,E-U) = -0.02688 -0.14709 0.02121

(meters)

(meters)

(meters)

(meters)

(meters)

(meters)

6.3.5

Evaluating the solutions

The primary indicator used in evaluating the quality of the solution is the "Postt nrms" which GAMIT writes to the Q-le after solve has been executed. If the data were randomly distributed and the a priori weights were correct, the solution usually produces a nrms of about 0.25. Anything larger than 0.5 means that there are cycle slips that have not been removed or 36

6.4. AUTOMATIC BATCH PROCESSING associated with extra bias parameters or that there is a serious modeling problem. If the nal solution of a batch sequence meets this criteria, there is usually no need to look carefully at any other output, though the rms of residuals in autcln.sum.post will show the relative quality of stations in the network.

6.4 Automatic Batch Processing

In the previous sections we have been discussing how to process a single session of data. In this section we focus on multiple session data processing. Once we have an insight of how to process a single session, understand the le structure of GAMIT, time can be saved by processing signicant data by using the automatic batch processing script sh_gamit. sh_gamit takes you, with a single command, from raw or RINEX data over a range of days to a solution and sky plots of phase data as a record of the GAMIT analysis. The only preparation required is setting up the control les, most of which are common to all analyses of a particular era, and assembling of data in one or more directories on your system. The rst step is to create an experiment directory. Parallel to this directory, we have the following directories: tables tables and templates brdc contains precise ephemeris glbf Binary h-les gplot For running plot of coordinate time series igs Precise igs orbits rinex Rinex les doy each day of the year processing directories We then execute the script sh_setup from the experiment directory which will invoke links.tables to link into experiment ./tables directory all of the standard data tables and also the following control and data les: Process.defaults: This le is edited to specify the computational environment, sources for internal and external data and orbit les, start time and sampling interval, and instruction for archiving results. The process.defaults le used in this experiment is as shown below:
## LOCAL DIRECTORIES # Directory for translation of raw data set rawpth = "" # Directory path for raw archives (search all levels); e.g. /data18/simon set rawfnd = "" # Input files for RINEX translators set mpth = "" # RINEX files directory set rpth = "$procdir/rinex" # Directory path for RINEX archives (search all levels); e.g.

37

CHAPTER 6. DATA PROCESSING AND ANALYSIS

/data18/simon set rnxfnd = "" # Broadcast orbit directory set bpth = "$procdir/brdc" # IGS files directory set ipth = "$procdir/igs" # G-files directory set gpth = "$procdir/gfiles" # GAMIT and GLOBK tables directory set tpth = "$procdir/tables" # Output gifs directory set gifpth = "$procdir/gifs" # Globk solution directory set glbpth = "$procdir/gsoln" # Globk binary h-file directory set glfpth = "$procdir/glbf" # Template files set templatepth = "procdir/templates" # Place to store temporary control files set cpth = "$procdir/control" # Archive root directory (cannot be null) set archivepth = "$procdir/archive"

## FTP INFO FOR REMOTE FILES # Raw data archive # set rawarchive = "" # set rawdir = "" # set rawlogin = "" # Addresses for CDDSI, SOPAC, IGSCB, UNAVCO, BKG, IGN, USNO are given in template/ftp_info ##GAMIT # Set sampling interval, number of epochs, and start time for processing set sint = 30 set nepc = 2880 set stime = 0 0 # Variables for updating tables set stinf_unique = "-u" set stinf_nosort = "-nosort" set stinf_slthgt = "2.00" # Set "Y" to use RINEX header coordinates not in lfile or apr file set use_rxc = "N" # Broadcast orbits set brdc = jplm # Minimum x-file size to be processed (Def. 300 blocks; most OS use 1 Kb blocks) set minxf = 300 # Set search window for RINEX files which might contain data for day - default check the previous day set rx_doy_plus = 0 set rx_doy_minus = 0 # Default globk .apr file set aprf = itrf00.apr # Set compress (copts), delete (dopts) and archive (aopts) options. (Dont forget to set the archivepth.) # Possible d-, c-, and aopts: D, H, ao, ac, as, b, c, d, e, g, h, i, j, k, l, m, o, p, q, t, x, ps, all" set dopts = ( c ) set copts = ( x k ao ) set aopts = # Set the rinex ftp archives (defined in ftp_info) you would like to look for data in. # (Default archives searched are: sopac, cddis and unavco). set rinex_ftpsites =""

38

6.4. AUTOMATIC BATCH PROCESSING

## RESOURCES # Minimum raw disk space in Kbytes set minraw = 0 # Minimum rinex disk space in Kbytes set minrinex = 300000 # Minimum archive disk space in Kbytes set minarchive = 0 # Minimum working disk space in Kbytes set minwork = 300000

Sites.defaults: We edit this le to specify the IGS stations which are to be used in the experiment and how the station log data are to be handled. Even though sh_gamit can automatically download RINEX and ephemeride les from SOPAC and precise orbits from IGS ftp-servers, it much simple when all the data is in the local system in each directory.

6.4.1 Processing Data from AAUC Site

In this subsection we describe the data processing from AAUC site and some of the results obtained. An experiment directory is created for each year and the rst 30-doy of each experiment is processed. Since all the experiments follow the same processing procedure, it will be a repetition to describe all, we therefore describe how year 2005 was processed even though the results for the other years will be shown. We create a directory called expt_2005 with a structure as described previously. In the rinex directory we execute the script:
sh_get_rinex -archive sopac -yr 2005 -doy 001 -ndays 30 -sites onsa pots

The above script download year 2005 RINEX les of IGS sites ONSA and POTS from SOPAC archives retrieving doy 001 to 30. In the case of AAUC, since the data is directly stored on the university server, we download as follows:
lftp http://gps.aau.dk/aauc/data/2005/

Which takes us to the ftp site. We then start downloading the data as follows:
mget AAUC00* mget AAUC01* mget AAUC02* mget AAUC0300.05.zip

The les downloaded from SOPAC are uncompressed and that of AAUC is unzipped and the ephemeris data is discarded. It must be emphasized that the RINEX le names of AAUC are in uppercase (i.e. AAUC) and during the single session processing, somehow, GAMIT fails to recognize it as a RINEX le but when the name was changed to a lowercase, it was processed. Therefore all the RINEX le names are changed into lowercases. In the igs directory, we execute the script 39

CHAPTER 6. DATA PROCESSING AND ANALYSIS

sh_get_orbits -archive sopac -center igs -yr 2005 -doy 001 -ndays 30 -ftp_prog ncftp

The above script downloads precise IGS orbits from SOPAC archives with year of 2005 from doy 001 to 030 using ftp program called ncftp. The data is uncompressed, we then run the script sh_sp3fit on each of the les to get G- and T-les. In the brdc directory we execute the script
sh_get_nav -archive sopac -yr 2005 -doy 001 -ndays 30 -ftp_prog ncft

This script downloads precise ephemerides from SOPAC archives of year 2005, from doy 001 to 030. The data is then uncompressed. In each doy directory we execute the script links.day to link to tables all the necessary tables and les needed for batch processing. sh_link_rinex2 is also executed from each doy directory to link to the rinex directory all the rinex les of each site and day needed for the data processing. We then start the batch processing for daily solution from the expt_2005 directory by executing the script
sh_gamit -expt aauc -s 2005 001 030 -orbit IGSF -noftp -aprfile itrf00.apr >&!sh_gamit.log

With start and stop days set to 1 to 30 of year 2005 specifying the use of nal ephemerides from IGS (i.e. IGSF). It writes to the screen every a record of each step which is redirected to a log le called sh_gamit.log. This is to enable identify the point and reason for failure should that occur. All the solutions and les are stored in the doy directories for analysis.

6.5

Data Preparation for GLOBK

The main input le for GLOBK is the GAMIT h-les which can be found in the doy directory after running sh_gamit successfully. All the ascii h-les are put in the the directory glbf described in the previous section. Firstly, we convert all the ascii h-les into binary h-les that can be read by GLOBK. This is accomplished via the program htoglb. Next, we run glred for all the binary h-les for the 30-day period to obtain time series of station coordinates, which we then plot and examine for outliers and appropriate scaling to obtain reasonable uncertainties. We then remove these outliers to obtain clean data set and then we repeat the processing by running globk instead of glred to combine the daily h-les into a single h-le that represent the estimate of the station positions. 40

6.5. DATA PREPARATION FOR GLOBK

6.5.1 Running Glred

The input les required to run glred are the binary h-les, le containing the a priori coordinates and velocities for the stations (itrf00.apr), table of the Earth orientation parameter values for the times of the data (pmu_bull_b), and globk and glorg command les, and the le containing the list of binary h-les to be processed with a gdl extent (global directory list). All these les are put in the gsoln directory for processing. We copy from the directory /af s/ies.auc.dk/sw/pack/gamit10.21/src/example/templates globk_comb.cmd and glorg_comb.cmd which are the globk and glorg command les respectively. We then edit this templates and set the constraints for the processing. Part of the globk_com.cmd command le is shown below:
* * Globk file for combination of daily h-files * x eq_file ../tables/renames x make_svs ../tables/sat1.apr com_file globk_comb.com x srt_file glb.srt x sol_file globk_comb.sol # x earth-rotation values (not needed if pmu free in final combination; pmu.bull_b begins only at 1 Jan 1992) in_pmu /afs/ies.auc.dk/group/07gr1049/expt_2005/tables/pmu.bull_b

# apr site file(s) apr_file /afs/ies.auc.dk/group/07gr1049/expt_2005/tables/itrf00.apr # Optionally use separate file for sites to be used and random or Markov noise x source ../tables/globk.uselist # (1) Max chi**2, (2) Max prefit diff, (3) Max rotation; defaults are 100 10000 10000 max_chi 30. 10000. 10000

* Apply the pole tide whenever not applied in GAMIT app_ptid ALL # Allow the network to be loose since using glorg for stabilization apr_neu all 100 100 100 0.0500 0.0500 0.0500

# Satellites are loose if combining with global SOPAC H-files # X Y Z XDOT YDOT ZDOT DRAD YRAD ZRAD BRAD XRAD DCOS DSIN YCOS YSIN BCOS BSIN apr_svs all 100 100 100 10 10 10 1 1 .02 .02 .02 . 02 .02 . 02 . 02 .02 . 02 # tight if not combining with global data apr_svs all .05 .05 .05 .005 .005 .005 .01 .01 F F F F F F F F F F F F F F F # Keep EOP loose x apr_wob 100. 100. 10. 10.0 0.0 0.0 x apr_ut1 100. 10. 0.0 0.0 0.0 0.0 x mar_wob 36500 36500 365 365 0 0 0 0 x mar_ut1 36500 365 0 0 0 0 # unless not using global data apr_wob .25 .25 .001 .001 0 0 0 0 apr_ut1 .25 .25 .001 .001 0 0 mar_wob 1 1 .001 .001 mar_ut1 1 .001

0.0

0.0

.0005 m**2/yr = 15 mm/half-yr * Estimate translation apr_tran .005 .005 .005 0 0 0

41

CHAPTER 6. DATA PROCESSING AND ANALYSIS

x mar_tran .0025 .0025 .0025 0 0 0 # Set minimal globk print options since using glorg output prt_opt blen brat vsum # Invoke glorg for stabilization org_cmd glorg_comb.cmd org_opt cmds psum gdlf blen vsum # Write out an h-file if needed for future combinations out_glb [email protected]

Part of the glorg_comb.cmd command le is shown below:

* * # # # x # Glorg file for repeatabilities apr site file(s) ITRF00 for global stabilization North Amercian frame apr_file /afs/ies.auc.dk/group/07gr1049/experiment/tables/itrf00.apr Regional stabilization for filtering spatially-correlated errors apr_file /afs/ies.auc.dk/group/07gr1049/expt_2005/tables/itrf00.apr

# Define the stabilization frame x source ../tables/stab_site.global x source /afs/ies.auc.dk/group/07gr1049/expt_2004/tables/stab_site.scand # Define the stabilization frame stab_site clear onsa_gps pots_gps

# Set parameters to estimate in stabilization x pos_org xrot yrot zrot xtran ytran ztran x no rotation if regional stabilzation pos_org xtran ytran ztran # Set height ratios x cnd_hgtv 10 10 2.0 2.0 x downweight heights in stablization x cnd_hgtv 1000 1000 2.0 2.0 # Iterations and editing x stab_ite 4 0.8 4.

Running daily solution for glred is not necessary since the data spans over several days and the analysis is performed for rst 30-day of each of the chosen years, drastic variations in coordinates and velocities are not expected just in a day. Therefore a strategy was adopted to process 5-day averages before computing deformation velocities, this implies that 6 gdl les were created, each containing 5 binary h-les to be processed for each year after which we combine the rst two 5-day averages and then scaling to the required weight in the gdl le. We then run glred as follows:
glred 6 globk_rep.prt globk_rep.log aauc05_fix.gdl globk_comb.cmd

where the value 6 means that the input should be sent to the current window, globk_rep.prt is the output print le with the solution in it, globk_rep.log is the log le containing the 42

6.6. RESULTS running time for the program and the pret 2 per degree of freedom value for each input matrix le, aauc05_fix.gdl is the gdl le for the year 2005, and globk_comb.cmd is the globk command le used to control the solution.

6.6 Results
After running the glred, two versions of the solution les are generated, one from the globk solution (the globk_rep .prt le in the globk command-line arguments) which contains summary of the nal solution of velocities and position estimates, baseline components and their sigmas as described earlier and one from the glorg solution (the globk_rep.org le in the globk command le). Some parts of the globk_rep.prt le obtained is shown below:
Globk Analysis SUMMARY VELOCITY ESTIMATES FROM GLOBK Ver 5.11S Long. Lat. E & N Rate E & N Adj. (deg) (deg) (mm/yr) (mm/yr) 13.066 52.380 20.74 -11.65 1.41 -25.89 11.926 57.396 28.33 16.65 11.08 3.06 9.987 57.014 -11.38 21.61 -11.38 21.61 GLOBK: BASELINE LENGTHS BASELINE POTS_GPS to ONSA_GPS POTS_GPS to AAUC_GPS ONSA_GPS to AAUC_GPS

E & N (mm/yr) 47.79 37.33 35.82

+-

RHO

37.95 0.100 34.34 -0.013 35.37 0.079

H Rate (mm/yr) -7.57 3.22 8.41

H adj. -6.27 0.68 8.41

+-

SITE

47.45 POTS 46.86 ONSA 46.81 AAUC

Length (m) 562966.02957 552471.03982 124603.10935

Adjust (m) -2.05201 -2.05102 -0.78300

Sigma (m) 0.00531 0.00524 0.00162

GLOBK: BASELINE COMPONENTS Baseline North East Est. Adj. +Est. Adj. +(mm) (mm) ONSA_GPS-POTS_GPS 558377728.4 -1004.9 25.8 -172551214.3 -548.0 6.1 -0.527 AAUC_GPS-POTS_GPS 515911933.7 -886.8 24.0 -282597970.1 -903.1 10.7 -0.803 AAUC_GPS-ONSA_GPS -42465794.7 118.2 3.4 -110046755.8 -355.1 6.1 0.565 GLOBK: BASELINE COMPONENT RATES OF CHANGE Baseline Length Est. Adj. +(mm/yr) ONSA_GPS-POTS_GPS AAUC_GPS-POTS_GPS AAUC_GPS-ONSA_GPS 26.42 42.25 35.33 27.06 49.15 14.56 38.58 38.78 24.63

Rne Est. -95947.4 -81299.5 14647.9 Adj. (mm) 2878.9 2672.5 -206.4

Est.

North Adj. +(mm/yr) 28.94 47.50 18.56 38.97 41.35 29.81

Est.

East Adj. +(mm/yr) 9.67 -12.79 -22.46

Rne Est.

Adj. (mm/yr) 6.95 14.68 7.73

28.30 33.27 4.96

7.60 -32.12 -39.71

59.72 0.099 56.76 0.253 24.43 -0.045

10.79 15.98 5.19

Examining the globk output is useful mainly if the glorg output indicates a problem with the solution and you want to determine if the source is in the data or the constraints. For the year 2000, we obtained plots of coordinate repeatabilities as shown in gures 6.1, 6.2, and 6.3 For the year 2001, we obtained plots of coordinate repeatabilities as shown in gures 6.4, 6.5, and 6.6 For the year 2004, we obtained plots of coordinate repeatabilities as shown in gures 6.7, 6.8, and 6.9 43

CHAPTER 6. DATA PROCESSING AND ANALYSIS

Figure 6.1: 5-day variation of the north-coordinate in year 2000.

Figure 6.2: 5-day variation of the east-coordinate in year 2000.

44

6.6. RESULTS

Figure 6.3: 5-day variation of the up-coordinate in year 2000.

Figure 6.4: 5-day variation of the north-coordinate in year 2001.

45

CHAPTER 6. DATA PROCESSING AND ANALYSIS

Figure 6.5: 5-day variation of the east-coordinate in year 2001.

Figure 6.6: 5-day variation of the up-coordinate in year 2001.

46

6.6. RESULTS

Figure 6.7: 5-day variation of the north-coordinate in year 2004.

Figure 6.8: 5-day variation of the east-coordinate in year 2004.

47

CHAPTER 6. DATA PROCESSING AND ANALYSIS

Figure 6.9: 5-day variation of the up-coordinate in year 2004. For the year 2005, we obtained plots of coordinate repeatabilities as shown in gures 6.10, 6.11, and 6.12 For the year 2006, we obtained plots of coordinate repeatabilities as shown in gures 6.13, 6.14, and 6.15 It can seen from the gures that there are some outliers in the plots of the years 2000, 2004, 2005, and 2006. But we do not jump right to delete them, we must rst have to assess what might be causing the problem. Probably they might be due to poor stabilization, or perhaps theres some poor bias xing that is degrading the results. This can also happen when constraints set in the command les are too tight. It must be emphasized that throwing out little bit of data may be signicant in estimating good station velocity and position since the data span is only the rst 30 days of each of the year specied in the experiment. In the gdl le, we comment on the the binary h-les causing the outliers and then we run globk to obtain the nal velocity of the station. These results will be shown later in the report. Table 6.1 shows the summary of the velocities of AAUC obtained in each year. The values of the velocities do not make any conclusive analysis. For instance, the Up rate in the year 2007 shows subsidence deformation, while that of 2005 show some uplift. This might be due to the fact that, there are some outliers which is degrading the results. These outliers must be removed before making any form of analysis on the nal results. Details of this analysis will be discussed in the next chapter.

48

6.6. RESULTS

Figure 6.10: 5-day variation of the north-coordinate in year 2005.

Figure 6.11: 5-day variation of the east-coordinate in year 2005.

49

CHAPTER 6. DATA PROCESSING AND ANALYSIS

Figure 6.12: 5-day variation of the up-coordinate in year 2005.

Figure 6.13: 5-day variation of the north-coordinate in year 2006. Table 6.1: Summary of velocities of site AAUC North rate (mm/yr) East rate (mm/yr) Up rate (mm/yr) -0.38 6.39 0.54 -4.33 32.84 -3.05 -13.79 25.05 -7.05 21.61 -11.38 8.41 16.48 -0.34 3.15 50

Year 2000 2001 2004 2005 2006

6.6. RESULTS

Figure 6.14: 5-day variation of the east-coordinate in year 2006.

Figure 6.15: 5-day variation of the up-coordinate in year 2006.

51

Chapter

Land Uplift in Fennoscandia

Fennoscandia, which comprises Norway, Sweden, and Denmark, has undergone numerous deformations through tectonic activities. In order to really understand the spatial characteristics of the present day land uplift and/or subsidence, knowledge of the past geodynamics will be very helpful. Figure 7.1. shows the map region covering Fennoscandia

Figure 7.1: Map of Fennoscandia. The land uplift, postglacial rebound (PGR) or glacial isostatic adjustment (GIA) (now commonly termed the latter) process in Fennoscandia has been a subject of scientic research for more than a century. It is now recognized to be part of the global process of GIA, which originates from the last glacial cycle culminating about 20,000 years ago. When the load from the ice (thickness of about 23 km) was removed, the Earth responded as a viscoelastic body, resulting in vertical and horizontal displacements towards a new equilibrium. A wide variety of geophysical and geodetic observable features and quantities have been exploited to study the GIA process. This includes tide-gauge records, repeated geodetic leveling, gravity anomalies, changes in gravity, and time series of ancient sea level elevations. These are, in a geometrical 52

7.1. BASELINE INVESTIGATION sense, primarily related to the vertical component. In recent years, especially in the early 1990s there has been a rapid development in space geodesy, evolving predominantly around GPS. In 1993, the Baseline Inferences for Fennoscandian Rebound Observations Sea Level and Tectonics (BIFROST) project was started, with a primary goal to establish a new and useful 3D measurement of movements of the Earths crust in this region, able to constrain models of the GIA process in Fennoscandia.[LJSD06]

7.1 Baseline Investigation

As mentioned earlier in the report two IGS stations were chosen as stabilization sites. These stations are POTS and ONSA. Figure 7.2 shows the location of the stations with respect to the analysis site of AAUC.

Figure 7.2: Map showing the stabilization stations. It can be seen from gure 7.2 that the baseline from AAUC to ONSA is shorter than AAUC to POTS and this produces a smaller sigma as it is evident in table 7.1 obtained from the processing of data in each year. Therefore all analysis will be with reference to ONSA. Also the station POTS is outside the region of Fennoscandia but was included in as reference stations to have two IGS stations as stabilization sites. This is to accept the idea that having a redundancy in choosing stabilization sites is better. GIA is expected to cause some relative land uplift between the sites AAUC and ONSA. This uplift should be evident in the baselines between the two stations, but from the gures and tables presented in this section, we cannot make any conclusive analysis as regarding possible 53

CHAPTER 7. LAND UPLIFT IN FENNOSCANDIA Table 7.1: Baselines and their sigmas and rates in 2000 Baselines Baseline Length [m] Sigma [mm] Rates of Change [mm/yr] ONSA to POTS 562965.99153 6.50 -2.69 AAUC to POTS 552470.99901 6.44 -11.92 AAUC to ONSA 124603.09579 2.73 10.48

Figure 7.3: 5-day variation of the baseline in year 2000. Table 7.2: Baselines and their sigmas and rates 2001 Baselines Baseline Length [m] Sigma [mm] Rates of Change [mm/yr] ONSA to POTS 562966.03426 5.17 25.22 AAUC to POTS 552471.03951 5.12 -15.85 AAUC to ONSA 124603.10661 2.19 -34.60 Table 7.3: Baselines and their sigmas and rates in 2004 Baselines Baseline Length [m] Sigma [mm] Rates of Change [mm/yr] ONSA to POTS 562965.98447 5.87 33.03 AAUC to POTS 552470.99370 5.82 -22.48 AAUC to ONSA 124603.09624 2.48 -15.08 Table 7.4: Baselines and their sigmas and rates in 2005 Baselines Baseline Length [m] Sigma [mm] Rates of Change [mm/yr] ONSA to POTS 562966.02957 5.31 26.42 AAUC to POTS 552471.03982 5.24 42.25 AAUC to ONSA 124603.10935 1.62 35.33 54

7.1. BASELINE INVESTIGATION

Figure 7.4: 5-day variation of the baseline in year 2001.

Figure 7.5: 5-day variation of the baseline in year 2004.

Table 7.5: Baselines and their sigmas and rates in 2006 Baselines Baseline Length [m] Sigma [mm] Rates of Change [mm/yr] ONSA to POTS 562965.90979 5.97 18.68 AAUC to POTS 552470.91643 5.89 25.70 AAUC to ONSA 124603.07114 1.75 15.44 55

CHAPTER 7. LAND UPLIFT IN FENNOSCANDIA

Figure 7.6: 5-day variation of the baseline in year 2005.

Figure 7.7: 5-day variation of the baseline in year 2006.

56

7.1. BASELINE INVESTIGATION differential motions especially with respect to land uplift as result of GIA in the Fennoscandia region. This is because the land uplift rates in this region is estimated to be within a fraction of a millimeter per year [LJSD06]. For instance results for the year 2000 of the baseline rates of AAUC-ONSA show some land uplift while that of the year 2001 shows some subsidence. The inconsistencies in rates of change of the baselines might be due to the fact that there is insufcient data to determine accurately the displacements. Ideally monitoring such displacements require campaigns for a long period of time ( 10 years) to get reliable estimates of vertical movement rates. [DCWM05]. It might also be caused by some systematic effect such as the GPS monument at AAUC. Though the building on which the GPS antenna is mounted has been installed long enough to have settled, we cannot a priori preclude some residual monument instability and a possible thermal dilatation effect. As a matter of fact, the latter factor is more plausible. But these effects contribute only few millimeters in error to the overall results.

57

Chapter

Final Results and Discussion

This chapter presents the nal results obtained after erroneous data have been identied and removed, its interpretation, signicance and implications will be thoroughly discussed, and spatially displaying the results. Figures 8.1, 8.2, and 8.3 display the nal results of the North, East, and Up rates respectively of the site AAUC.

Figure 8.1: Yearly variation of the north coordinate. These results were obtained from the nal globk analysis. It can be seen from the gures that nrms and the wrms of the nal results are unusually high even though the rates are very encouraging. One of the reasons is due to the fact that, the parameters set in the nal globk command le are loosely constrained. Also, although the zenithal tropospheric path delays are treated as unknowns by gamit, their mismodelling in some stages of the processing could also involve errors of a few millimetres in the up component of the site. Indeed, the zenithal delays are estimated from slant-range delays, which may be strongly biased by specic meteorological conditions, in particular for low elevation signals. When these are carried to the globk processing stage, they cannot be resolved since globk assumes a linear model and these errors could propagate to the nal solution. Another possibility of error could be inherent in data editing be it weighting, deleting, etc. It is also known from the comparison of time series of permanent GPS sites computed by various centres that the analysis noise in a general sense may produce an articial, steadily increasing or decreasing velocity component [DCWM05]. These errors are not accounted for 58

Figure 8.2: Yearly variation of the east coordinate.

Figure 8.3: Yearly variation of the up coordinate.

59

CHAPTER 8. FINAL RESULTS AND DISCUSSION in my nal results since it is outside the scope of this project to model signal noise affecting velocity of GPS sites. The up component rate shows some subsidence since the estimated deformation velocity is negative. The east component rate is also negative implying that site is moving westwards with respect to its actual position. Since north rate is positive, therefore we can sum all and say that GPS resultant velocity is in north-west direction with subsidence. Finally we spatially show the direction of the movements of the site AAUC relative to the other IGS sites used in the experiment of this project.

Figure 8.4: Map displaying the direction of movement of sites. Figure 8.4 shows the direction of movements of site AAUC. It can be seen that the site is moving almost in the same direction as the site ONSA.

60

Chapter

Conclusion and Recommendations

The processing and analysis of the GPS data used for the project has demonstrated that the GPS reference station at AAUC is in continuous motion as a result of the Fennoscandian glacial isostatic adjustment (GIA) process, based on continuous GPS observations obtained using only the rst 30 days of data of January 2000, 2001, 2004, 2005, and 2006. But the amount of data used for the analysis is insufcient to determine accurately the magnitude and direction of motion of the site. This is evident in the high rate of the northern and easting values of the nal results obtained. The poor repeatability of AAUC with respect to its high northern and eastern rates estimate of this project makes it unreliable for looking for accurate tectonic displacements.

9.1 Recommendations
Due to the fact that during this project certain hinderance were encountered, it is very imperative that they are mentioned so that if further investigation is required in the near future, having fore knowledge of these impediments could be very valuable assets. I therefore recommend the following: 1. The GAMIT software used for this project should be updated regularly since some of the le formats and tables become obsolete in no time, making it incompatible with new version of the software 2. Knowledge in C programming and LINUX will be a priceless asset since all the codes that come with GAMIT are in C language and runs only under LINUX system. Editing and customizing for ones use requires some knowledge in C programming language. 3. For reliable accuracy of nal results to investigate tectonic displacements, time span of data to be processed should be at least 10 years. 4. Data to be processed should be chosen such that the campaigns were carried out in the same season in order to minimize the unknown impact of seasonality on GPS results, thus reducing aliasing problems in the incipient time series.

61

Bibliography
[Bru93] W.M Brunner, F.K. & Welsch. Effects of The Troposhere on GPS Measurements. GPS World, 4(1):4251, 1993. [DCWM05] A. Demoulin, J. Campbell, A. De. Wulf, and A. Muls. GPS monitoring of vertical ground motion in northern Ardenne-Eifel:ve campaigns of the HARD project. International Journal of Earth Sciences, 94(3), 2005. [1] Danish GPS Centre (DGC). About (dgc). Website, Aalborg University, Available at http://gps.aau.dk/about.htm, August 2003 (accessed September, 2006). [Eng01] Pratap Misra & Per Enge. Global Positioning System. Signals, Measurements, and Performance. Ganga-Jamuna Press, 2001. [Fre96] Gregory T. French. An introduction to GPS; What it is and how it works. GeoResearch, Inc., 1996. [HKM06] T. A. Herring, R. W. King, and S. C. McClusky. Documentations for the gamit/globk gps analysis software. Introduction to gamit/globk, Department of Earth, Atmospheric, and Planetary Sciences Massachusetts Institute of Technology, (accessed November 6, 2006). [HLC01] B. Hoffman-Wellenhof, H. Lichtenegger, and J. Collins. GPS Theory and Practice. Springer-Verlag Wien New York, 2001. [Klo91] J.A. Klobucha. Ionospheric Effects on GPS. GPS World, 2(4):4851, 1991. [Lei04] Alfred Leick. GPS Satellite Surveying. John Wiley & Sons, Inc., third edition, 2004. [LJSD06] M. Lidberg, J. M. Johansson, Hans-Georg Scherneck, and J. L. Davis. An improved and extended gps-derived 3d velocity eld of the glacial isostatic adjustment (gia) in fennoscandia. GPS Solutions, 81(3):213230, 2006. [pbtNWUC06] Report prepared by the NIMA WGS 84 Update Committee. Department of defense world geodetic system 1984, its denition and relationships with lokal geodetic systems. Technical Report Nima TR8350.2, National Imagery and Mapping Agency, Available at: http://earthinfo.nga.mil/GandG/publications/tr8350.2/wgs84n.pdf, 4 July 1997 (accessed October 31, 2006). [Riz99] Chris Rizos. Basic GPS Notes. Technical report, University of New South Wales, Australia, 1999. 62

BIBLIOGRAPHY [SB97] Gilbert Strang and Kai Borre. Linear algebra, Geodesy and GPS. WellesleyCambridge Press, 1997. [Sic01] John Van Sickle. GPS for Land Surveyors. Taylor & Francis Group, second edition, 2001.

63

List
2.1 6.1 7.1 7.2 7.3 7.4 7.5

of

Tables
12 50 54 54 54 54 55 68

Overall error budget . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Summary of velocities of site AAUC . . . . . . . . . . . . . . . . . . . . . . . Baselines and their sigmas and rates in 2000 Baselines and their sigmas and rates 2001 . Baselines and their sigmas and rates in 2004 Baselines and their sigmas and rates in 2005 Baselines and their sigmas and rates in 2006 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

A.1 GAMIT Scripts with their inputs and outputs . . . . . . . . . . . . . . . . . .

64

List
6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 6.10 6.11 6.12 6.13 6.14 6.15 7.1 7.2 7.3 7.4 7.5 7.6 7.7 8.1 8.2 8.3 8.4

of
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figures
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 44 45 45 46 46 47 47 48 49 49 50 50 51 51 52 53 54 55 55 56 56 58 59 59 60

5-day variation of the north-coordinate in year 2000. 5-day variation of the east-coordinate in year 2000. . 5-day variation of the up-coordinate in year 2000. . . 5-day variation of the north-coordinate in year 2001. 5-day variation of the east-coordinate in year 2001. . 5-day variation of the up-coordinate in year 2001. . . 5-day variation of the north-coordinate in year 2004. 5-day variation of the east-coordinate in year 2004. . 5-day variation of the up-coordinate in year 2004. . . 5-day variation of the north-coordinate in year 2005. 5-day variation of the east-coordinate in year 2005. . 5-day variation of the up-coordinate in year 2005. . . 5-day variation of the north-coordinate in year 2006. 5-day variation of the east-coordinate in year 2006. . 5-day variation of the up-coordinate in year 2006. . . Map of Fennoscandia. . . . . . . . . . . . . Map showing the stabilization stations. . . . 5-day variation of the baseline in year 2000. 5-day variation of the baseline in year 2001. 5-day variation of the baseline in year 2004. 5-day variation of the baseline in year 2005. 5-day variation of the baseline in year 2006. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Yearly variation of the north coordinate. . . . . . . Yearly variation of the east coordinate. . . . . . . . Yearly variation of the up coordinate. . . . . . . . . Map displaying the direction of movement of sites.

65

Appendix
These are the main le formats of GAMIT and their description A - le: ASCII version of the T-le (tabular ephemeris) B - le: controls the batch mode of data processing C - le: observed - computed (O-Cs), partial derivatives D - le: driver le of sessions and receivers E - le: broadcast ephemeris, in RINEX navigation le or FICA Blk 9 format G - le: orbital initial conditions and non-gravitational parameter values H - le: adjustments and full variance-covariance matrix for input to GLOBK I - le: receiver clock polynomial input J - le: satellite clock polynomial coefcients K - le: values of receiver clock offset during observation span, from pseudorange L - le: station coordinates M - le: controls merging of data (C-) les for solve and editing programs N - le: data-weight overrides for solve created from autcln.sum.postt O - le: record of the analysis (reduced form of Q-le) for post-processing analysis P - le: record of a model run Q - le: record of the analysis (solve run) S - le: no longer used T - le: tabular ephemeris U - le: loading and meterological data for model V - le: editing output of SCANRMS W - le: meteorological data in RINEX met-le format 66

GAMIT File Formats

A.1. SUMMARY OF GAMIT PROCESSING X - le: input observations Y - le: satellite yaw parameters Z - le: output meteorological data

A.1 Summary of GAMIT Processing

1. Create an experiment directory with subdirectories for tables and each session (day) 2. Run links.tables in the tables subdirectory 3. Assemble RINEX (or FICA) les from diskettes, tape, and/or ftp 4. Prepare an L-le 5. Enter antenna heights and receiver software versions in station.info 6. Enter scenarios in session.info (optional if running MAKEXP) 7. Edit the template sestbl. and sittbl. to set the controls for the analysis 8. For each day: Run links.day Run MAKEXP to get a D-le and batch le for MAKEX 9. Run MAKEJ to get a J-le 10. Run MAKEX with [day].makex.batch to get X- and K-les 11. Create a G-le and T-le from IGS or broadcast information 12. Run FIXDRV with the D-le as input 13. Run b[expt].bat which executes ARC, MODEL, AUTCLN, CFMRG, and SOLVE to obtain a solution for station and orbital parameters, written to q[expt]a.day, o[expt]a.day, and h[expt]a.day ), and updates g[expt]b.[day] , l[expt]a.[day] , and m[expt]a.[day] 14. Check q[expt]a.[day] for good nrms (<0.3), small station and orbital parameter adjustments, and good ambiguity resolution; and autcln.post.sum for number of observations and rms of residuals for each station. 15. Delete the C-les to save disk space; compress and/or backup RINEX and X-les

A.2 GAMIT Input and Output Files

Here, we present the various programs run in GAMIT with their input and output les whiles processing. 67

APPENDIX A. GAMIT FILE FORMATS

Program makexp

makej makex

arc xdrv

model

autcln cfmrg solve

cview scandd

Table A.1: GAMIT Scripts with their inputs and outputs INPUT OUTPUT - RINEX (or X-)les -D-le - station.info - session.info (optional) - session.info - RINEX nav le - J-le (satellite clock le) - C-le (optional) - raw observations (RINEX or FICA) - K-le (receiver clock) - station.info (rcvr, ant, rware, HI) - X-le (input observations) - session.info (scenario le) - RINEX nav le - J-le (satellite clock le) - L-le (coordinates of stations) - arc.bat (batch input le) - arcout.ddd (output print le) - G-le (orbital initial conditions) - T-le (tabular ephemeris for all sat.ses.) - D-le (list of X-, J-, L-, T-les) - B-le (bexpy.bat : primary batch le) - sestbl. (session control) - B-le (bexpy.nnn : secondary batch les) - .sittbl. (site control) - I-le (rcvr clock polynomials) - T, J, L, X (or C) input - L-le (site coordinates) - C-le ( residuals and partials ) - station.info (ant heights) - P-le (documentation of models) - X-le - I, J, T-les - antmod.dat (PCV models) - RINEX met le - otl.list/grid, atml.list/grid -C-le -C-le (cleaned) -C-le -M-le (points to the C-les) - C-le - Q-le - M-le - G-le - H-le - L-le - M-le and C-les - C-les (cview only)

68

Appendix

GLOBK Processing Summary, File Formats, and Examples

In this chapter, we write briey about GLOBK processing summary and some examples of input and output les

B.1 Summary of globk analysis of GPS data

Convert the experiment ASCII h-les into binary h-les (readable by globk) using the program htoglb. Run glred/glorg for all the (binary) h-les from a survey (or period of continuous observations) to obtain a time series of station coordinates, which can then be plotted and examined for outliers. Remove outliers from the corresponding h-les by renaming the outlier station in an earthquake le. Run globk to combine daily h-les into a single h-le that represents the averaged station coordinates for the chosen time range (e.g., monthly averages). Run glred/glorg and globk/glorg again (using the combined h-les) to obtain a time series from glred/glorg, and estimates of station velocities from globk/glorg for the entire period spanned by your data

B.2 Example of globk/glorg output

The le below shows part of the output of globk obtained by processing data of the rst 30 days of 2006 with loosely contraints.
--------------------------------------------------------GLOBK Ver 5.11S, Global solution --------------------------------------------------------Solution commenced with: 2006/ 1/ 1 0: 0 (2006.0000) Solution ended with : 2006/ 1/30 23:59 (2006.0822) Solution refers to : 2006/ 1/30 11:59 (2006.0808) [Seconds tag 45.000] Satellite IC epoch : 2006/ 1/30 12: 0 0.00 GPS System Information : Time GPST Frame J2000 Precession IAU76 Radiation model BERNE MODELS Used in Analysis: SD-WOB SD-UT1 RAY-MOD E-Tide K1-Tide PoleTideOC-Load MeanPTD Run time : 2007/ 7/16 22:34 21.00 SUMMARY VELOCITY ESTIMATES FROM GLOBK Ver 5.11S Long. Lat. E & N Rate E & N Adj. E & N +RHO H Rate H adj. (deg) (deg) (mm/yr) (mm/yr) (mm/yr) (mm/yr)

+-

SITE

69

APPENDIX B. GLOBK PROCESSING SUMMARY, FILE FORMATS, AND EXAMPLES

13.066 11.926 9.987 52.380 57.396 57.014 20.50 16.88 -0.34 -3.54 15.16 16.48 1.18 -0.37 -0.34 -17.78 1.57 16.48 47.67 37.27 35.66 37.39 0.107 34.14 -0.010 35.27 0.086 -1.26 1.26 3.15 0.04 -1.28 3.15 47.05 POTS 46.53 ONSA 46.57 AAUC

PARAMETER ESTIMATES FROM GLOBK Vers 5.11S # PARAMETER Estimate Adjustment Sigma Int. POTS_GPS 3800689.76800 882077.25900 5028791.24600 -0.01610 0.01610 0.00770 1997.002 Loc. POTS_GPS N coordinate (m) 5830872.10311 35.32594 0.07886 Loc. POTS_GPS E coordinate (m) 887872.74079 -0.33742 0.09292 Loc. POTS_GPS U coordinate (m) 129.87850 -14.64732 0.26423 NE,NU,EU position correlations -0.0922 0.1633 -0.0422 Loc. POTS_GPS N rate (m/yr) -0.00354 -0.01778 0.03739 Loc. POTS_GPS E rate (m/yr) 0.02050 0.00118 0.04767 Loc. POTS_GPS U rate (m/yr) -0.00126 0.00004 0.04705 NE,NU,EU rate correlations ONSA_GPS N coordinate (m) ONSA_GPS E coordinate (m) ONSA_GPS U coordinate (m) NE,NU,EU position correlations Loc. ONSA_GPS N rate (m/yr) Loc. ONSA_GPS E rate (m/yr) Loc. ONSA_GPS U rate (m/yr) NE,NU,EU rate correlations Loc. Loc. Loc. 0.1070 0.0053 -0.0267 6389249.75971 34.29040 0.07945 715321.47523 -0.93611 0.09334 34.11634 -11.57684 0.26538 -0.0830 -0.1692 -0.0108 0.01516 0.00157 0.03414 0.01688 -0.00037 0.03727 0.00126 -0.00128 0.04653 -0.0103 -0.0108 0.0013

Int. AAUC_GPS 3427985.87000 603660.70000 5326788.90000 0.00000 0.00000 0.00000 1997.002 Postion of AAUC_GPS referred to 2006.0882 XYZ offsets -0.0121 0.0794 -0.0451 years Loc. AAUC_GPS N coordinate (m) 6346783.95363 34.40801 0.07881 Loc. AAUC_GPS E coordinate (m) 605274.68140 -1.31290 0.09373 Loc. AAUC_GPS U coordinate (m) 48.75524 -11.79060 0.26529 NE,NU,EU position correlations -0.0849 -0.1438 0.0492 Loc. AAUC_GPS N rate (m/yr) 0.01648 0.01648 0.03527 Loc. AAUC_GPS E rate (m/yr) -0.00034 -0.00034 0.03566 Loc. AAUC_GPS U rate (m/yr) 0.00315 0.00315 0.04657 NE,NU,EU rate correlations 0.0863 -0.0227 -0.0142 GLOBK: BASELINE LENGTHS BASELINE POTS_GPS to ONSA_GPS POTS_GPS to AAUC_GPS ONSA_GPS to AAUC_GPS

Length (m) 562965.90979 552470.91643 124603.07114

Adjust (m) -2.17114 -2.16749 -0.84202

Sigma (m) 0.00597 0.00589 0.00175

GLOBK: BASELINE COMPONENT RATES OF CHANGE Baseline Length North Est. Adj. +Est. Adj. +(mm/yr) (mm/yr) ONSA_GPS-POTS_GPS 18.68 19.32 37.17 18.71 19.35 AAUC_GPS-POTS_GPS 25.70 32.60 37.44 20.02 34.26 AAUC_GPS-ONSA_GPS 15.44 -5.33 24.46 1.32 14.91

Est. 37.63 40.33 29.76

East Adj. +(mm/yr) -3.62 -1.55 -20.84 -1.52 -17.22 0.03

Rne Est. 59.53 0.112 56.40 0.277 24.06 -.029

70