/ i -x- /D -._A ......

i " 'RF Project 761054/711055 //-/O_-_SReport No. 384 _fE/r_

l

BASELINE ESTIMATION FROM SIMULTANEOUS SATELLITE

i LASER TRACKING

George C. Dedes

I Department of Geodetic Science and Surveying

l

{WASA-CR-182750_ B&SELIHE ESTIBITIOH YROB H88-22052

i SIBgLT&HEOgS S&TELLITE IASER TE&CKIHG (Ohio

State Univ. l, 204 p CSCL 22& Onclas

i G3/13 0140235

i NATIONAL AERONAUTICS AND SPACE ADHINISTRATION

Goddard Space Flight Center

i Greenbelt, Maryland 20771

i Grant No. _SG 5-265, Supp. No. I0

i

I October 198

i

i _ rl The Ohio State University

I tt ).,_'_"_1II Research Foundation

v_ 1314 Kinnear Road

Columbus, Ohio 43212 I

I
|

I

i

I

I

I

I

i

I

I

I

I

!

I

I

I

I

I

I

BASELINE ESTIMATION FROM SIMULTANEOUS SATELLITE

LASER TRACKING

by

George C. Dedes

Report No. 384

Department of Geodetic Science and Surveying

The Ohio State University

Columbus, Ohio 43210-1247

October, 1987
IiIIIIIIIIIIII

DEDICATION

To Alexandra

...who never ceases

to wonder...

ii
l

I

I

I

l

I

l

I

I

I

I

I

I

I

I

I

I

I

I

PREFACE

This project is under the supervision of Professor Ivan I. Mueller,

Department of Geodetic Science and Surveying, The Ohio State

University. The science advisor is Dr. David E. Smith9 and the technical

officer is Dr. Gilbert B. Mead, both at Code 601, Crustal Dynamics

Projectj Space and Earth Sciences Directorate, NASA Goddard Space

Flight Center, Greenbelt, Maryland 20771. The work is carried out

under NASA Grant No. NSG 5265, OSU Research Foundation Project No.

711055.

This report is an expanded version of a dissertation submitted to

the Graduate School of The Ohio State University in partial fulfillment of

the requirements for the PhD degree.

0,o

111
I

I

I

I

I

I

I

I

I

I

I

I

I

I

i

i

I

I

ABSTRACT

Simultaneous Range Differences (SRD's) to Lageos are obtained by

dividing the observing stations into pairs with quasi-simultaneous

observations. For each of those pairs the station with the least number

of observations is identified, and at its observing epochs interpolated

ranges for the alternate station are generated. The SRD observables

are obtained by subtracting the actually observed laser ranges of the

station having the least number of observations from the interpolated

ranges of the alternate station. On the basis of these observables

semidynamic single baseline solutions have been performed. The aim of

these solutions is to further develop and implement the SRD method in

the real data environment, to assess its accuracy, its advantages and

disadvantages as related to the range dynamic mode methods, when the

baselines are the only parameters of interest.

Baselines, using simultaneous laser range observations to Lageos,

have also been estimated through the purely geometric method. These

baselines formed the standards of comparison in the accuracy

assessment of the SRD method when compared to that of the range

dynamic mode methods. On the basis of this comparison it was

concluded that for baselines of regional extent (i.e., up to 3700 km) the

SRD method is very effective, efficient and at least as accurate as the

range dynamic mode methods, and that on the basis of a simple orbital

modeling and a limited orbit adjustment.

The SRD method is insensitive to the inconsistencies affecting the

terrestrial reference frame and simultaneous adjustment of the Earth

Rotation Parameters (ERP's) is not necessary. Therefore, this method

offers an inexpensive alternative for projects designed to study regional

plate tectonic motions.

iv
I

I

I

i

I

I

I

I

I

I

I

I

I

l

i

I

I

I

ACKNOWLEDGMENTS

Although it is difficult to acknowledge all those individuals who

contributed to the successful completion of this work in the the long

journey of studies, acknowledgments are given but by no means is this

list meant to be complete.

My deepest appreciation and thanks go to my wife Alexandra for her

endless love and her continuing patience, support and encouragement

tirelessly shown during the difficult but very fruitful years of my

graduate studies. Special thanks and appreciation are due to my

parents who, with their tireless guidance and their continual support in

my high school and undergraduate years, provided me with courage,

faith and love for knowledge making the graduate studies a long and

very challenging experience. The support and understanding of my

parents-in-law during the last years of my studies is also greatly

appreciated.

My gratitude and appreciation are also due to my supervisor and

mentor Dr. I.I. Mueller for his expert guidance, constructive criticism,

continual encouragement and tireless support during the course of this

work and throughout my graduate studies at the Department of Geodetic

Science and Surveying. His patience and continual involvement to

successfully provide financial support even in very difficult times is

greatly appreciated. His well organized and vast personal library, ready

at my disposal, was at times of invaluable help and for this I am indeed

very thankful.

I wish also to express my thanks to Dr. U. A. Uotiia and Dr. C. C.

Goad for their constructive criticisms and their valuable suggestions

which enhanced the quality of this work. Special thanks are also given

to Dr. Walter Mitchell for reading the dissertation, in which form this

report was first presented.

V
The overall guidance of the faculty members provided in the

coursework throughout my studies at Ohio State is also greatly

appreciated.

I would especially like to thank former student Dr. E.C. Pavlis (now

of EG&G Analytical Services) for his valuable suggestions in regard to

the successful implementation of GEOSPP software, which he developed

during his studies at Ohio State. His continuing support in providing

information on Lageos' orbit is also greatly appreciated. The

discussions with former student Dr. Brent Archinal (now at the U.S.

Naval Observatory) helped to clarify many aspects in regard to the

efficient use of the computer facilities which proved to be of great

value in the course of this study. Acknowledged are also the valuable

suggestions of David Rowlands, a former student of this Department,

related to the efficient use of GEODYN programs.

Special thanks are extended to my fellow students past and present,

Mike Baker, Lucia Tsaoussi, Nikos Pavlis, Vasilios Despotakis, George

Priovolos, Dr. Y. Bock (now of the Massachusetts Institute of

Technology) and Dr. H.B. Iz (now of ST Systems Corporation), who,

although deeply occupied with their own work, were always ready to

listen, discuss and offer valuable advice for my problems.

Thanks are also extended to Dr. D. C. Christodoulidis (now of Jet

Propulsion Laboratory, Pasadena) for providing the data sets employed

at the initialstage of this study and also to Dr. Steve Klosko at EG&G

for furnishing valuable information about the quality, distribution and

simultaneity of the laser range observations to Lageos. Dr. Henry

Linder at GSFC, with his valuable comments, clarified many aspects

related to the contents of the tapes containing the laser range

observations to Lageos.

The financial support for this work provided by NASA Research

Grant NSG 5265, OSURF Project 711055, is also greatly appreciated. Also

acknowledged is the financial support of the American Association of

Geodetic Surveying (AAGS)-American Congress on Surveying and

Mapping(ACSM) Graduate Fellowship (1984) and the Wild Heerbrugg-ACSM

Geodetic Fellowship(1986).

vi

I

I

i

I

I

!

I

I

I

I

I

I

I

I

i

I

l

I
I

I

I

I

I

I

I

I

I

l

I

l

I

I

I

I

I

I

I

The extensive computer funds necessary for the completion of this

work were provided by the Instruction and Research Computer Center

(IRCC) of The Ohio State University.

Last but not least I wish to thank Irene Tesfai for her excellent

typing and valuable suggestions as well as Karen Wasielewski and Tracy

Runyon for the excellent typing of parts of this manuscript.

vii
I

I

I

I

I

i

I

I

I

I

I

I

I

I

I

I

I

I

TABLE OF CONTENTS

Page

Preface ........................................................... iii

Abstract .......................................................... iv

Acknowledgments ................................................... v

List of Tables .................................................... xi

List of Figures ................................................... xii

Abbreviations and Acronyms ........................................ xiii

I. INTRODUCTION .................................................. 1

i.i Baseline Estimation in the Dynamic and Semidynamic

Environment .............................................. 1

1.2 Scope of This Investigation .............................. 2

2. ESTIMATION METHODS ............................................ 6

2.1 Geometric Methods ........................................ 6

2.1.1 Mathematical Model ................................ 7

2.1.2 Normal Equations .................................. 8

2.1.3 Critical Configurations ........................... 12

2.2 Dynamic and Semidynamic Mode Methods ..................... 17

2.2.1 Simultaneous Range-Difference Semidynamic Mode

Method ............................................ 20

2.2.2 Mathematical Modeling ............................. 20

2.2.3 Orbit Determination with the Method of Special

Perturbations ..................................... 23

2.2.4 Reference Frames and Systems ...................... 26

2.2.5 Orbital Model ..................................... 30

2.2.6 Normal Equations .................................. 41

3. GENERATION OF THE OBSERVABLES ................................. 46

3.1 Satellite Laser Ranging .................................. 46

3.2 Satellite Laser Ranging System, Its Components and

Their Contribution to the Total Error Budget ............. 47

3.2.1 Hardware of the Ground Segment .................... 47

3.2.2 Atmospheric Channel ............................... 51

3.2.3 Space Segment ..................................... 52

3.2.4 Instrument Origin ................................. 53

3.3 Systematic Corrections of the Observations External to

the SLR System ........................................... 55

3.4 Description of the Data Set Utilized in This

Investigation ............................................ 57

Lx

P_CF_)_G PAGE BLANK NOT FILMED
3.5

3.6

Data Editing ............................................. 60

3.5.1 Data Snooping Procedure ........................... 61

3.5.2 Effectiveness of the Data Snooping Procedure in

Editing the Laser Range Observations .............. 67

Generation of Simultaneous Ranges and Simultaneous

Range Differences ........................................ 75

3.6.1 Chebychev Polynomials and Spline Functions in the

Context of Global and Piecewise Interpolation ..... 78

3.6.2 Chebychev Polynomials vs. Cubic Spline Functions

in the Functional Representation of Laser Ranges .. 84

3.6.3 Data Selection for the Generation of the

Simultaneous Ranges and the Simultaneous Range

Differences ....................................... 92

4. BASELINE ESTIMATION ........................................... 119

4.1 Baseline Estimability .................................... 119

4.2 Steady State Response of the Geometric and SP_ Methods ... 119

4.3 Baseline Estimation via the Steady State Response of

the Geometric Method ..................................... 123

4.3.1 Geometric Strength of the Available Observations .. 125

4.3.2 Baseline Results .................................. 130

4.4 Baseline Estimation via the Steady State Response of

the SP_ Method ........................................... 135

Pass Baseline Geometry and Its Manifestation in

the Design Matrix of the SRD Observable ........... 137

Baseline Results .................................. 141

4.5 Baseline Comparison ...................................... 161

4.6 Response of the SRD Method to the Simplifications of

the Orbital Model ........................................ 164

5. CONCLUSIONS AND RECOMMENDATIONS ............................... 171

Appendix A:

Appendix B:

Partial Derivatives of the SED Observable ............ 178

Sensitivity of the Partial Derivatives with Respect

to the Initial State Vectors ......................... 181

References ........................................................ 183

X

I

I

I

I

I

i

I

I

I

I

I

I

I

I

I

l

I
!

l

!

!

!

I

i

!

!

!

!

!

I

I

I

!

i

TABLE

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

LIST OF TABLES

Page

Lageos Along-Track Acceleration and Its Reflectivity

Coefficients ................................................ 39

Station Location, Laser Instruments and Precision Estimates . 55

Monthly MERIT Releases ...................................... 60

Monthly Precision Estimates (cm) for the American Stations

(Chebychev Fitting/After Data Snooping) ..................... 107

Precision of Chebychev Interpolation, Before and After

Data Snooping ............................................... 115

Degree of Chebychev Polynomials ............................. i17

Global Statistics of the Geometric Adjustments .......... :.... 127

Baseline Steady State Response of the Geometric Solution .... 128

Steady State Response of Baseline 7109-7886, Parallel Passes 146

Steady State Response of Baseline 7110-7220, Parallel Passes 147

Steady State Response of Baseline 7110-7265, Passes Within

*30" - *50" . ................................................ 148

Steady State Response of Baseline 7109-7110, Parallel Passes 149

Steady State Response of Baseline 7110-7886, Parallel and

Perpendicular Passes ........................................ 150

Steady State Response of Baseline 7110-7086, Passes Within

*20" - *60" . ................................................ 154

Steady State Response of Baseline 7110-7122, Long-Arc Mode,

Parallel Passes ............................................. 156

Steady State Response of Baseline 7109-7105, Long-Arc Mode,

Passes Within *20" - _30" . .................................. 157

Steady State Response of Baseline 7110-7105, Long-Arc Mode,

Passes Within *20" - *60" . .................................. 160

Baseline Differences ........................................ 163

Baseline Differences with Respect to Those Computed Using

an Orbital Model Including a 12x12 Gravity Field + (i) ...... 167

Steady State Response of the SRD Method ..................... 173

xi
FIGURE

1

2

3

4

5

6

7

8

10

11

12

13

14

15

16

17

18

LIST OF FIGURES

Residuals of Chebychev interpolation, before data snooping ..

Residuals of Chebychev interpolation, after data snooping ...

Residuals of Chebychev interpolation, after data snooping

Page

71

73

and editing residuals greater than 0.4 m ................... 76

Orbit residuals (cubic spline SRD's) ........................ 88

Successive data gaps for station 7115 ....................... 89

Orbit residuals (Chebychev SRD's) ........................... 90

Differences between cubic spline and Chebychev SRD's ........ 91

Arc overlap of two Lageos passes recorded by American

stations and European stations .............................. 95

Recovery errors and distribution of ground truth observations

(dense data) ................................................ 97

Recovery errors and distribution of ground truth observations

(sparse data) ............................................... 99

Chebychev residuals ......................................... i00

Location of the American stations used in the present study 104

Arc overlap of three Lageos passes (American stations) ...... 105

Chebychev residuals for stations 7907, 7105, 7086, 7109,

7110, 7886, after data snooping ............................. 108

Chebychev residuals for stations 7112, 7265, 7121, 7210,

7220, 7122, 7062, after data snooping ....................... IIi

Arc overlap of a Lageos pass involving four American

stations .................................................... 118

Subsatellite track-baseline geometry ........................ 138

Station location and Lageos groundtracks .................... 145

xii
IIIIIIIIIII

i

IIIIII

AAGS

ACSM

AU

BIH

BTS

CCR

CCRS

CDDB

CDIS

CEP

COTES

CRF

CSR

CTRS

EG&G

ERP

FR

GAST

GEML2

GEOM

GMST

GPS

GSFC

GTLN

IAU

IERS

IUGG

LAGEOS

bEP

MERIT

M_C

NASA

OTIU

PC

PLN

P_ff

PPN

RMS

SLR

SR

SRD

ABBREVIATIONS AND ACRONYMS

American Association of Geodetic Surveying

American Congress on Surveying and Mapping

Astronomical Unit

Bureau International de l'Heure

BIH Terrestrial Reference System

Cube Corner Reflectors

Conventional Celestial Reference System

Crustal Dynamics Data Bank

Crustal Dynamics Information System

Celestial Ephemeris Pole

IAU/IAG Joint Working Group on the Establishment

and Maintenance ofa Conventional Terrestrial

Reference System

Celestial Reference Frame

Center for Space Research

Conventional Terrestrial Reference System

EG&GWashington Analytical Science Center

Earth Rotation Parameters

Full Rate

Greenwich Apparent Sidereal Time

Goddard Earth Model L2

GEOMetric solution

Greenwich Mean Sidereal Time

Global Positioning System

Goddard Space Flight Center

Goddard Laser Tracking Network

International Astronomical Union

International Earth Rotation Service

International Union of Geodesy and Geophysics

LAser GEOdynamic Satellite

MicroChannel Plate

Monitor Earth Rotation and Intercompare the Techniques

of observation and analysis

Main MERIT Campaign

National Aeronautics and Space Administration

Optical Time Interval Unit

Personal Computer

Participating Laser Network

PhotoMultiplier Tube

Isotropic Parametrized Post-Newtonian n-body metric

Root Mean Square

Satellite Laser Ranging

Simultaneous Ranges

Simultaneous Range Differences

xiii
SSD

TAI

TDB

TDT

TIU

TLRS

TRF

USNO

UT

UTC

UTX

UTI

UTIR

Z IPE

SeaSat Decimal format

International Atomic Time

Barycentric Dynamical Time

Terrestrial Dynamical Time

Time Interval Unit

Transportable Laser Ranging System

Terrestrial Reference Frame

United States Naval Observatory

Universal Time

Universal Coordinated Time

University of Texas

Universal Time corrected for polar motion

UTI with the effects of zonal earth tides up to 35

days removed

Central Institute for Physics of the Earth at Potsdam

xiv

I

I

I

I

I

I

I

I

I

I

I

1

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

l

I

I

Chapter 1

INTRODUCTION

1.1 BASELINE ESTIMATION IN THE DYNAu_IC AND SEMIDYNA_C

ENVIRONMENT

In the dynamic environment accurate baseline estimation requires a

highly sophisticated orbital modeling and a baseline-pass geometry

leading to near cancellation of the accumulated along-track and

cross-track orbital errors caused by the erroneous constraints imposed

on a large number of estimable quantities (Rao, 1973), the recovery of

which is not possible due to their reduced data sensitivity. In this

environment proper implementation of the Terrestrial Reference Frame

(TRF) requires simultaneous recovery of the Earth Rotation Parameters

(ERP) or utilization of a consistent set of ERPs obtained through a

separate step.

Although fulfillment of these requirements makes it possible to

effectively recover baselines of global and regional extent, it results in

low temporal resolution of baseline recovery.

In the semidynamic environment (Section 2.2), and on the basis of

simultaneous observations, only regional baselines can be recovered with

an accuracy compatible to that of the observed laser ranges. The

maximum regional baseline length effectively recovered in this

environment depends on whether the simultaneous observations collected

by the baseline end stations are enough to result in a steady state

response (Section 4.2.1). This, however, is a function of the satellite

altitude, and for the Lageos satellite the effective regional extent may

include baselines of up to 3703 km (Section 4.4.2).

In the semidynamic environment a relatively simplified orbital model

is required and only the position and the orientation of each of the arcs

involved is adjusted to "best" fit the available observations (Section
2.2.6). Adjustment for the ERP parameters is not necessary since

proper implementation of the TRF frame is warranted by the use of only

simultaneous observations. The relaxing data requirements and the

limited orbit adjustment make it possible to increase the resolution of

baseline recovery without any loss on the achieved accuracies, and at

the same time to substantially decrease the required computations,

thereby making it possible to effectively implement the semidynamic

methods with limited computer facilities as, for instance, in the personal

computer (PC) environment.

The sophistication of the orbital modeling in the semidynamic

environment can be further simplified by appropriately transforming the

observed laser ranges to bring them "closer" to the estimable quantities

being recovered (i.e., baselines). The term "closer" indicates that on

the basis of the same orbital model the errors affecting the computed

value of the transformed observations, referred to from now on as

"observables," are smaller than those affecting the computed value of

the observations themselves. By bringing the observations closer to the

estimated baselines, the sophistication of the orbital modeling could be

further reduced if the performed transformation cancels out the errors

caused by the model simplifications.

Transformation of the laser range observations to Simultaneous

Range Difference (SRD) observables brings them closer to the estimated

baselines (Pavlis, 1982). The potential of using SRD observables to

estimate baselines was studied with simulated data by Pavlis (1982). The

results of this study were very promising not only with respect to their

accuracy but also with respect to their simplicity.

1.2 SCOPE OF THIS INVESTIGATION

For the reasons mentioned in the previous section, it was considered

appropriate and worthwhile to pursue the present investigation, the aim

of which is to further develop and implement the SRD method (Section

2.2.1) in the real data environment, to assess its accuracy, its

advantages and disadvantages as related to range dynamic mode

methods, when the baselines are the only parameters of interest.

2

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

i
Since during the MERIT Main Campaign many stations collected

simultaneous observations (Section 3.4), it was decided to proceed not

only with the development of the SRD method but also with baseline

estimation through the geometric method (Section 4.3). These baselines

formed the standards of comparison in the accuracy assessment of both

the SRD and the range dynamic mode methods (Section 4.5).

In pursuing this study, the geometry does not take part in the

physical events, as happens in general relativity, but rather it is used

in its deductive form either purely (i.e., geometric methods) or in

combination with an inductive form (i.e., dynamic and semidynamic

methods), both of which were employed to formulate the spatial-temporal

relationships of the observing stations and the observed satellite

positions.

In a pure deductive mode the geometric method is entirely based on

Euclidean geometry, without any reference to the inductive inference

that the satellite moves along the path chosen by its physical

environment. The implied Euclidean geometry is revealed, in the

arithmetic framework, on the basis of Cartesian coordinates (Section

2.1.I). There exist, however, configurations in which some of the

vertices of the resulting figures are free to move along a locus of

points, thereby forming configurations that are not unique and are

referred to as "critical configurations" (Section 2.1.3).

In the inductive mode, the differential form of a satellite's motion,

modeled on the basis of Newton's inductive laws, is expressed through

its equations of motion (Section 2.2.5). These equations are integrated,

on the basis of the deductive mathematical framework, to reveal the

geometric path of the satellite (i.e., orbit), chosen by the satellite's

physical environment (Section 2.2.3). Having the geometric path of the

satellite, the observing stations are connected with this path through

the Euclidean distance formulated in terms of Cartesian coordinates.

Since the simplest geometry of the satellite orbit reveals itself in an

inertial reference frame, the satellite'sequations of motion are expressed

with respect to such a frame (Section 2.2.4). However, the estimated

Cartesian coordinates of the observing stations are referenced with
respect to an earth-fixed frame (Section 2.2.2}.

The deductive and/or inductive formulation described above contains

a large number of slow varying quantities which can be considered

constant for the time span of the observations and a subset of which

constitutes a set of quantities that are estimable only if the necessary

units, constants and/or constraints have been properly adopted for

their unique determination (Rao, 1973).

There are three types of estimable quantities that can be

differentiated from their relation to the physical environment, or the

observing environment, or the links of these two environments.

In satellite geodesy the interstation distances (i.e., baselines) are

estimable quantities related to the observing environment. The estimable

quantities of the physical environment are those related to its cause

(i.e., Anm and Bnm potential coefficients) and those to its effect (i.e.,

geometric and dynamic characteristics of the orbit). The estimable

quantities molded in the link of the physical and the observing

environments are those resulting from the latter as relates to the cause

and effect duality of the former (i.e., station geocentric distances,

latitudes, and longitude differences}.

The baselines are computed from the earth-fixed station coordinates

which are recovered through an inversion process, such as

least-squares adjustment, on the basis of both the geometric and the

SRD methods (Sections 2.1.2 and 2.2.6). The input to this inversion

process are the Simultaneous Ranges (SR) and SRD observables. These

observables were generated through an interpolation of the observed

laser ranges because it is quite unlikely, if not impossible, to record

simultaneous observations to a passive satellite (Section 3.6). Because of

the peculiarities of the SLR system (Section 3.1 I, it is very likely that

the recorded laser ranges will be affected not only by white noise but

also by large blunders (Section 3.21. For this reason and since the SR

and SRD observables will be generated via an interpolation, it is

important to edit the laser ranges before proceeding with the geometric

and the orbit adjustments.

4

I

1

I

I

I

I

!

!

i

I

I

I

I

I

I
!

!

I

I

I

!

I

I

I

I

I

I

I

I

I

I

!

I

I

In the geometric method the recovered baselines will only be

affected by the errors resulting from an improper or from the not yet

reached steady state response (Sections 4.2). The steady state response'

of the geometric methodj on the basis of a minimum least-squares

solution_ is affected only by the observational errors. Such a response9

however_ was not possible for longer baselines (Section 4.3.1)_ and

therefore an overconstrained solution was adopted to form the standards

of comparison with the anticipation that it is the least erroneous when it

is compared either to SRD or to dynamic solutions (Section 4.3.2).

The accuracy of the baselines_ recovered via the SRD (Section 4.4)

is assessed from the comparison with the baselines obtained via the

geometric methods and the range dynamic mode methods (Section 4.5).

The response of the SRD method to the simplification of the orbital

model has also been investigated (Section 4,6).

5
Chapter 2

ESTIMATIONMETHODS

In this chapter an attempt is made to briefly describe the mathematical

models and the principles involved in the geometric and dynamic mode

methods as they are applied to satellite geodesy. Although the

geometric solutions performed in this study are only used as standards

of comparison (Sections 4.3 and 4.5), their mathematical formulation is

presented first because historically the geometric methods were the first

ones to result in accurate differential positioning.

2.1 GEOMETRIC METHODS

Geometric methods are based on the analysis of the relative geometry of

the observations without any reference to the physical processes

creating the problem under question. On the contrary, some or all of

the systematic corrections applied to the observations are computed with

the use of physical models.

In the geometric approach of satellite geodesy (Veis, 1960; Mueller,

1964a), the observed satellite positions are treated as auxiliary

independent points in space, and they are only used to relate the

observations geometrically. This in turn leads to the generation of

space networks. These networks manifest not only the relative geometry

of the observations, but also any a priori information which is necessary

for their realization. Thus each observation relates the position of the

observing station with the observed satellite position. The unknown

parameter vector includes the Cartesian coordinates of the observing

stations together with the Cartesian coordinates of the observed satellite

positions at each of the observing epochs. Since the coordinates of the

unknown satellite positions constitute an independent set of unknowns,

6

!

II

!

!

!

I

I

!

!

I

I

!

!

!

!

!

!

!
I

I

I

1

I

I

I

I

I

I

I

I

I

I

I

I

!

I

it is necessary to have a sufficient number of observations at each

observing epoch. The number of observations should be sufficient not

only to eliminate the unknown satellite coordinates at each of those

epochs, but also to solve for the unknown station coordinates.

The process described, however, necessitates the usage of

simultaneous (referenced to satellite time) observations without any

reference to the fact that the motion of the satellite obeys the physical

laws of dynamics. These two distinct features create the advantages

and disadvantages of the geometric approach to satellite geodesy.

2.1.I Mathematical Model

The mapping of the parameter space into the observational space is

referred to as observational modeling. The analytical expression

responsible for the realization of this mapping is referred to as either

the mathematical or the observational model. The mathematical model

employed in the geometric approach is that of the Euclidean range from

a ground station to an observed satellite position expressed in terms of

Cartesian coordinates:

Fij = [(uj-ui) 2 + (vj-vl) 2 + (wj-wl)2] _ - rij = 0 (2-1)

The quantity rij is the true value of the range observable from the

ground station i to the satellite position j (see Section 3.1), while the

quantities ui, vi, wi and uj_ vj_ wj denote the true values of the

Cartesian coordinates of station i and satellite position j. These

coordinates may be referred to any arbitrarily chosen Cartesian

coordinate system since ranges are invariant under any rigid body

rotations.

The linearized form of equation (2-1) forms the basis for the

generation of the observation equations when four or more stations are

observing simultaneously (see Section 2.1.3). With these observation

equations the normal equations are derived on the basis of a weighted

constrained least squares adjustment (Uotila, 1967; Krakiwsky et al.,

1967). The resulting normal equations are reduced by eliminating the

unknown satellite coordinates. The reduced normals are then solved to

7
estimate the stations' Cartesian coordinates which are finally transformed

to interstation distances.

Inverting the normals on the basis of the minimum required

information (i.e, minimum constraints) leads to baseline errors that

depend solely on the errors of the observed ranges and on their

geometric strength as well. This is true only when the scale has been

properly incorporated into the solution either through the observations

and their geometry, or if this is not possible, through some other

additional constraints (see Section 4.3).

2.1.2 Normal Equations

The observation equations used to derive the normal equations are

obtained from the linearization of equation (2-1). Linearization is

achieved with a Taylor series expansion about the approximate values of

the station coordinates, the satellite coordinates and the observed

ranges as well. The resulting linear equations have the form:

A

Aij Xij - Vij + L|j = 0

where

(2-2)

Aij = [ alj _ -aij]

aFij

aiJ = a(uj, vj, wj)

(2-3)

(2-4)

Xij =

ou,] (IV4

_W_

(IU

dv i

. dwl

(2-s)

Vlj = residual vector corresponding to the range observable

Lij = the computed minus the observed range

8

!

I

I

!

!

I

!

!

I

I

I

I

!

!

B

I

!

D
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

In cases when the geometric strength of the observations is not

good enough to warrant a steady state response, it could still be

possible to reach such a response if estimated or observed interstation

distances are incorporated into the solution (see Section 4.3). This is

accomplished by introducing the interstation distances as fictitious

observations into the adjustment. Appropriate weights should be

applied to these (fictitious) observations in order to avoid any scale

conflicts that might contaminate the solution (Uotila, 1967).

The mathematical model of the interstation distance between stations

k and _ has the following form

Gk_ = [(U_--Uk)2 + (V_--Vk)2 + (W_--Wk)2]_ -- Lk_ = 0 (2-6)

where Lk_ is the true value of the fictitiously observed range between

these two stations. Linearization of equation (2-6) about the

approximate station coordinates and the fictitiously observed distance

results in

(2-7)

A

Ck_ Xk_ - Vk_ + Dk_ = 0

where

Ck_ = [ T_ _ -T_ ] (2-8)

aGk_

T_ = a(u_, v_, w_) (2-9)

Vk_ = residual of the (fictitiously) observed interstation distance

Dk_ = the computed minus (fictitiously) observed interstation

distance

So far the observation equations have been developed on the basis

of one range observation and one interstation distance. Considering

many range observations and many fictitiously observed interstation

distances, equations (2-2) and (2-7) take the following form:

9
(2-1o)

(2-11)

These equations can be rewritten as

[A] 0 liDIvire (2-12)

and with some obvious substitutions they take the following form:

A*X + L* = V* (2-13)

Equation (2-13) forms a set of observation equations which are used to

derive the normal equations on either deterministic or statistical

grounds.

Deterministically the principle of least squares requires that the

quantity (VZTPZV * + _TPx_ ) assumes minimum value, subject to the

condition AZX - V* + L* = 0. The matrix Px is the weight matrix

associated with the coordinates of the ground stations and of the

observed satellite positions, while the weight matrix P* takes the

following form

(2-14)

where P and Pc are the weight matrices associated with the range

observables and the (fictitiously) observed interstation distances

respectively.

Statistically the maximum likelihood principle requires maximization of

the a posteriori conditional density function of the parameter vector

given that the observations L* have been obtained. The resulting

estimator X is referred to as Bayesian least squares estimator. Both of

the above principles lead to the same estimator X only if normality is

assumed not only for the a priori density function of the estimated

I0

!

!

!

II

!

!

!

Ii

!

!

!

!

!

II

!

!

II

!
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

parameter vector X, but also for the conditional density function of the

observational vector L* given that the parameter vector ){ has been

estimated. These assumptions should only hold for the maximum

likelihood principle since any least squares estimator is a distribution

free estimator. The estimator X is obtained from the resulting normal

equations (Uotila, 1967; Krakiwsky et al., 1967; Cappellari et al., 1976):

] + [A]= 0C -P_c x k c D

(2-1s)

where k c are the Lagrange multipliers associated with the (fictitiously)

observed interstation distances. Elimination of the Lagrange multipliers

from equation (2-15) leads to the following form for the estimator X:

= -(ATpA + CTPc C + Px) -x (ATpL + CTPcD ) (2-16)

Substitution of equations (2-2) through (2-5) and (2-7) through

(2-9) into equation (2-16) followed by elimination of all the satellite

coordinates leads to the following equations (Krakiwsky et al., 1967):

- The 3x3 diagonal matrix associated with the kth ground station

Nkk = [j_ akj T Pkj akj}- j_ {akj T Pkj akj [_ at j

akj T Pkj akj} + Pk + [ 0 ] TkT Pk_ Tk 1

T Pij alj} -z

(2-17)

- The 3x3 off-diagonal matrix corresponding to the k th and _th ground

station

J

0 ] TkT Pk_ T_ 1(2-18)

- The 3xl constant vector associated with the k th ground station

11
T Pij aij)-*

i 1

In equations (2-17) and (2-19) the j summation is performed over all

the satellite positions observed from station k, while in equation (2-18)

the j summation is performed over all the satellite positions observed

simultaneously from both stations k and _ respectively. The summation

i, on the other hand, is performed over all the stations observing the

satellite position j simultaneously. In these three equations the number

1 is used only when the interstation distance between the stations k

and _ is involved in the solution, otherwise the value of 0 is used.

Equations (2-17) through (2-19) form the basis for the estimation of

the ground station coordinates by inverting the normals through a

procedure accredited to Banachiewich. This procedure is carried out in

two steps.

- The first step involves the representation of the normal matrix as

a product of right and left triangular matrices with the left

triangular matrix having unit diagonal elements.

-The second step involves the computation of the inverse normal

matrix on the basis of only the inverted diagonal elements of the

right triangular matrix and the off-diagonal elements of the left

triangular matrix.

The above procedure is very closely related to the Cholesky algorithm

(Uotila, 1967).

2.1.3 Critical Configurations

In the geometric approach of satellite geodesy the observed satellite

positions (targets) are treated as auxiliary independent points in space.

They are used only to relate the observations geometrically through the

resulting range space networks (see Section 2.1). In certain cases the

ground stations and/or the targets which form the vertices of this

12

II

I

II

II

II

I

I

II

I

I

II

i

I

II

II

II

II

II
I

II

It

I

li

II

I

I

II

I

!

I

II

It

It

I

II

I

network are involved in a kind of configuration for which a unique

adjustment is impossible although the number of observations is

sufficient and the coordinate system well defined. These configurations

are referred to as "critical."

With range observations the critical configurations have been

extensively studied in the past (Rinner, 1966; Blaha, 1971; Tsimis, 1972,

1973). The critical configurations have been traditionally analyzed

according to whether all of the observing stations are either in a plane

or generally distributed in space. For both of these cases the resulting

singularities are divided into three categories:

1. Singularity A, resulting from the relative geometry of an individual

station connected to its observed targets.

2. Singularity B, resulting from the relative geometry of the observing

stations only.

3. Singularity C, resulting from the relative geometry of all the

observing stations connected to their observed targets when

singularity A and singularity B are not present.

2.1.3.1 Critical configurations when all of the observing stations

are in a plane. When all of the observing stations are in a plane the

singularity problem is analyzed according to the number of stations

observing all the targets. This number may be three or more, or less

than three. The number three is important since ranges from three

stations are needed to eliminate the coordinates of one target, provided

that this target is not located on the plane of these three stations

(Blaha, 1971).

If the number of ground stations observing all the targets is three

or more, singularity A occurs when an individual station---excluding the

stations used to define the coordinate system, since for these three

stations singularity A cannot occur--is either observing less than three

distinct targets or is in the same plane with all of its observed targets.

Furthermore, singularity B occurs when all the observing stations but

one are lying in a straight line or more generally when all the stations

13
are lying on a second-order curve. Since at least five stations are

needed to determine a second-order curve, singularity B can only be

avoided if six or more stations are involved. In the absence of

singularities A and B singularity C occurs when the stations making

off-plane observations (i.e., observed targets are not in the same plane)

are not themselves off-curve stations (i.e., not lying in the same

second-order curve). To avoid singularity C at least three off-curve

stations should make off-plane observations (Blaha, 1971; Mueller et al.,

1975). In the case when all the stations observe all the targets,

singularity A loses its importance because it always implies singularity

Cp since off-plane observations are necessary to avoid singularity C.

When there do not exist three stations observing all the targets,

elimination of the coordinates for all the targets using the same three

stations cannot be achieved. Thus, in the elimination process one, two

or all of these three stations will have to be replaced. This leads to

the first, second and third replacement respectively, and therefore to at

least four-station events.

We first denote with k the station used in the first replacement. In

this replacement singularity A for all the stations but k, or singularity

B for all stations would occur as though there were three stations

observing all the targets. For station k, however, singularity A occurs

if any new stations coming into play are lying either on the x-axis of

the local coordinate system (e.g., line formed by two of the three

stations used to define the local coordinate system) or in the

intersection line (denoted by _) of the plane _ (see below) with the

plane of the ground stations. The plane _, if it exists, is the plane

containing station k together with the satellite positions (denoted by Jk)

that were observed by this station (e.g., k) up to the epoch at which

the new station(s) started observing. Singularity C, in the absence of

singularities A and B, is further analyzed according to whether Jk are

off-plane or in-plane targets. If Jk are off-plane targets, singularity C

occurs whenever stations making off-plane observations are not

themselves off-curve stations. The case when Jk form in-plane targets

14

II

II

II

II

II

I1

I

I

I

I

I

i

II

II

II

I

II

II
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

is not discussed here because it is very unlikely to encounter in

practice. This case, of rather academic interest, is discussed in (Blaha,

1971, page 63).

Next we denote by s" the station used in the second replacement.

In this replacement, if there is no other station s'" to start observing

for the first time after s" has started, singularity A occurs as though

only the first replacement would have taken place. However, with

station s'" present, singularity A for station k occurs if, in addition to

the above, the station s'" is lying either in the x-axis of the local

coordinate system or in the line _. Furthermore, singularity A for

station s" occurs if in addition to being in the plane _" {defined below),

the station s'" is also lying either in the line defined by the station

used as the origin of the local coordinate system and the station k or in

the intersection line {denoted by _') of the plane _" with the plane of

the ground stations. The plane _', if it exists, is the plane containing

the station s" together with the targets {denoted by is') observed by

the station s" up to the epoch at which the station s'" started

observing. Singularity B occurs as though three stations observing all

the targets exist. If is" and Jk are off-plane targets, in the absence of

singularity A and B, singularity C occurs, when no other stations

besides k and s" exist to make off-plane observations. The case when

is" or Jk form in-plane targets is not discussed here because of the

unlikelihood of encountering it in practice.

The singularities resulting from three or more replacements are

similar to the ones described above. By avoiding singularities A, B and

C a nonsingular network can only be formed if at least six stations in at

least four station events are involved. This is so because five stations

are needed to define a second-order curve and only the sixth station is

possible to serve as an off-curve station.

Once a nonsingular network has been realized any extension of it

will result in a nonsingular one if for any additional station singularity

A is eliminated and if no target is lying in the plane of the ground

stations.

15
2.1.3.2 Critical configurations when the observing stations are

generally distributed in space. When the ground stations are in general

configuration, singularity B loses its meaning because the effect of

ground stations cannot be separated from that of the targets.

Consequently, singularity B will not be considered here. However,

another type of singularity called "reverse singularity B" is the

singularity B if one assumes that the satellite points (targets) observe

the ground stations. Therefore, this singularity occurs when all the

targets are in a plane in a second-order curve. This in practice could

approximately occur when two short passes of about the same altitude

have been observed.

Having the observing stations in a general configuration, a

nonsingular network can be formed if at least six targets are being

co-observed by at least four stations. Accordingly, the analysis of

critical configurations proceeds by grouping the ground stations in

tetrads ("quads"). With four stations observing all the targets,

singularity A occurs only with respect to the fourth station because

singularity A never occurs for the three stations that have been used

to define the local coordinate system (Blaha, 1971). With respect to the

fourth station, singularity A occurs if all the targets are lying on a

plane through this station, or if all the targets are on the plane formed

by the stations used to define the local coordinate system.

Furthermore, in the absence of singularity A, singularity C occurs when

all the observing stations and all the targets are lying on a

second-order surface.

With more than four stations observing, the singularity problem is

analyzed by grouping the observing stations in quads. If the number

of stations observing all the targets is three or more, singularity A

occurs, as though all the observing stations were lying on a plane (see

Section 2.1.3.1). In the absence of singularity A, singularity C occurs

either when all the observing stations with their observed targets are

located on a second-order surface or when each tetrad of stations

together with its observed targets are located on specific second-order

critical surfaces. These surfaces intersect each other in one

16

I

!

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

second-order curve containing the three stations used to eliminate the

coordinates of each target and to define the local coordinate sys.tem.

Furthermore, when all the stations are co-observing, these second-order

critical surfaces coincide.

When three stations observing all the targets do not exist, then the

concept of station replacement should be utilized in exactly the same

way as described in the previous section. Proceeding with this concept

it is found that singularity A occurs as though three stations observing

all the targets exist. As for singularity C, it is again associated with

other specific second-order surfaces in addition to the ones resulting in

singularity C when three stations observing all the targets exist.

By avoiding singularities A and C and reverse singularity B a

nonsingular network can be formed. What is important to keep in mind

is that when the ground stations are generally distributed in space a

nonsingular network can be formed if at least six targets are

co-observed by at least four stations. In fact four stations and five

targets can uniquely define a second-order surface, and the sixth target

could make the network nonsingular if it is not located on this surface.

Once a nonsingular network has been realized, an extension of it

will result in a nonsingular one if singularity A is eliminated for any

additional station and if no targets are on the plane of the three

stations defining the Cartesian coordinate system.

2.2 DYNAMIC AND SEMIDYNAMIC MODE METHODS

In contrast to the geometric methods, the observed satellite positions in

the dynamic and in the semidynamic methods are not treated as auxiliary

independent points in space but rather they are constrained to lie in a

space curve (Schwarz, 1969). This curve should resemble within the

required degree of accuracy the satellite orbit under question. The

satellite orbits, on the other hand, are modeled either empirically or

dynamically or by combining both empirical and dynamical modeling.

Empirical modeling of satellite orbits was extensively used in the

early years of satellite geodesy since many of the model parameters

17
entering the equations of motion were not known with the degree of

accuracy required for precise geodetic work (Mueller, 1964b). These

methods are used today in very special circumstances and only in

combination with dynamic modeling (Tapley et al., 1985a).

Dynamic modeling results in three second-order differential

equations or six first-order differential equations referred to as

equations of motion of the satellite. These differential equations are

integrated either analytically (i.e., general perturbation methods) or

numerically (i.e., special perturbation methods) to generate the satellite

orbits.

In the general perturbation methods, the equations of motion of the

satellite are reformulated in terms of a set of orbital elements leading to

a set of differential equations which can be integrated analytically.

Unfortunately, a closed form analytical solution for the equations of

motion of the perturbed two-body problem does not exist. It is

possible, however, to obtain approximate analytical solutions either by

restricting the complexity of the perturbation model or by truncating

high power expansions (Kaula, 1961, 1966; Mueller, 1964b; Goad, 1977;

etc.) These solutions are approximate and in many cases cannot be used

for precise geodetic work. They are extremely useful, however, in order

to gain a keener insight into the effects of various perturbing forces on

the satellite orbits.

In the special perturbation methods, the equations of motion of the

satellite are integrated numerically (see Section 2.2.3). The main

advantage of these methods is that all the perturbing forces can be

accommodated to a high degree of accuracy. The special perturbation

methods, on the other hand, have proven to be computationally more

efficient, if one takes into consideration the high repetition rate of

recent geodetic observations (Rizos et al., 1985; Krakiwsky et al., 1985).

A combination of empirical and dynamic modeling is usually employed

either when the satellite orbits are integrated continuously over long

periods of time (i.e., two months or more) or when unexplained

perturbing accelerations are present. In the latter case the dynamic

models are supplemented with empirical models for the as yet not fully

18

II

li

II

I

I

I

I!

II

II

I

I

I

I

I

II

I

I

II
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

understood perturbing accelerations while in the former the empirical

models are employed to account for the accumulated along-track,

cross-track and radial errors (Tapley etaL, 1985a). For instance, the

draglike acceleration of Lageos' orbit which causes a decay of the

semimajor axis at a rate of 1 mm/d is modeled empirically (see Section

z.2.s).

The dynamic mode methods are further subdivided into semidynamic

(short-arc) and dynamic (long-arc) methods. There is no clear

distinction between these two terms and their exact meaning depends on

the investigator and on the kind of problem being analyzed.

In the present study the estimation of the baselines is performed in

the semidynamic mode environment. In this environment the lengths of

the arcs employed are relatively short (i.e., mostly up to three days and

very r rarely up to seven days) (see Section 4.4). A relatively short arc

is defined as having a maximum length over which the total modeling

error of a simple dynamic model is well below the noise level of the

observations (i.e., an order of magnitude or less). Consequently, with

this definition one may select a relatively simple dynamic model and then

determine the length of the arc, or one may choose the length of the

arc and then determine the required sophistication for the dynamic

models. With such a procedure the systematic errors caused by model

imperfections cannot accumulate up to a level that may corrupt the

semidynamic mode solution. The relatively short arcs, however, are not

stable in the sense that their position and orientation in space depends

primarily on the geometry of the observations. This instability may also

cause ill-conditioning of the normal equations (see Section 4.4).

Furthermore, relatively short arcs cannot be well tracked to bring

tracking sites of a global extent into a consistent satellite reference

frame. This implies that it is not possible to use semidynamic mode

methods for absolute position determination. Instead these methods can

be effectively used for baseline determinations (Latimer et al., 1977;

ChristodouUdis et al., 1981; Pavlis, 1982; Section 4.4). Baseline estimates

are even more accurate if the observables are insensitive to the

position, orientation, and the shape of the trajectory as is the case, for

19
instance, with the Simultaneous Range Difference observables (see

Section 4.4; Pavlis, 1982). With range observables, which are sensitive

not only to the position and orientation of the trajectory in space but

also to its shape, it is still possible to obtain accurate baselines if a

local support tracking network is available (Christodoulidis et al., 1981).

2.2.1 Simultaneous Range-Difference Semidynamic Mode Method

On the basis of the discussion presented in the previous section and

keeping in mind that our aim is to achieve highly accurate differential

positioning, we have chosen to use in this investigation the semidynamic

(short-arc} method formulated in the context of special perturbations

(see Section 2.2.3). Furthermore, the laser range observations have

been transformed to Simultaneous Range-Differences (SRDs), and

although differencing is a noise generating operation it is anticipated

that these observables are less affected by the biases in the orbit, the

reference frame and the observations themselves (Pavlis, 1982).

Using laser range observations to Lageos it is impossible to obtain

strictly simultaneous observations not only because Lageos is a passive

satellite but also because there will always exist synchronization errors

among the variou¢ observing stations. Therefore Simultaneous Range

Differences can only be obtained through an interpolation (see Chapter

3). More specifically, the observing stations are divided into pairs of

simultaneously observing stations. For each pair the station with the

most observations is interpolated to generate ranges at the observed

epochs of the alternate station. Finally the interstation distances for

each of the pairs involved are estimated by processing the generated

SRDs through a least squares adjustment formulated in the context of

the special perturbation methods as applied in the semidynamic mode

environment.

2.2.2 Mathemtical Modeling

The mathematical model for the Simultaneous Range Difference (SRD}

observable 6pj is obtained by subtracting the Euclidean ranges from

station 2 and station 1 to the simultaneously observed satellite position j

2O

I

I

I

I

I

I

l

I

I

I

I

I

I

I

I

I

I

I
(2-20)

I

I

!

I

I

I

I

!

I

where the vectors Sj - (Uj, Vj, Wj) T, Xl - (Ul, Vs, Wl) T and X2 = (u2,

v2, w2) T denote the Cartesian coordinates of the satellite position j,

ground station 1 and ground station 2 respectively. Since the SRDs are

invariant under any rigid body rotation the above vectors could be

expressed in any arbitrarily chosen Cartesian coordinate system. In the

present study, the vectors §j, Xs and X2 at epoch j are expressed in a

Cartesian coordinate system whose origin is conveniently chosen to

coincide with the center of mass of the earth, and its orientation is

aligned to that of the true-of-date system (Mueller, 1969).

The adjusted parameters, in any estimation procedure involving a

dynamic process, are transformed to a reference frame in which they

can be considered constant for a certain period of time. This period

should be long enough to allow for collection of a sufficient number of

observations needed for a reliable recovery of the adjusted parameters.

For this reason the ground station coordinates are transformed to a

terrestrial reference frame (TRF) while the coordinates of the satellite at

epoch j are transformed to a celestial reference frame (CRF) with the

help of the following formulas (see Section 2.2.4)

§j = N-PRj ; XI = STYi , i = i, 2 (2-21)

I

I

I

I

I

The quantities S, N and P designate the earth rotation, the nutation and

precession matrices respectively, while the vectors R j, Ys and Y2 denote

the inertial position vector of the satellite at the epoch j and the

earth-fixed position vectors of stations 1 and 2 respectively. The

inertial and earth-fixed frames correspond to the CRF and TRF frames

respectively as they are described in Section (2.2.4). Substituting

equation (2-21) into (2-20) one obtains

_: {Dj2T0j2} _-{DjTOjI}_Dj_i10jll (2-22)

where

21
Dj2 = I_Rj - STy2 and Djl = NPRj - STy1 (2-23)

The satellite position vector Rj is a function of an initial state

vector and a large number of parameters affecting the motion of the

satellite (i.e., potential coefficients, reflectivity, etc.). The choice of the

model parameters to be estimated depends on the data coverage and

distribution which in turn dictates the adopted lengths for continuous

integration of the satellite orbit (see Section 2.2.1). In the present

study, the shape of each of the satellite arcs involved is assumed known

and only its position and orientation in space is adjusted to "best" fit

the available data (see Section 4.4). Thus, the only adjusted parameters

inherent to the satellite position vector Rj are the components of the

initial state vector of the corresponding arc.

In the derivation of the observation equations, on the basis of the

equation (2-20), one needs the satellite position vectors at each of the

observing epochs together with their partial derivatives with respect to

the corresponding initial state vector. The former is obtained by

integrating the equations of motion of the satellite while the latter is

obtained by integrating the variational equations of state (see Section

2.2.3). The partial derivatives with respect to earth-fixed station

coordinates, also needed in the derivation of the observation equations,

are easily obtained by differentiating equation (2-20) (see also Pavlis,

1982). The resulting observation equations are used to obtain the

normal equations through a weighted least-squares adjustment (see

Section 2.2.6). The normal equations are subsequently solved to estimate

the initialstate vectors for each of the arcs involved together with the

earth-fixed coordinates which are finally transformed to interstation

distances.

In the present study the initial state vectors are treated as

"nuisance" parameters, and therefore one is not concerned with how well

each of those initial state vectors has been recovered as long as the a

posteriori variance of unit weight is close to unity. In fact, the reason

for using SRDs instead of ranges is to reduce the need for accurate

knowledge of the satellite orbits and yet to increase the potential for

22

I

ii

it

I

It

I

I

II

I

I

II

I

i

I

I

I

I

I
I

I

I

l

I

I

I

I

I

I

i

I

I

I

I

I

I

I

baseline estimation with an accuracy compatible to or even better than

that of the observations. This is possible because SRDs have the

potential to reduce the effects of biases caused not only by the orbital

model and the reference frames but also by eliminating uncorrectable

systematic errors affecting the laser range observations (Pavlis, 1982).

2.2.3 Orbit Determination with the Method of Special Perturbations

Dynamic and semidynamic methods, based on special perturbations,

require integration of the satellite's equations of motion. The degree of

sophistication in formulating these equations depends on the integration

length and the required accuracy.

Following the MERIT standards the relativistic perturbations are

ignored from the equations of motion (Melbourne et al., 1983).

Accordingly, ephemeris time (t) constitutes the independent variable in

the equations of motion. Up to 1983, ephemeris time was used as an

independent variable in the planetary equations of motion and therefore

in the construction of all the almanacs. Since January Oh 1984,

ephemeris time has been replaced by Terrestrial Dynamical Time (TDT)

and Barycentric Dynamical Time (TDB) (The Astronomical Almanac, 1984).

This was a necessity since data collected in interplanetary missions are

routinely processed in the context of the relativity theory (Moyer, 1971).

In this context TDT time corresponds to proper time (i.e.,time measured

by the observer's clock) while TDB corresponds to coordinate time (i.e.,

time measured at the barycenter of a motionless solar system in the

absence of all gravitational fields).

At each observing epoch j, the ephemeris time (tj} is computed from

the Universal Coordinated Time (UTCj} with the help of the following

formula

tj = 32_184 + [TAI - UTC]j + UTCj = TDT (after 1984) (2-24)

where

tj = ephemeris time at the epoch j

23
[TAI - UTC] j = no. of leap seconds at the epoch 3

Almanac, p. B5)

UTCj = Universal Coordinated Time at the epoch j.

(The Astronomical

Using the ephemeris time (t) as an independent variable, the Lageos

equations of motion take the following form (Cappellari et al., 1976;

Pavlis, 1982):

R - RpN + RNS + RTD + RSR + RAT (2-25)

Each of these accelerations is expressed relative to the center of mass

of the earth. More specifically

= d2 ( IT R = d-_ uj, v j, wj = total quasi-inertial Lageos acceleration at

the epoch j

RpM = gravitational acceleration due to point masses

RNs = gravitational acceleration due to nonsphericity of the gravitational potential

RTD = acceleration due to solid earth tidal effects

RSR = acceleration due to solar radiation pressure

RAT = Lageos empirical drag-like acceleration

The acceleration vector R is a function of an initial state vector and a

large number of parameters affecting the motion of the satellite. These

parameters pertain to the gravitational potential, to solar radiation

pressure, etc. As it was described in the previous section, the only

adjusted orbital parameters considered in this study are the initialstate

vectors of all the arcs involved. With these orbital parameters the

variational equations assume a very simple form (Pavlis, 1982):

Y(t) = A(t) Y(t) (2-26)

24

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

i

l

I

I

I

I

I

I

I

I

I

with initial conditions

Yo • [ I I 0 ].x6 (2-27)

where

and

Y(t) = [ aR(t)].

a (_(to) R(to)}

' 3X6

A(t) [

(2-28)

(2-29)

The matrix Y(t) is referred to as the state transition matrix and is used

to map the variations of the initial state into variations of the current

state.

Equations (2-25) and (2-26) can be integrated either by one-step or

by multi-step methods. In each integration step, the multistep methods

require fewer derivative evaluations than the one-step methods of

compatible accuracy. Fewer computations, on the other hand, not only

reduce the round-off errors but also require less computing time.

Furthermore, since these methods possess a larger number of parasitic

solutions they are more susceptible to instability problems.

The multistep algorithm used in our study employs a self starting,

variable step, variable order predictor-corrector mode of operation.

This mode selects the order automatically while the stepsize is subject

to accuracy requirements and numerical stability. Keeping the stepsize

constant, the predictor-corrector is reduced to an Adams-Bashforth

predictor of order q and to an Adams-Moulton corrector of order q + 1.

With this algorithm the second-order differential equations are

integrated directly without reducing them to a first-order system,

because a second-order set exhibits better numerical stability

characteristics. The described algorithm was developed and implemented

in computer coded form by Krogh (1969a, 1969b, 1973a, 1973b, 1974).

25
2.2.4 Reference Frames and Systems

Reference frames constitute realizations of reference systems. The

reference frames are used to describe the spatial relationships and the

temporal variations of objects on the Earth (i.e., terrestrial frames} and

in space including the Earth (i.e., celestial frames} (Kovalesky and

MueUer, 1981). A reference system consists of an underlying principle

and all those elements (e.g., physical environment, theories and

constants) that are necessary to accomplish its realization. The elements

of a system, depending on the application and the accuracy

requirements, are selectively chosen and therefore the term

"conventional system" is often used to designate the selection process

that is usually involved in any realization of a reference system.

In this context the IAU/IUGG MERIT and COTES Joint Working

Group recommended the following concepts in regard to reference

systems and frames (Wilkins and Mueller, 1986):

The Conventional Terrestrial Reference System (CTRS) be

defined by a set of designated reference stations, theories and

constants [necessary elements], chosen so that there is no net

rotation or translation between the reference frame and the

surface of the earth [underlying principle]. The frame is to be

realized by a set of positions and motions of the designated

reference stations.

The Conventional Celestial Reference System (CCRS) be defined

by a set of designated extragalactic radio sources, theories and

constants [necessary elements], chosen so there is no net rotation

between the reference frame and the set of the radio sources

[underlying principle]. The frame is to be defined by the

positions and motions of the designated radio sources. The origin

of the frame is to be the barycentre of the solar system.

The above concepts are to be incorporated in the operation of the new

International Earth Rotation Service (IERS). This service is scheduled

to start operating as of January Oh, 1988 (Mueller and Wilkins, 1986}.

In the Newtonlan framework, the reference frame implied by a CCRS

can be considered as being an ideal inertial frame in the sense that the

time is homogeneous and the space described by this frame is

homogeneous and isotropic (Landau et al., 1960}. In the general

26

I

I

l

l

l

l

l

l

II

I

l

l

I

I

I

I

I

I
I

I

I

l

I

i

I

I

I

I

I

l

I

I

I

I

I

I

relativistic framework, on the other hand, the reference frame realizing

a CCRS is aimed to describe the curved space time for which a global

inertial reference frame does not exist (Moritz, 1979; Fukushima, 1986).'

This is the reason why, in the above recommendations, the term

"inertial" has been dropped and the term "celestial" has been used

instead.

For precise geodetic work, these seemingly conceptual differences

manifest themselves when 10-8 or 10 -9 accuracies are sought.

Therefore, when working at such accuracy levels, care should be taken

to account for relativistic effects either by using Newtonian formalism

with small corrections (Moritz, 1979) or by formulating the problem

entirely in the general relativistic framework {Fukushima, 1986). In the

present study, in accordance with the MERIT standards and since the

obtained accuracies hardly reach the 10-8 level, we have used the

Newtonian formalism to formulate the equations of motion of the satellite.

The status today in terms of reference systems and frames is

confusing because the user community employs a variety of different

celestial systems (i.e., extragalactic radio source systems, stellar

systems, dynamical systems, etc.) and a variety of different terrestrial

systems as well (i.e., BIH terrestrial reference system, CSR terrestrial

reference systems, etc.). Investigations, however, are currently

underway with the objective of linking all of the available terrestrial

systems into a unified terrestrial system referred to as "BIH Terrestrial

Reference System" {BTS) (Boucher and Altamimi, 1985, 1986). Linking

different celestial systems, through their frames, into an ideal celestial

frame is not an easy task not only because of lack of collocations but

also because daily polar motion resolution is necessary {Mueller, 1985).

This kind of resolution is not achievable by the satellite related systems

due mainly to the deficiencies in nutation theory (Himwich and Harder,

1986) and to inadequate observational coverage.

The choice of ideal terrestrial and celestial frames is not important

in baseline estimation. It is important, however, to consistently link the

involved terrestrial and celestial reference frames, by choosing the

27
appropriate set of transformation parameters. Effective choice of the

transformation parameters would only assure a reliable recovery of the

relative geometry of the observations since these parameters cannot be

effectively recovered in a semidynamic mode environment. The relative

geometry of the observations manifests the way the satellite arcs are

related to the observing stations. Reliable recovery of the relative

geometry, on the other hand, results in accurate baseline estimation

simply because baselines are estimable quantities. Estimable quantities

are molded by the geometric and dynamic characteristics creating the

problem under question.

In the present study, the Terrestrial Reference Frame (TRF) is

implied by the gravity field used to integrate the equations of motion

and by the adopted polar motion series. The origin of the TRF frame,

relative to the center of mass of the earth, is defined by the potential

coefficients C1o, C11 and S lz, while its orientation is primarily defined

by the potential coefficients C21 and $2_ as well as C22 and $22. More

specifically, the coefficients C2_ and $2_ define the orientation of the

third axis, while C22 and $22 define the orientation of the first axis.

The orientation of the first axis, however, is weakly defined because the

Earth's equatorial moments are nearly equal (i.e., C22 _ $22).

The modified GEML2 gravity field, proposed by the MERIT standards,

has been replaced in this study by the PGS1680 gravity field

(Christodoulidis et al., 1985). This violation of the MERIT standards was

necessary to make the gravity field consistent with the adopted BIH

polar motion series. The modified GEML2 gravity field has all its

coefficients but S2t and C21 equal to the coefficients of the GEML2

gravity field (Lerch et al., 1985). C2t and $2, have been modified to be

consistent with the mean BIH polar motion values, computed over a

complete wobble cycle which lasts from 6.5 to 7 years (Melbourne et aL,

1983). Modifying only the C2_ and $2, has caused inconsistencies as to

what frame the coefficients of the modified GEML2 field refer. As a

result, the PGS1680 gravity field was developed in order to avoid these

inconsistencies and the resulting confusion as well. In the development

of this field the coefficients $21 and C2x were constrained to the BIH

28

I

I

I

I

I

I

I

i

I

I

I

i

I

I

I

I

i

I
I

I

I

I

I

I

I

I

i

I

I

I

I

I

I

I

I

I

implied values while the remaining coefficients were free adjusting

(Christodoulidis et al., 1985).

The coefficients C,o, C** and S** of the PGS1680 gravity field are

zero thereby imposing the origin of the TRF frame to coincide with the

center of mass of the earth. The computed UT1 time, on the other

hand, is assumed to be consistent with the x axis of the implied TRF

frame (see Section 4.3).

The Celestial Reference Frame (CRF) employed for the integration of

the equations of motion, is realized from the implied TRF frame through

the following transformation (Mueller, 1969}:

m m

(CRF) : (SNp)T(TRF)

with

(2-30)

S : R2 (-xp) R, (-yp) R3 (GAST) (2-31)

where according to the MERIT standards the following quantities have

been used.

P Precession matrix based on the IAU (1976) system of astronomical

constants (Lieske, 1979)

N Nutation matrix based on the 1980 IAU theory of nutation (Wahr,

1981a). This matrix implies a pole whose nearly diurnal space-fixed

and earth-fixed motions vanish. This pole is referred to as the

Celestial Ephemeris Pole (CEP) (Mueller, 1981; Moritz and Mueller,

1987)

GAST _ GMST (OhUTI) + f(UTI) + EQ.E (2-32)

GMST(OhUT1) = 6 h 41 m 50_54841 + 8640184_812866 Tu + 0_093104 Tu 2

- 6_2 x 10 -6 Tu 3 (2-33)

where

conversion factor from Universal time to sidereal time

1.002737909350795 + 5.9006xI0 -.I Tu - 5.9xi0 -*s Tu2

Tu =

UTI =

Julian centuries elapsed from J2000.0

UTC(USNO) + (UTI - UTIR) + [UTIR - UTC(BIH)] +

+ [UTC(BIH) - UTC(USNO)] (2-34)

29
UTC(USNO)--USNO Universal Coordinated time (time scale used to

time tag the observations of the GLTN stations)

UT1- UT1R = Tidal variations of UT1 caused by zonal tides with

periods up to 35 days (BIH Annual Report 1981 onwards, Table B1)

UT1R- UTC(BIH) = Variations of the regularized UT1 (i.e., UT1R)

from the UTC(BIH) (BIH Annual Reports, Table 8)

UTC(BIH) - UTC(USNO) --Variations

UTC(BIH)

of UTC(USNO) in relation to

EO.E = A@ • cos (_ + A_) (2-35)

A_ = Nutation in longitude computed from the 1980 IAU nutation

theory

= Mean obliquity of the ecliptic

= 23" 26' 217448 - 4678150Tu - 0700059Tu 2 + 07001813Tu 3 (2-36)

_ = Nutation in obliquity computed from the 1980 IAU nutation

theory

The CEP pole positions Xp and yp in equation (2-31) have been taken

from the smoothed values of Circular D (BIH Annual Report, 1983, 1984,

Table 7). These pole positions are referenced to the 1979 BIH system

during the first period of the MERIT Main Campaign (Sept. 1983 - Dec.

1983), while during the remaining period of the campaign (Jan. 1984 -

Dec. 1984) they are referenced to the BIH Terrestrial System (BTS). Our

study in not affected by this transition because in shifting from the

1979 BIH system to the BTS system a nonrotation constraint was applied

to assure the continuity of the BIH system (Boucher and Feissel, 1984).

2.2.5 Orbital Model

The set of elements necessary to determine a satellite orbit constitutes

the orbital model of the satellite. Thus, an orbital model consists of all

those elements that are essential to formulate and integrate the

3O

I

I

I

I

I

I

I

I

l

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

equations of motion (i.e.,perturbations to be considered and associated

assumptions, initialconditions, etc.).

The Lageos state vectors required in the evaluation of the

observation equations (Section 2.2.6) are obtained by numerically

integrating the equations of motion (2-25). In this equation the inertial

accelerations are expressed relative to the geocenter. The components,

however, of Lageos' inertial acceleration caused by the nonsphericity of

the earth are evaluated in the TRF frame and subsequently are

transformed into the corresponding CRF frame, while the components of

the remaining inertial accelerations are directly expressed in the

corresponding CRF frame. The two-step procedure used to evaluate the

Lageos inertial acceleration caused by the nonsphericity of the earth is

necessary because the gravity potential coefficients are conveniently

expressed in a TRF frame. Expressing these coefficients directly in a

CRF frame would make them time dependent and therefore a potential

source of unnecessary complications.

The aim of the present study is not to estimate the Lageos orbit

with the highest accuracy but rather to model it as simply as possible

and yet be able to recover the baselines with an accuracy compatible to

or even better than that of the observations. With this in mind the

MERIT standards have been violated whenever the proposed model is

complicated and cumbersome to incorporate into the solution. In such

cases a simpler model has been adopted. It turns out, however, that in

some cases the employed orbital model could be further simplified

without affecting the accuracy of the recovered baselines (Section 4.6).

2.2.5.1 Point mass gravitational acceleration. The point mass Lageos

gravitational acceleration based on the effects of the three major

perturbing bodies (Earth (E), Moon (M) and Sun (S)) and expressed

relative to the geocenter takes the following form (Cappellari et al.,

1976; Pavlis, 1982)

(2-37)

31
where

P_, _, P_ =: Lageos position vectors relative to the center of

mass of the Earth, the Moon and the Sun respectively

=: Earth position vectors relative to the center of mass

of the Moon and the Sun respectively

MM Ms

ME ' ME

=: ratios of lunar and solar masses to the mass of the

Earth

GM E =: geocentric gravitational constant

For the evaluation of equation (2-37), one needs the Lageos geocentric

position vector as well as the geocentric position vectors of the Moon

and the Sun respectively. The Lageos geocentric position vectors are

obtained from the numerical integration of the equations of motion while

the heliocentric position vector of the Earth and the geocentric position

vector of the Moon are calculated from the information supplied by the

DE/LE200 lunar planetary ephemeris (Standish, 1981). This ephemeris is

disseminated in terms of Chebychev coefficients. These coefficients can

only be used to calculate the geocentric positions of the Moon and the

barycentric positions of remaining planets and the Sun. With this

information, however, one can very easily calculate the position of any

desired planet with respect to any of the remaining planets and to the

sun as well.

The reference frame implied from the computed coordinates of the

planets has been accurately adjusted to the dynamical equinox J2000.O

(ibid.) The Chebychev coefficients of the DE/LE200 lunar planetary

ephemeris are based on the planetary coordinates estimated through the

numerical integration process involved in the adjustment of

interplanetary observations collected over a long period of time (ibid.).

In this adjustment, the planetary equations of motion were formulated on

the basis of the isotropic, parametrized post-Newtonian (PPN) n-body

metric (Moyer, 1971). The independent variable in the PPN metric is the

Barycentric Dynamical Time (i.e., coordinate time), and therefore this

time scale should be used as an entry to the DE/LE200 ephemeris.

32

1

I

I

1

I

1

1

I

1

1

I

1

I

1

1

I

1

I
I

I

I

i

I

I

I

I

I

I

I

I

I

I

I

I

I

I

The TDB time at any epoch j is computed from the ephemeris time of the

same epoch via the following formula

TDBj = tj + AT

where

(2-38)

TDBj = Barycentric Dynamical Time at the epoch j

tj = ephemeris time at epoch j, obtained from equation (2-24)

In the PPN framework the ephemeris time coincides numerically but not

conceptually with the Terrestrial Dynamical Time (i.e., proper time).

Thus, TDB time at any epoch is obtained by adding to the ephemeris

time of the same epoch a small correction AT. This correction accounts

for the general relativistic effects involved in the transformation of

proper time (i.e., TDT time) to coordinate time (i.e., TDB time). An

approximate value for the correction AT is given by the following

formula (Astronomical Almanac, 1984)

AT = 0_001658 sin (g) + 0.000014 sin (2g) (2-39)

where

g = 357:53 + 35999.05 Tu (2-40)

In both of these equations (2-39 and 2-40), higher-order terms have

been neglected, g designates the mean anomaly of the Earth in its orbit,

and Tu designates the Julian centuries elapsed since J2000.0.

To complete the evaluation of equation (2-37), one still needs the

ratios of the lunar and solar masses to the mass of the earth as well as

the geocentric gravitational constant. For the mass ratios, we have

used the values recommended by the MERIT standards, but for the

geocentric gravitational constant the value estimated simultaneously with

the potential coefficients of the PGS1680 gravity field has been used:

G_4E = 3.986004359 x 1014 mS/s 2 (2-41)

33
The scale in the range dynamic mode methods is implied not only by the

adopted value of the geocentric gravitational constant through the

modified Kepler's third law but also through the speed of light used to

convert time measurements to range measurements. Thus, the adopted

value of the geocentric gravitational constant should also be consistent

with the speed of light implicit in the range observations. In the

present study we have used the speed of light proposed by the MERIT

standards (i.e.,c - 299,792,458 m/s) (Lerch et al., 1985; Christodoulidis

et al., 1985).

The associated partial derivatives of equation (2-37) contributing to

the variational equations of state (i.e., to matrix Y(t), equation (2-28))

are given in (Cappellari et al., 1976, eq. 4-21; Pavlis 1982, eq. 13).

2.2.5.2 Gravitational acceleration due to nonsphericity of the

gravitational potential. The inertial acceleration induced on the satellite

by the nonsphericity of the earth is obtained via the gradient of the

perturbing potential. The perturbing potential is a scalar function

describing the nonspherical part of the geopotential in terms of an

infinite spherical harmonic series (Heiskanen and Moritz, 1967):

VNS(r,¢,X) = (_dE[ _ {a--_Er}n _. [Cn. cos(m_) + Sn. sin(m_)]Pnm(Sin¢)}

r l_n= 2 m=O

(2-42)

The zero-degree harmonic has been modeled in equation (2-37) and

therefore is not included in the above equation. The first-degree

harmonics are also not included because the origin of the PGS1680

implied (TRF) coincides with the the center of mass of the earth (see

Section 2.2.4). With the gradient of the perturbing potential the

components of Lageos' inertial acceleration, caused by the nonsphericity

of the earth, are expressed in the PGS1680 implied (TRF) frame.

Incorporation of this acceleration into the equations of motion (2-25)

requires transformation of its components from the (TRF) frame into the

corresponding (CRF) frame via the transformation equation (2-30)

34

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

(2-43)

where

= iaVNs aVNs aVNs]

r : L_x ' a---y- ' -'_-z J (2-44)

is the gradient of the perturbing potential function VNS (equation 2-42).

The expressions for the partial derivatives of the perturbing potential

function VNS are given in (Cappellari et al., 1976). The radius (aE) of

the reference sphere, also needed in the evaluation of the perturbing

potential, is the same with the radius employed in the estimation of the

PGS1680 gravity field (i.e., a E = 6378144.11 m). For our study we have

truncated this field at degree and order 12 because perturbations

caused by higher harmonics over a two-week period contaminate the

computed SRD observables with errors having magnitude well below the

noise level of the SRD quasi-observables (see Section 4.6). Furthermore,

a nonvariant nature of the coefficients C2z and S2z has been adopted,

although it is well known that these two coefficients are largely affected

by the forced diurnal motion of the figure axis caused by the

nonrigidity of the earth (Moritz and Mueller, 1987).

The associated partial derivatives of equation (2-42) contributing to

the variational equations of state are given in {Cappellari et al., 1976,

eq. 4-54; Pavlis, 1982, eq. 21).

2.2.5.3 Lageos tidal inertial acceleration. According to the MERIT

standards, the tidally induced variations in the earth's external potential

should be incorporated in the orbital model as variations in the

geopotential coefficients (Melbourne et al., 1983; Eanes et al, 1983). In

order to save computing time a two-step procedure is proposed to carry

out the implementation of these variations (ibid). In the first step the

variations of the geopotential coefficients are computed on the basis of a

nominal frequency independent Love number kn, while in the second

step these variations are corrected to account for the frequency

dependent nature of the nominal Love number k n.

35
The frequency dependent variations of the Love number k, have been

estimated for an elliptical, rotating, elastic fluid outer core and solid

inner core, oceanless earth {Wahr, 1981b). Although the proposed

two-step procedure is computationally more effective than any one-step

procedure, it is not appropriate for our investigation because one still

needs to evaluate many trigonometric functions at each of the observing

epochs (Melbourne et al., 1983}. This however, not only would make the

orbital model more complicated but also it would make the SRD method

computationally less efficient. Therefore, it was decided to compute the

tidally induced space potential by assuming a solid earth (i.e., oceans

not included} which exhibits the same elastic response over all possible

orders within a certain degree (Diamante et al., 1972; Pavlis, 1982). With

such an earth model the tidally induced potential on the surface of the

earth takes the following form (Diamante et al., 1972; Goad, 1977; Pavlis,

1982}:

UTD = n=2_ kn UTn (aE) (2-45)

where k, is the nominal Love number of degree n and UT (aE) is the n th

surface harmonic of the tidal potential. Solving the Dirichlet problem

the tidally induced space potential is obtained:

UT0 = aE}n+z ,==_ _ k, UT,(a[) (2-46)

where I1_1 is the norm of the Lageos position vector expressed relative

to the center of mass of the earth. To the first order, terms with n > 2

in equation (2-46) can be neglected (Diamante et al., 1972), and therefore

this equation takes the following form

a_

UT o :_'_ k2 UT2 (a E)

I_1

(2-47)

where

36

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
i

i

I

I

I

I

I

i

I

I

I

I

I

i

I

I

I

o- [(, ]

UT= (aE) = a_ 3 l_b _ 2 _ I (2-48)

21_.1" _1

The quantity k= denotes the second-degree Love number while Mb

denotes either the lunar or the solar mass, or for that matter the mass

of any planet that is considered to be a disturbing body. In this study

the Moon and the Sun have been considered as the only disturbing

bodies. The vectors I_ and l_b designate the geocentric position vectors

of the satellite and the disturbing body respectively. The components

of these vectors are expressed in the corresponding (CRF) frame. In

this frame the tidally induced acceleration on Lageos takes the following

form (Diamante et al., 1972; Pavlis, 1982):

= 3 _._.A_" a_ =)

._o-_k= I_.1" I-_ [{1-5(m" _) u+2CUb U)Ub] (2-49)

where

Ub = _ and 5 = (2-5o)

To account for a phase lag produced by the earth's dissipative forces,

the vector Rb in equation (2-50) has been replaced by another vector

1_.*b. This vector is obtained from the vector Rb via the following

transformation

(2-51)

where 6 L (-0*.35) is the phase lag. The R3 rotation is performed about

the third axis of the corresponding CRF frame (see Section 2.2.4). In

equations (2-47) and (2-49) the value 0.29 was adopted for the

second-degree Love number k=. This value is different from that

proposed by the MERIT standards (k= - 0.30). This deviation, although

not of much importance, is justified since the tidal corrections applied

in the estimation of the PGS1680 gravity field were based on the altered

value of k= (i.e.,k= -- 0.29). The permanent tidal deformation affecting

37
I

the potential coefficient C2o is inherently present in equation (2-49).

Consistent incorporation of this equation in the equations of motion

requires that the permanent tidal deformation is not included in the

PGS1680 C2o value. This, however, seems to be the case for the

PGS1680 gravity field (Melbourne et aL, 1983; Christodoulidis et al.,

1985). Furthermore, the ocean tidal perturbations are not included in the

orbital model of Lageos not only because they are small (i.e.,one order

of magnitude smaller than the solid earth tidal effects (Section 4.6)) but

also because their evaluation would increase the bulk of the

computations considerably.

The contribution of the tidally induced acceleration to the

variational equations of state is given in (Pavlis, 1982, eq. 30). In that

equation the vector Rb should be replaced by the vector R*b from

equation (2-51).

2.2.5.4 Lageos solar radiation pressure acceleration. The

acceleration induced on Lageos due to photon momentum transfer is

referred to as solar radiation pressure acceleration, and it is given by

the following formula (El'Yasberg, 1967; Cappellari et al., 1976; Pavlis,

1982).

•" s[Au}

I

I

I

I

I

I

I

I

I

I

A R- RS

C, " _ " IRs - RI (2-52) I

The eclipse factor 7 assumes the values zero, one, or any other value in

between depending on whether the satellite is in complete shadow

(umbra), in sunlight, or in partial shadow (penumbra) respectively. In

our study the eclipse factor 7 is determined by a simple cylindrical

model (Cappellari et al., 1976; Pavlis, 1982). This model is easy to

incorporate into the equations of motion, but it does not differentiate

between umbra and penumbra regions. A full model for the earth's

shadow, as proposed by the MERIT standards, would increase the bulk

of computations thereby complicating the solution. It is rather doubtful

38

I

I

I

I

I

I
I

!

I

I

I

I

I

!t,

I

I

I

I

I

I

I

I

I

i

that this complication would make any difference. The mean solar flux S

is the amount of photon energy flow through a unit surface per unit

time at a distance of one astronomical unit (AU) (i.e., AU = 1.4959787066

x 10 .1 meters) (Melbourne et al., 1983). The ratio (S/C) is the photon

momentum transfer to a unit surface per unit time at a distance of one

astronomical unit. The value (4.5605 x 10 -6 N/M 2) was used in our

study for the ratio (S/C) as proposed by the MERIT standards. The

position vectors I_ and t_s in equation (2-52) designate the geocentric

position vectors of Lageos and the Sun respectively. The reflectivity

coefficient (C,) depends not only on the mechanism of light reflection

but also on the thermal emission distribution of the satellite surface.

The monthly values for the reflectivity coefficient (C,) are shown in

Table 1 for the entire MERIT campaign. These values have been

estimated together with other parameters in the adjustment of the MERIT

laser range data performed by the GEODYN II programs (Pavlis, 1986,

private communication). In the present study we have used the

reflectivity coefficient values listed in Table 1 instead of using the

value proposed by the MERIT standards.

Table 1 Lageos Along-Track Acceleration and Its Reflectivity

Coefficients

Magnitude of Lageos Along-Track

Acceleration x 10-*2 m/s a

Lageos Heflectivity

Coefficients

Sep. 1983 -2.909 1.141

Oct. " -3.549 1.136

Nov. " -3.893 1.135

Dec. " -3.825 1.133

Jan. 1984 -4.343 1.136

Feb. " -4.319 1.134

Mar. " -4.065 1.109

Apr. " -3.550 1.057

May " -3.189 1.096

Jun. " -3.393 1.139

Jul. " -3.928 1.132

Aug. " -2.946 1.126

Sep. " -2.301 1.126

Oct. " -2.524 1.126

39
This deviation of the MERIT standards, although plausible, will not affect

the accuracy of the estimated baselines. The effective area (A) of the

surface normal to the incident light is given, for a spherical satellite

like Lageos, from (El'Yasberg, 1967)

A = _Rs= (2-53)

where R s is Lageos' radius (i.e., R s = 0.30 m) (Melbourne et al., 1983).

Finally, the Lageos mass (M) of 407 kg has been used in the evaluation

of equation (2-52).

The associated partial derivatives of the solar radiation pressure

acceleration contributing to the variational equations of state are given

in (Cappellari et al., 1976, eq. 4-161 and 4-162; Pavlis, 1982, eq. 25 and

26).

According to the MERIT standards the inertial acceleration induced

on the satellite due to the diffused reradiated light from the earth

(earth a|bedo) is not included in the orbital model of the Lageos

satellite.

2.2.5.5 Lageos alon_-track empirical acceleration. Ever since the

launch of the Lageos satellite it has been observed that its semimajor

axis decreases at a rate of 1 mm/day. This has been traced to an

unexpected and still physically unmodelled along-track acceleration

acting on the Lageos satellite. Attempts to explain the origin of this

mysterious acceleration have either totally or partially failed. These

attempts are based on a variety of possible causes ranging from

assuming helium concentrations at satellite altitudes (Rubincam, 1980) to

considering the solar eclipses (Rubincam et al., 1985). Although all of

these attempts have partially failed, it is quite clear that this

acceleration is the result of a combined effect caused by the

asymmetries of the earth's albedo and by the charged particles traveling

in the vicinity of Lageos (Smith et al., 1985; Alfonso et al., 1985). Since

the physical process producing this acceleration is unknown, its

modeling has been accomplished with an empirical model. This model

4O

l

I

I

!

i

!

I

I

I

I

I

I

l

I

i

I

!

I
I

!

I

I

I

i

I

I

I

!

i

!

I

i

I

I

I

I

assumes a "drag-like" force acting against the satellite motion (Pavlis,

1982):

"" R

RAT = - o{ _'_ (2-54)

where _ is the magnitude of the along-track acceleration. The monthly

magnitudes of this acceleration are also listed in Table 1 for the entire

MERIT campaign. These values have been estimated with the GEODYN II

program (Pavlis, 1982, private communication).

Contributions to the variational equations of state due to this

acceleration are neglected because of their small magnitudes.

2.2.6 Normal Equations

The observation equations for the SRD quasi-observables are readily

obtained through a Taylor series expansion of equation (2-22). In this

expansion only the zero- and first-order terms are retained while all of

the remaining higher-order terms are neglected. The expansion is

performed about the approximate earth-fixed station coordinates and the

celestial initialstate vector of the corresponding arc:

Vj = AjXj + Lj (2-55)

where

(Aj),.,= = [ Bj i cj ],.,= (2-5s)

a6pj ]

(sj),x, = a(_,,_=) ,.,

(2-57)

[Yj(t)]3x6

(2-58)

41

I
Xj = (dY, , dY2 , dl_o , dRo),x,2 _ T (2-59)

I.,j = the computed minus the generated SRD quasi-observable

Vj = residual corresponding to the jth SRD observable

The matrices (Bj),x6 and (Sj),x3 are readily obtained by differentiating

equation (2-22), while the state transition matrix [Yj(t)]3x6 is obtained

in the numerical integration process of the variational equations of

state. The celestial satellite coordinates at the epoch j needed to

evaluate the vectors (Bj),x6 and (S3),x_ as well as the scalar Lj are

obtained from the numerical integration of Lageos' equations of motion

(see Sections 2.2.3 through 2.2.5). The adjusted parameter vector Xj

contains corrections to the earth-fixed approximate coordinates of

stations 1 and 2 (i.e.,dY, and dY2) and to the corresponding celestial

initial state vector (i.e., ch_o and dRo). Extension of equations

(2-55)-(2-59) to include all the available SRD observables and all the

observed satellite arcs leads to the following equations

= A*X + L (2-60)

where

A*-- [B* i C*] (2-61)

B* Fa6_

= LWJ

nx3k

(2-62)

= (S)nx._ [Y*(t)]3_,6m (2-63)

(2-64)

42

I

!

I

!

!

I

I

!

I

!

I

I

I

I

I

I

I

i
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

= vector containing the computed minus the SRD quasi-observables

= residual vector

Vector 6p contains all the available SRD observables, vector Y the

Cartesian coordinates of all the observing stations, vector 7. the

Cartesian coordinates of the observed satellite positions, and vectors

Iwo, I_o the initialstate vectors of all the arcs involved. The adjusted

A

m

parameter vector X contains the corrections to the approximate

earth-fixed coordinates of the observing stations (i.e.,dY) together with

the corrections to the initialstate vectors of the observed satellite arcs

(i.e., d(IWo, IITo)). The integers _ and n denote the number of the

observed satellite positions and the number of observations, while the

integers k and m denote the number of the observing stations and the

number of the observed satellitearcs respectively.

A close examination of equations (2-62) and (2-63) reveals that the

submatrice8 (B_,xsk) and (Snx_) would be exactly the same even if the

satellite positions were treated as auxiliary independent points in space

(i.e., geometric approach). The constraints imposed on the observed

satellite positions to lie in the corresponding satellite arcs are applied

through the state transition matrix (Y*(t)sSx6m).

Singularity A (see Section 2.1.3) affects the dynamic and geometric

solutions in exactly the same way because the submatrix B* in equation

(2-61) is the same for both the geometric and the dynamic approach.

With SRD observables singularity A occurs not only from the resulting

geometry of one station and its observed targets (see Section 2.1.3) but

also from the geometry of two coobserving stations and their observed

targets. Singularity B or singularity C cannot exist in a dynamic

solution because the structure of the matrix (S, x3_) is altered by its

multiplication with the state transition matrix (Y*{t)3_x_m). The

alteration of the matrix {Snxsl} not only differentiates the dynamic from

the geometric approach but also furnishes the dynamic approach with

better stability characteristics {see Sections 4.2 and 4.4).

Taking into consideration that the state transition matrix is different

from epoch to epoch, one can readily prove that in the absence of

43

I
singularity A the design matrix A* (apart from the ill-defined origin of

longitudes) is nonsingular if the adjusted parameter vector consists only

of corrections to the approximate station coordinates and to the initial

state vectors. This is not surprising because the TRF frame, with an

ill-defined origin of longitudes, is implied by the PGS1680 gravity field

while the CRF frame is subsequently realized via the transformation

equation (2-30). Including polar motion and/or variations of UT1 in the

adjusted parameter vector results in an extremely ill-conditioned design

matrix A* because polar motion and station coordinate are nearly

inseparable, while variations in UT1 and in the satellite node are

inseparable parameters as well (Van Gelder, 1978; Pavlis, 1982).

In the present study, polar motion and the variations in UT1 are

not included in the adjusted parameter vector. Thus, after resolving

the problem of the ill-defined origin of longitudes (see Section 4.4) we

proceed with the formation and the solution of the normal equations.

Using the same arguments as in Section 2.1.2 and the observation

equations (2-60), the normal equations take the following form (Uotiia,

1987)

(A*TPA* + Px)X + A*TPL = 0 (2-65)

where P and Px are the weight matrices associated with the SRD

observables and the adjusted parameter vector respectively. Since the

weight matrix P is diagonal the normals are formed sequentially through

the following formula

A*TPA* + Px = _ A]A] + Px (2-66)

where Aj (i.e., equation (2-56)) is the jth row of the design matrix A*, n

is the total number of the SRD observables and _ is the variance for

the jth SRD observable. These variances are computed via the following

formula

44

I

I

I

I

I

I

I

I

I

I

i

I

I

I

I

I

I

I
I

I

I

I

I

|

I

I

I

I

I

I

I

I

I

I

I

I

I

_r_ = _r_1 + _2 {j = l, , n) (2-67)

where ¢_1 and ¢_2 are the variances of the actually observed and the

interpolated ranges respectively. The inversion of the normal equation

matrix (2-66) was obtained with the Cholesky algorithm (Uotila, 1967).

45
Chapter 3

GENERATION OF THE OBSERVABLES

This chapter starts with a description of the SLR system in an attempt

to identify and understand the systematic errors affecting the laser

ranges. It continues with a description of the data set employed in this

investigation and finally ends with the generation of the simultaneous

range and SRD observables. These two observables constitute the input

to the geometric and SRD methods respectively.

3.1 SATELLITE LASER RANGING

A satellite laser ranging system consists of three basic components:

(i) the ground segment,

(ii) the atmospheric channel, and

(iii) the spaceborne segment.

The ground segment consists of a global network of fixed and

highly mobile satellite laser ranging stations forming a network

configured to allow measurements of the plate tectonic motions (Coates et

al., 1985). Tectonic plate motions are essential in understanding the

geodynamic processes necessary for earthquake and volcano erruption

predictions. Each of the stations in the network is equipped with the

necessary hardware to produce, emit, receive and measure the roundtrip flight time of very short laser pulses to a retroreflector equipped

artificial satellite such as LAGEOS.

The atmospheric channel is the optical path followed by a laser

pulse in its round trip from the station to the satellite.

The spaceborne segment consists of approximately 14 retroreflector

equipped satellites (Degnan, 1985). For geodesy and geodynamics Lageos

46

I

I

I

I

I

|

I

I

I

I

I

I

I

I

I

I

l

I
I

I

I

I

I

i

I

I

I

I

I

i

I

I

I

I

I

I

is an example of such a satellite in orbit high enough not to be

influenced by the difficult to model high frequency variations of the

gravity field and the atmospheric drag but yet at low enough altitude to

assure good signal returns to the tracking stations. Therefore, the

propagation of the orbital errors in the estimated geodetic parameters is

substantially reduced. This error reduction is very important because

variations in certain geodetic parameters such as baselines, polar motion,

and length of day are routinely used in understanding the mechanisms

driving geodynamic processes.

In the operational environment, depending on the technology

employed and the models used, each component of the satellite laser

ranging system will contribute in part to the total error affecting the

inferred geometric range. The next section contains, for each component

of the SLR system, a brief discussion of its operational principles, the

error sources, their status during the MERIT Main Campaign and the

future possibility of either reducing or eliminating them.

3.2 SATELLITE LASER RANGING SYSTEM, ITS COMPONENTS AND THEIR

CONTRIBUTION TO THE TOTAL ERROR BUDGET

3.2.1 Hardware of the Ground Segment

For each satellite ranging system, the hardware of the ground segment

consists of the laser transmitter, the laser receiver, their transmittimg

and receiving optics, the timing subsystem and the computer.

The /aser transmitter in most of the modern laser systems consists

of a mode-locked Nd:YAG laser oscillator followed by one or more Nd:YAG

laser amplifiers. The name Nd:YAG is derived from the crystal used in

the light amplification by stimulated emission of radiation {i.e., lasing

process) which is a YAG crystal {Yttrium aluminum garnet : Y3AlsO2)

"doped" with Neodymium {Johnson et al., 1978). Since the mode-locked

Nd:YAG lasers operate in a single spatial mode they are not affected by

"wavefront-distortion" errors (Degnan, 1985). The biases introduced by

47
the mode-locked transmitters are of the order of subcentimeter level.

The crystals used in the lasing process, for the stations participating in

the MERIT Main Campaign, are reported in the SLR coordinator's report

and its updates (Schutz, 1983b).

The /aser receiver is designed to measure the round-trip flight time

of the laser pulse to a retroreflector equipped satellite. This time

interval is multiplied by the speed of light and divided by two to infer

the optical range from the station to the satellite. The basic elements of

a laser receiver are the photomultiplier, the discriminator and the time

interval unit.

The photomul_iplier is a device used to detect the incoming laser

pulse. Its principle of operation is based on the photoelectric effect

(Halliday et al., 1962; Drain, 1980). Most of the SLR systems

participating in the MERIT Main Campaign made use of the conventional

type photomultipliers referred to as dynode-chain photomultipliers

(Degnan, 1985). The time it takes for the photoelectrons to propagate

from the photocathode to the anode via the dynodes is called transit

time. If the transit time were constant it could be completely accounted

for through either calibration or common channel procedures (ibid).

Variations in the transit time, referred to as transit time jitter,

influence the inferred ranges by as much as 15 cm (ibid). This error is

mainly caused by the motion of the satellite image within the

photocathode when the instrument tracks the satellite. However,

successful focusing of the satellite image onto the photocathode reduces

this error to the 1 cm level. Other factors such as the impulse

response of the PMT's, the amplitude of the input signals and the

background radiation also contribute to this error. These problems are

currently being solved with the replacement of the conventional PMT's

with the so-called microchannel plate photomultiplier tubes (MCP/PMT)

recently introduced on the market. These photomultipliers are

characterized by well-defined photoelectron path lengths with much

shorter transit times, and greatly reduced sensitivity in the image

position effects, the strength of the input signal and the background

radiation. When using the MCP/PMT photomultipliers, the resulting

48

I

i

I

i

I

I

I

I

I

I

!

I

i

I

I

I

B

i
I

I

I

i

I

I

I

I

i

I

I

I

I

I

I

I

I

I

errors in the inferred range appear to be below the 5 mm level (ibid).

The purpose of a discriminator is to define on the photomultiplier's

output waveform a timing point and subsequently to generate a

rectangular logic pulse that starts or stops the time interval unit. The

output waveform has a quasi-Gaussian form with a randomly varying

amplitude. This amplitude variation introduces, in the determination of

the timing point, a time bias which is highly repeatable and can be

estimated if the amplitude of the input pulse is measured and recorded

along with each observation. In practice the amplitude-dependent time

bias is determined experimentally and is compensated for by

incorporating a hardwired circuitry into the discriminator. The degree

of success in implementing this circuitry is determined experimentally

and is shown in the time walk characteristic of the discriminator. The

time walk characteristic is a curve obtained by plotting the signal

amplitude dependent time biases versus the input signal amplitudes.

The time walk characteristic shows the time bias introduced by

amplitude variations, while its RMS deviation from the zero horizontal

line characterizes the efficiency of the discriminator (ibid). For the

discrimators used during the MERIT Main Campaign, this RMS deviation

reached the value of 1.5 cm, while for discrimators currently appearing

on the market this value has been reduced to the 0.2 cm level. The

latter discriminators are currently being tested for implementation in the

continuously upgraded SLR systems.

The purpose of the Time InLerval Unit (TIU) is to measure the

round-trip flight time of the laser pulse. The rectangular logic pulse

generated by the start discriminator activates the time interval counter

while the corresponding logic pulse from the stop discriminator

commands the counter to stop. The basic component of the TIU is an

oscillator which determines the stability and accuracy of the TIU. The

oscillators used by the SLR stations which participated in the MERIT

Main Campaign were either cesium beam type or rubidium type (Schutz,

1983b). To achieve maximum accuracy, the measurement of the round-trip

flight time is split up in three parts (i.e., T = T1 + T12 - T2), where T1

is the time elapsed from the starting epoch to the first following

49
positive crossover of the oscillator, T2 is the time elapsed from the

ending epoch to the next following positive crossover, and T12 is the

time interval between the aforementioned positive crossovers. T12 is

obtained by multiplying the number of intervening positive crossovers

(N12) by the period (TO) of the master oscillator (i.e., T12 = N12 x TO).

The fractional times T1 and T2 are accurately measured either by

charging and discharging a capacitor with constant but different

currents or by using a second oscillator which is slightly off from the

master oscillator. Common biases introduced in measuring T1 and T2 are

canceled out because the times T1 and T2 are subtracted in the

computation of the round trip flight time (i.e., T). Residual errors at

the cm level are still present. These errors can be reduced at the mm

level with the use of streak-cameras employed in the Optical Time

Interval Unit (OTIU) currently being investigated for implementation in

the new two-color laser receivers (Abshire et al, 1985).

The transmitting optics are used to align the laser pulse towards

the satellite being tracked, while the receiving optics are used to

receive and focus the reflected laser pulse onto the cathode of the

photomultiplier. Unsuccessful focusing introduces the image position

effects previously mentioned. A substantially reduced single

dual-purpose telescope performs both of the above functions, thanks to

the technological advances in the field of signal detectors,

photomultipliers and discriminators. These advances introduced a

substantial reduction in the station design which in turn triggered the

construction of the highly transportable laser ranging systems TLRS-1

(Silverberg, 1982; Shelus, 1983), TLRS-2 (Transportable Laser Ranging

System No. 2), etc. These systems are extremely valuable in the study

of geophysical processes because of their ability to make observations of

limited time in remote areas, in a hostile environment.

The time receiver is either a LORAN-C or a GPS receiver and is

used to time-tag each observation in the Universal Coordinated Time

(UTC) scale. The synchronization accuracy with a LORAN-C receiver is

of the order of 1 #s, while with a GPS receiver this accuracy is of the

order of 50 ns. Synchronization errors affect the inferred optical range

5O

I

I

I

1

I

I

I

I

I

I

!

I

i

I

I

I

I

i
I

i

I

I

I

I

I

I

I

I

I

I

I

I

I

!

I

by about 4 mm/_s. The laser stations used in this study belong either

to the Goddard Laser Tracking Network (GLTN) or to the Participating

Laser Network (PLN) {Shawe et al., 1985), and their observations are

time-tagged with the UTC scale kept by USNO. According to the MERIT

standards and for reasons explained in Chapter 2, the UTC(USNO) has

been transformed to UTC(BIH).

The computer is used as an auxiliary equipment to control satellite

tracking operations, to assist the operator with such functions as data

quality and quantity assessment, maintenance, testing procedures, etc.

3.2.2 Atmospheric Channel

As the laser pulse propagates through the atmospheric channel it

experiences a continuously varying refractive index. This variation

depends primarily on the variations of the local pressure with only a

weak dependence on the local temperature and local humidity {Degnan,

1985; Abshire, 1985}. A varying refractive index, on the other hand,

bends the laser pulse according to Snell's law and also decreases the

group velocity of the laser pulse as it travels through lower pressure

layers at higher altitudes. The error due to bending of the laser pulse

is relatively small and reaches a maximum value of 3-4 cm at 10 degrees

elevation while the error due to the decrease of the group velocity is

very large, reaching the value of about 13 m at the same elevation

(Abshire, 1985). A great number of formulas have been developed to

correct the inferred optical length of the laser pulses (i.e.,the inferred

laser ranges} for atmospheric refraction effects. In the present study

and according to the MERIT standards, the Marini and Murray formula

has been used to correct for these effects {Marini and Murray, 1973).

This formula is based on the assumption of a spherically symmetric

atmospheric refraction and it uses only the pressure, temperature and

relative humidity taken at the ranging site. This formula is in error at

the 4-6 cm level as the satellite reaches an elevation angle of 20

degrees (ibid). Today, use of two-color laser ranging systems equipped

with streak-camera receivers promises an atmospheric refraction

correction with an accuracy down to the mm level, thanks to the weak

51
dependency of the group refractivity on the water vapor at optical

wavelengths {Abshire, 1985).

3.2.3 Space Segment

Although the space segment of the SLR system consists of several

satellites, the LAser GEOd ynamics Satellite ( LAGEOS ) is devoted

exclusively to geodynamic and geodetic applications (Moritz and Mueller,

1987). Lageos is a passive sphere with a 59.988 cm diameter orbiting

the earth at an altitude of about 5900 km. Its mass-to-area ratio of

1.44x103 kg/m 2 effectively minimizes the solar radiation pressure and

atmospheric drag perturbations. The high altitude of Lageos' orbit not

only reduces the effects of the poorly modeled high frequency

variations of the gravity field but also warranties good simultaneous

tracking of continental extent. The altitude of the orbit, however, is

low enough to assure the geometric strength necessary for successful

implementation of simultaneous laser satellite tracking methods.

Consequently, only laser range observations to Lageos were facilitated in

order to investigate the effectiveness of the SRD and geometric methods

in baseline determinations.

The surface of the Lageos satellite is speckled with 422

solid-cube-corner reflectors (CCR's) made of fused silica and four made

of germanium (Cohen et al., 1985). When the direction of the incoming

laser beam relative to the normal of each individual CCR reaches the

value of 25 degrees, reflection ceases to take place (Degnan, 1985). As

a result, 10 to 15 CCR's contribute to the laser pulse detected at the

receiver. Therefore, it is difficult to locate for each returning pulse its

reflection point which constitutes the ending point of the inferred

optical laser range. The location of the reflection point is needed to

compute the correction necessary to transform the ending point of the

optical laser range to the center of mass of the satellite. This

correction is referred to as center of mass correction, and its value has

been determined experimentally for different pulse widths prior to the

Lageos launch. The standard deviation in estimating this correction is

about 2 ram. In the current investigation the value of 24 cm has been

52

I

i

I

I

I

I

I

I

I

I

I

I

I

I

I

I

!

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

adopted for the center of mass correction. This value is supplied on

each data record in the tape containing the observations (see Section

3.4). Furthermore, the interference of the individual CCR returns, at

the receiver's level, may introduce a random error in the inferred laser

range the standard deviation of which reaches the value of 1.15 cm

(Fitzmaurice et al., 1977). This error is referred to as the coherent

fading effect.

3.2.4 Instrument Origin

Effective use of laser range observations necessitates a clear

identification of the starting and ending points of the inferred ranges.

As already mentioned, the ending point is identified with the center of

mass of the satellite. This is a natural choice because the equations of

motion of the satellite are conveniently expressed relative to this point.

The starting point is identified with a fixed reference point within the

laser instrument and is referred to as instrument origin. The

instrument origin usually coincides with the intersection of the

telescope's azimuth and elevation axes, but other points within the laser

instrument may be used as well. Realization of the instrument origin is

achieved either through calibration or through the common channel

receiver approach (Abshire et al., 1984; Degnan, 1985). During the

MERIT Main Campaign, calibration procedures were employed to identify

the instrument origin through the estimation of the system delay

{Schutz, 1983a}. The system delay is measured by making repeated

observations to a calibration target of known distance, usually before

and after each satellite pass. With this information, the system delay

introduced by the instrument's electronics can be readily estimated.

The distance to the calibration target is measured with a geodimeter

located very close to the laser ranging instrument. Additional surveys,

therefore, are necessary to determine the position of the geodimeter

relative to the instrument's origin. In this process, errors of the

order of 2 cm may be introduced. In order to reduce these errors,

some laser instruments are equipped with fiberoptics allowing for

self-calibration (Silverberg, 1982). The common channel-receiver

53
approach, on the other hand, eliminates the need for calibration since

the electronic system delay, except for a calibratable signal-amplitude

effect, cancels itself out (Abshire et al., 1984; Degnan, 1985). This

approach is currently being tested for implementation in the laser

ranging systems.

From the above discussion it is obvious that on the basis of the

technology employed and the models used, one could come up with a

standard deviation depicting the accuracy of the laser range

observations recorded by a certain station. This approach, however,

would not take into consideration errors resulting from improper

calibration, from operator errors or from any other errors not being

accounted for. Bearing this in mind, it was considered appropriate to

estimate, for every station used in the present study, a standard

deviation that would reflect the station's overall performance during the

MERIT Main Campaign. Such an estimate can be obtained by taking an

average value of the monthly precision estimates determined for every

station and for the entire MERIT Main Campaign by the University of

Texas (Analysis of Lageos Laser Range Data, Sept. 1983-Oct. 1984). For

the stations involved in our study, these estimates are shown in Table 2

along with the station ID's and the kind of laser instruments with which

they were equipped during the MERIT Main Campaign.

54

I

II

11

i

I

I

I

I

i

I

li

I

i

I

I

!

|

i
I

i

I

I

I

l

t

I

NAME

IQUINC2

_r_PEAK

MAZTLN

GRFI05

PLATVL

MCDON

TL0126

OTAY

QUINC3

MONPK2

HOLLAS

HUAHIN

ARELAS

ID

7109

7110

7122

7105

7112

7086

7265

7062

7886

7220

7210

7121

7907

Table 2. Station Location, Laser Instruments

and Precision Estimates.

LOCATION

Quincy, CA

Mount Laguna, CA

Mazatlan, Mexico

LASER

INSTRUMENT

MOBLAS-8

MOBLAS-4

MOBLAS-6

!OBSERVATIONAL

PRECISION (m)

0.028

0.033

0.12/0.05'

Greenbelt, MD

Platteville, CO

McDonald Obs., Ft. Davis, TX

Barstow, CA

Otay Mr., San Diego, CA

Quincy, CA

Mt. Laguna, CA

Lure Obs., Maul, HI

Huahine, Society Is., Pol.

Arequipa, Peru

MOBLAS-7

MOBLAS-2

MLRS

TLRS-I

TLRS-2

TLP_-I

TLRS-I

HOLLAS

MOBLAS-I

AREfixed

0.034

0.125

0.084

0.080

0.060

0.070

0.060

0.042

0.094

0.145

I

I

I

I

I

I

I

I

I

*before and after upgrading.

3.3 SYSTEMATIC CORRECTIONS OF THE OBSERVATIONS EXTERNAL TO

THE SLR SYSTEM

Effective application of the least-squares adjustment assumes constant

adjusted parameters over at least the time span of the observations.

Thus, baseline estimation from station coordinates requires corrected

station coordinates for their temporal variations. Besides shocks and

regional deformations, the temporal variations of station coordinates are

caused either by tectonic plate motions or by earth tides. The regional

deformations are ignored because either they are unknown or their

effects are very small. For instance, the ocean loading effects for the

stations used in this investigation are very small (Melbourne et al.,

1983). Since the time span of the observations covers about one year,

the plate tectonic motions have also been ignored. Thus we have

considered only the temporal variations of station coordinates caused by

the earth tides. The tidal corrections are accounted for by correcting

either the observations or the station coordinates. In the present study

the traditional way of correcting the station coordinates has been

55
adopted. This correction has been conveniently formulated through the

station displacement vector caused by the tidal deformation (ibid.). This

formulation is based on the same elastic response of a solid earth over

all orders within a certain degree, and if only the second degree is

considered it takes the following form:

where

GMj

GM E

: R_(-_L) • _j (3-2)

= phase lag caused by the earth's dissipative forces

= gravitational parameter of the attracting body. In the

present study only the moon (j = 2) and the sun (j = 3)

have been considered

= geocentric gravitational constant

Rj,Rj = unit vector and the magnitude of the geocentric vector

of the moon (j = 2) and the sun (j = 3) respectively

?,r = unit vector and the magnitude of the stations' geocentric

vector

h2 = nominal second-degree Love number

_a = nominal second-degree Shida number

R_,R_ = unit vector and the magnitude for the geocentric vector of

of the moon (j = 2) and the sun (j = 3) in the absence

of dissipative forces

The way the tidal displacements have been incorporated in our study

differs in two aspects from what was suggested by the MERIT

standards. The phase lag caused by the dissipative forces has been

modelled in equation (3-2), while the second-degree Love and Shida

numbers have been assumed to be frequency independent. The latter

assumption results in a maximum error of 1.3 cm in the stations' height

(ibid.). This assumption is well justified not only because the resulting

error is well below the noise level of the observations but also because

56

I

!

I

I

I

I

i

I

I

I

|

I

I

I

I

I

I

I
I

l

i

I

i

I

I

I

I

II

I

I

I

i

I

It

I

i

it reduces the bulk of the computations considerably.

3.4 DESCRIPTION OF THE DATA SET UTILIZED IN THIS INVESTIGATION

Both the SRD and the geometric methods require strict simultaneity.

Although baseline estimations based on exclusively simultaneous

observations are largely insensitive to the orbital errors (Christodoulidis

et al., 1981) and to reference frame model errors (Pavlis and Mueller,

1983), no specific campaign was ever devoted to coordinate simultaneous

tracking. This, together with the inability of the SLR systems to track

a satellite through a cloudy atmosphere, makes it even more difficult to

achieve extensive simultaneous laser tracking.

Fortunately enough, as early as 1978 the IAU Symposium No. 82 on

"Time and the Earth's Rotation" recommended setting up a working

group to organize a program of international collaboration to Monitor the

E_arth's Rotation and Intercompare the Techniques of observation and

analysis (MERIT). The proposed techniques of observation included

laser ranging and radio interferometry (Wilkins, 1980). As early as 1980

(August-October) the MERIT Short Campaign was undertaken to test and

develop the organizational arrangements that would be required for a

realistic coordination and successful implementation of the MERIT Main

Campaign which very successfully took place during the 14-month period

of September 1, 1983, to October 31, 1984.

The MERIT Working Group in collaboration with the Conventional

Terrestrial Reference System (COTES) Working Group, decided to extend

the objectives of the MERIT Main Campaign in order to include the

preparation of a catalog with a precise and consistent set of station

coordinates (Wilkins et al., 1986). Consistency of the station coordinates

is achieved by accurately linking together the reference frames realized

by each of the techniques involved. This is accomplished either

through collocations or by estimating for each technique the diurnal

differences of the earth rotation parameters. Towards achieving this

goal, it was decided to have an intensive campaign, during which, in

addition to other requirements, all the MERIT stations were asked to

57
observe as frequently as possible and in full capacity (ibid). The SLR

technique_ one of the techniques tested during the MERIT Main

Campaign9 reached its full potential. It would never have reached this

potential if it weren't for the NASA Crustal Dynamics Project whose SLR

network after several years of buildup approached during the MERIT

Main Campaign its full capacity (Coates et al., 1985). It was this peak in

the operation of the SLR system that resulted in an extensive SLR

simultaneous tracking which in turn sparked the initiation of the

present study. Consequently, the SLR observations collected during the

MERIT Main Campaign were used in this study and herein this data set

is referred to as the MERIT Main Campaign (MMC) data set.

According to the SLR organizational arrangements for project MERIT,

each of the observing stations was obligated to submit its Full Rate (FR)

observations to the Crustal Dynamics Information System (CDIS) located

at Goddard Space Flight Center (GSFC). The FR observations would

have to be submitted within three months after their collection (Schutz,

1983b). The MERIT FR data format is a character oriented format

referred to as Seasat Decimal Format (SSD) (ibid.).

The MMC data set is available to any investigator from the Crustal

Dynamics Data Bank (CDDB). This data set was sent to us upon request

in nine-track magnetic tapes. Each record in the tape is stored in SSD

format and contains the observed range, the epoch of the observation

and a number of indicators pertaining to the corrections that have been

applied and to those that have yet to be applied (i.e., atmospheric

refraction corrections, center of mass corrections, etc.). The observed

ranges stored in those records are corrected for system delay, signal

amplitude dependent effects and any other effects pertaining to the

laser instrument (see Section 3.2). The center of mass correction has a

negative sign, and therefore it should be applied to the computed range.

Each record also contains the atmospheric refraction correction computed

with the Marini and Murray formula together with the pressure,

temperature and relative humidity recorded at the observing site. The

meteorological data is included for the analysts who might prefer to

compute the atmospheric refraction correction with a different formula

58
than that of Marini and Murray. The Marini and Murray formulation was

suggested for use with the MMCdata set as an extension of the MERIT

standards (ibid.).

At the time when this study was initiated the MMC data set was not

available. Thus, at the initial stage of this work, we employed the

Lageos laser ranges recorded during the last three months of the year

1979 by stations 7114 (Owens Valley) and 7115 (Goldstone) located in

California. The initial stage of this work was primarily devoted to the

development and testing of the software necessary to edit the laser

range observations and to generate the simultaneous range and SRD

observables which constitute the input to the geometric and SRD

methods respectively. Baseline estimation, however, was based solely on

the MMC data set (see Chapter 4).

The MMC data set is not a unique data set because for each month

there have been several releases issued due to data problems or to

missing data (Section 3.6.3). Table 3 shows the monthly releases of the

MMC data set used in our study. These releases, besides having

erroneous observations, also contain data records with unacceptable

characters such as asterisks, plus and minus signs, etc. These records

should not have been there since the received data set was supposed to

have been preprocessed by the Bendix Field Engineering Corporation

(Schutz, 1985). Editing this data set was a tedious process and a

considerable amount of time was spent for this purpose.

59
Table 3.

Monthly MERIT Releases

Final Month Release Release

Sep. 1983

Oct. 1983

Nov. 1983

Dec. 1983

Jan. 1984

Feb. 1984

Mar. 1984

Apr. 1984

May 1984

Jun. 1984

Jul. 1984

Aug. 1984

Sep. 1984

Oct. 1984

C

B

E

E

E

C

C

D

D

B

B

C

B

C

yes

no

yes

yes

yes

no

no

no

no

no

no

no

11o

yes

3.5 DATA EDITING

The brief description of the error sources affecting the SLR systems

(see Section 3.2) reveals that erroneous timing of the returned laser

pulse is possible, especially when the observations are made during

daylight time with a single photon detection laser instrument (i.e., TLRS

1, 2, 3 or 4) (Shawe et al., 1985). Erroneous timing is even worse for

the laser instruments equipped with single stop Time Interval Units

(TIU). With single stop TIU's it is not possible to detect the returned

laser pulse if a noise pulse with an energy level exceeding the stop

discriminator threshold enters the receiver prior to the returned laser

pulse. Multistop TIU's, on the other hand, have the potential to reduce

erroneous detection substantially since they are designed to detect more

than one returning pulse for each of the emitted pulses. The multistop

TIU's available on the market today have time resolution at the

nanosecond level. Such low resolution makes them inadequate for use

with centimeter-level accuracy laser instruments. Erroneous laser range

observations may also result from operator errors, inadequate

maintenance and from any other sources affecting the proper operation

of the laser instrument.

6O
The erroneous observations should be detected and rejected before

the generation of the Simultaneous Range (SR) and the Simultaneous

Range Difference (SRD) observables. In common practice_ editing of the

erroneous observations is incorporated in the final adjustment of the

observations when station coordinates and baselines are estimated. In

the present study_ however_ the editing of the laser ranges should

precede the final adjustment simply because the SR and SRD observables

are obtained through an interpolation of the laser rangesD and therefore_

the presence of erroneous ranges will affect the entire set of the

generated SR and SRD observables. Thus_ effective generation of the

SR and the SRD observables requires early editing of the observed laser

ranges.

3.5.1 Data Snooping Procedure

Any kind of editing procedure requires a functional representation of

the observed ranges. The estimation of the parameters involved in this

representation allows the prediction of the observed ranges and

therefore the estimation of the errors associated with each of those

ranges. The difference of the estimated error vector (i.e._ residual

vector) from the true error vector accounts for the projection of the

true error vector into the model space generated by the column space

of A (see equation 3-12). This is the component of the true error that

is lost in the estimation process_ and therefore it cannot be recovered.

Based on the statistical properties of the estimated component of the

true error (i.e._ residual)p statistical tests may be derived to allow for

the detection and rejection of erroneous observations.

This section proceeds with the development of the statistical

formulation necessary for the derivation and implementation of those

tests. The next section contains a description and elaboration on how

the base functions employed in the functional representation of the laser

ranges should be chosen to effectively edit the laser range

observations.

Let's consider the linear adjustment model employed in the

functional representation of the observed laser ranges (Uotila_ 1986)

61
-e = AX - Lb (3-3)

where

e = true error vector of the observations

L b = vector of the observed laser ranges

X = true parameter vector

A = design matrix of the experiment

The minimum variance unbiased estimate of the parameter vector X has

the following form (ibid.)

= (ATpA)-IATpLb (3-4)

The true errors, under the null hypothesis, have a multivariate normal

distribution with 0 mean and _e variance-covariance matrix

e ~ N(O,Xe)

where

(3-5)

_e = #_P-' = variance-covariance matrix of the observations

P = weight matrix of the observations

#_ = the a priori variance of unit weight

(3-6)

The unbiased estimate of the a priori variance of unit weight takes the

following form (ibid.)

=

where

(n - u)

(AX - Lb)Tp(AX - Lb) (3-7)

n - u = degrees of freedom of the adjustment

The projection of the true error vector e into the orthogonal space

complement to the model space constitutes the estimated error (i.e.,

residual) vector

V = AX - Lb = [A(ATpA)-XATp - I]L b (3-8)

62
Substitution, in this equation, of the weight matrix P with the identity

matrix leads to the following equation

V = [A(ATA)-*A T - IlL b (3-9)

This substitution is well justified because the editing of the

observations is carried out on a station-by-station basis and for each of

the stations involved the variance-covariance matrix of their

observations is assumed to have a diagonal form and equal diagonal

elements. Thus, the weight matrix P, apart from a constant factor, is

equal to the identity matrix.

Since the true error vector e has zero mean, equation (3-3) takes

the following form

Lb = e + E(Lb) = e + AX (3-10)

where the symbol E denotes the expected value. Substitution of

equation (3-10) into (3-9) leads to the following equation

V = Me (3-11)

where

M = [A(ATA)-*A v - I] (3-12)

The matrix M is a symmetric idempotent matrix with rank (n - u), and

since MA = 0 it represents the orthogonal complement operator of the

model space (i.e., projects any vector into the orthogonal complement of

the model space) (Pope, 1976).

With the help of equation (3-11) and the law of covariance

propagation, the variance-covariance matrix of the residuals takes the

following form

Zv = H_ (3-13)

A

_v = M_ (3-14)

63
The first factor of the matrix M (eq. 3-12) depends on the design

matrix A or, in other words, on the chosen experiment. From an

estimation point of view the experiment (i.e., design matrix) should be

chosen to minimize the projection of the true observational error into

the orthogonal complement of the model space (i.e., to minimize the effect

of the true observational error on the predicted ranges}. From an

editing point of view in which we are interested, the M matrix should be

nearly diagonal. This is preferred because each individual residual (v i )

will be primarily affected by the corresponding true error ei. This

makes it easier to identify erroneous observations with one-dimensional

residual testing (see next section).

Using an independent set of analytic base functions for the

functional representation of the laser ranges results in a severe

limitation in regard to the choice of the model space. This limitation

arises because any analytic base function can be approximated by a

partial sum of monomials up to degree k. Thus, the space spanned by

the monomials (1, t, . . . , t k) closely resembles the model space spanned

by any independent set of base functions. The basis, however, selected

to span the model space will determine the conditioning of the normal

matrix (i.e., ATA) and the distribution of the approximation errors in the

interval of approximation. Choosing Chebychev polynomials as base

functions results not only in a well-conditioned normal equations matrix

but also in an even distribution of the residuals over the interval of

approximation (see next section}.

Having a functional representation for the observed laser ranges we

can proceed with the editing of those ranges using one-dimensional

statistical testing of the residuals. This is accomplished with the

one-dimensional data-snooping procedure originated by Baarda (1968).

In this procedure, under the null hypothesis the true observational

errors have a multivariate normal distribution with zero mean and _e

variance-covariance matrix (equation 3-5}. Thus, under the null

hypothesis Ho the residuals have a multivariate normal distribution

V ~ N(O, Zv) (3-15)

64
I

I

I

I

I

I

I

I

i

I

I

I

I

I

I

I

I

I

with zero mean and Y.v variance-covariance matrix (equation (3-13)).

This leads to a normal marginal distribution for each individual residual

with zero mean and variance _vi 2

vi ~ n(0,_ i) (i : 1, . , n) (3-16)

Thus, the original null hypothesis Ho is now replaced by a sequence of

null hypotheses Hot (i = 1, , n):

Ho : v i ~ n(0, _i) (i = 1, , n) (3-17)

where

v i = the residual of the ith observation

_Vt = _o

#_ = a priori variance of unit weight

mil = diagonal elements of the matrix M (equation 3-12)

(3-18)

Under the sequence of the hypothesis Ho. (i = 1, .

l

Wi = vi/_vi has a standard normal distribution

., n) the statistic

W i ~ n(O, i) (i : i, , n)

The capital letter N in the above equations denotes multivariate normal

distribution while the small letter n denotes one-dimensional normal

distribution.

The theoretical value Go in equation (3-18) is unknown and therefore

the estimated _o (i.e.,a posteriori variance of unit weight) is used to

evaluate the W i statistic. This new (Wi) statistic, under the sequence

of null hypotheses Ho. (i = I, ... n), and as the degree of freedom

1

becomes larger and larger tends to have a student's t distribution

(Pope, 1976)

65
^ v t

Wt = _ ~ tn-u, _o

where

(3-19)

_o = Type I error = P (rejecting Ho I when Ho is true)

(3-20)

n - u =: degrees of freedom

_o2 =: a posteriori variance of unit weight (equation (3-7))

The critical region C of a statistical test based on the Wi statistic takes

the following form

(3-2i)

while the Type I error becomes

.o.,[.v,.t_u,.o}

Therefore, the data-snooping procedure is carried out in the following

steps:

1.

2.

Choose a probability level (O_o = 0.005)

Take from a student t table the critical value (c) for rejection

(e.g., c = t,-u, C_o, with n - u _ 200 = 2.84)

3. Compute the individual statistic W_

4. Reject each observation Lb. leading to W1 > c

1

If the observations are correlated a slightly different data snooping

is required (Pope, 1976). The data-snooping procedure has been

effectively applied in photogrammetry (Gruen, 1979), in deformation

analysis (FIG Deformational Analysis Working Group, Heck 1982) and in

the present study to edit laser range observations.

66

I

I

I

I

I

I

I

!

1

I

I

I

I

I

I

!

!

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

3.5.2 Effectiveness of the Data Snooping Procedure in Editing the

Laser Range Observations

Implementation of the data-snooping procedure requires a functional

representation of the observed laser ranges. Before proceeding with

the choice of the base functions, we should try to single out those

properties of the base functions that would make the data-snooping

procedure more effective. To accomplish this, we should realize that the

data-snooping procedure is based on a sequence of one-dimensional

tests carried out for each residual individually. These one-dimensional

tests would be effective ff the residuals are uncorrelated and evenly

distributed over the interval of approximation. It is not possible to

obtain uncorrelated residuals with the linear adjustment model (3-3); if

this were possible it would imply that rank (M) = n, which is a

contradiction since rank (M} = n - u (see equation (3-12)). Thus we can

only look for base functions yielding residuals with reduced

correlations. This, however, is not possible because a set of k

independent analytic base functions can be uniquely mapped into the

linear space spanned by the monomials up to degree k (i.e., 1, t, . . . ,

tk). Thus we are left only with the choice of base functions that yield

an even distribution of the residuals over the interval of approximation

and presumably a well-conditioned normal equations matrix.

To understand how an uneven distribution of the residuals may

result, we choose the monomials to represent the observed laser ranges

in the interval Its, tel. Since any such arbitrary interval can be

transformed to the interval [-1, 1] by the change of variables,

2t - t E - t B tB ¢ t ( t E

(3-22)

= tE - tB -1 ( T • i

it is sufficient to examine the behavior of the monomlals (1, T, T2 , ... Tk )

in the interval [-1, 1]. In this interval each monomlal assumes the same

maximum absolute magnitude 1 at r = *1 and the same absolute minimum

magnitude 0 at r - 0. If the observed ranges would be approximated

with the polynomial

67
R.(_) =_T . T (3-23)

where

= (1, _, T 2, ... Tk) t (3-24)

= estimated monomial coefficient vector (equation (3-3))

the errors in the parameter vector X will produce small residuals for

small T (i.e.,r near zero) and large residuals for T close to 1 or close

to -1. Such uneven distribution will apparently place limits on the

effectiveness of the data-snooping procedure. In addition to this, use

of monomials as base functions gives rise to numerical problems

associated with the inversion of a nearly-singular normal matrix when k

is moderately large (Carnahan et al., 1969; Pavlis, 1982). These

numerical problems further deteriorate the effectiveness of the

data-snooping procedure. Thus, the monomials cannot be effectively

used with the data-snooping procedure.

To avoid uneven distribution of the residuals it seems reasonable to

look for functions having evenly distributed extreme values of equal

magnitude in the interval [-1, 1]. The Chebychev polynomials appear to

be good candidates since their cos(nS) origin fulfills the above

requirements. These polynomials are defined with the following

equations

T,(T) : cos(n8) ; n = 0, I, ... k (3-25)

where

e = cos-*(T) (3-26)

From this definition, one readily obtains

TO(T) : I (3-27)

T,(T) : T (3-28)

and

68

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
!

!

I

I

I

!

!

S

!

!

!

!

I

I

I

!

R

I

Tn(T) = 2T • Tn--,(T) --Tn-2(T) (3-29)

The n real roots of the polynomial T,(r) occur in the interval [-I, 1] at

the points

.,=co,{<_ _o'_'} ; i = I, ... k (3-30)

Using equations (3-25) through (3-29) it is a matter of simple exercise

to prove that the extreme values of Chebychev polynomials have the

same absolute magnitude of (1) and are evenly distributed in the

interval [-1, 1]. In this interval the Chebychev polynomials are

orthogonal with respect to a weighted integral operator with weight

function w(r):

w(.) = i /41- _.2 (3-31)

We desire, however, that these functions be orthogonal with respect to

summation as well. Fortunately, this is true but only if the summation

is carried out over a specific set of points in the interval [-1, 1]

(Pavlis, 1982)

m

Z

_=0

T,(._) Tj (_Z) =

0 i=j

m+l

i=j=O 2

m+l i:j:O

(3-32)

where r_ denote the roots of the polynomial Tm÷I(T ) given by equation

(3-30). Thus, in the formation of the normal equations matrix (ATA)

there exists a tendency for cancellations among the products of

different degree Chebychev polynomials according to the equation (3-32).

This tendency prevents the numerical problems associated with the

inversion of the normals and it justifies the term "nearly orthogonal"

often used when the Chebychev polynomials constitute the base

functions in a least-squares adjustment.

69
The effectiveness of the data-snooping procedure to edit laser range

observations is shown below for only two passes. Extensive

experimentation_ howeverj with sparse and dense data sets recorded by

many different stationsj indeed confirms that the results for these two

passes are indicative for the overall performance of the data-snooping

procedure. Fig. 1 shows the residuals_ computed with equation (3-9)_

for two passes observed by stations 7114 (Owens Valley) and 7115

(Goldstone) on October 31_ 1979. These residuals indicate that erroneous

observations with blunders as large as 350 m do exist in the original

data set. Rejection of the erroneous observations is carried out by the

data-snooping procedure (see Fig. 2). This figure shows the

distribution and the magnitude of the residuals after the application of

the data-snooping procedure. It is evident that a rejection of less than

10% of the the observations not only eliminates the blunders and makes

random the residuals but also reduces the RMS from about 46 m down to

0.11 m.

A close inspection of Fig. I and 2 also reveals that observations

having 20 m residuals not only survived in the data-snooping process

but also reduced their residuals down to the 0.10 m level. This not

only demonstrates the dependence of the residuals on the number and

the magnitude of the blunders affecting the observations but also shows

how difficult it would be to detect and reject erroneous observations by

testing the individual residuals v i alone. The data-snooping process

overcomes this difficulty by testing the normalized residual (i.e._

vl/_vl) instead. Furthermore, Fig. 2 shows that about three percent of

the residuals have an absolute magnitude greater than 0.40 m while the

rest of them are randomly distributed about zero with an RMS of 0.12 m.

Thusj it was decided to edit out the observations with residuals of

absolute magnitude greater than 0.40 m. A plot of the residuals for the

remaining observations is shown in Fig. 3. These residuals have a

random distribution about zero with an RMS ranging from 0.I0 m to 0.12

m. These RMS values are consistent with the expected accuracy of the

observing stations thereby confirming our claim that these observations

form a clean (i.e.jno presence of blunders) data set. Although the

70
IIIiI

|

4,,--

4e

0_

li

y-

(_1 :,

(_C_ ¢------"

_('r)J

Ol1.,,. | _ I

!

E$'9[| 8h';h 9S'L_ l- |9"86 I- 99"IL[ Io [L'Eh_- N NI S7UITOIS]_

'b

8

-__ua

'b

"b

9L'SI£- --

IIIIIIII

II

-sn--

O'J(.o

(Z:

(._

)-00

,--,¢_ !

SO" 9_

I

LS'S((

t_Z I

4----,--q,,m

k

.._, __.

-;,_,

,_.-_--,

6o"s s('_- Ls"s1'- s_'9_-

_OIX)N NI S7UNOIS]U.j | I

I

'b

)...

U')

't_ t---

T-

)-.

_8" 8(-Q

W

,.Qm0

:>©0

.,Q

q.1 00

"cs

°_..¢

71
p,.

II

n,-"

Or')

C) ("_,,I _ I"",-

n4

)--CO I

9h'£8

u')

(%,I

II

n- i

(I")

C_u-)

GD

cr_ OD

f_

r_

rmuj

_-nm I

G6"LG_

I

W NI

.I

-I "o • I

U')

$9"L9L _h'SO[ =-

s7unois3u

I'__w

OZ'£h[--I

SO" 96il _! "bOllH

"b

I,--.

rJ1

'1:) _,-... _ ,.,..

-$_

I,,,,-

'1=, I

6l'_I hL'6_- 19"ILl I- 19"(9Z-_

NI s7unols3u

A

=0o

_4

IIIIIIII

72
I o

•-. fJ=l'

O} _"

i_, nI

I I i I I i

I o_', o,'o .'o oo-oL o_'o- o,'o'- I,I NI s7BnoIsqH

C3

.,-,,

rr

nl,m

_" r_

lZ3

i

OZ" I 08 "IO o,,"o oo'o'- o_,'o'-" oe-oL

W NI S"IurloIs3u

"/3

4_

l...

'1:3_

"_ua

'Il,-,.-

"b

"b

,R

OZ" [ - -

'tD

,pp.-

'l_b-

_,-,-

't}

'b

0_" -o

o_..w

¢w

8

OD

¢I

"o

¢.,

v

.o

,==.s

o

_w

¢)

>

.¢

_3

o

m

'q0

m

g)
;o

r_. (Z:

u')

CO

u.J

,-. ¢_

I

O_" i 08 "iO

Pd

C_

C)

0") 'I"

P.- O:

u'_

CO

¢:_ l.¢J

•r .¢,j_

>.- OQ

I

OZ'I

J I

Oh "0

W

Og °JO

' i

Oh °0

W

i

! I I

00"0- Oh'O- 08"0-

N[ _[ S7UI'IOIS3U

"o

u_

u'_

I--

¢,_,-

-¢ ,,,

,s-

'b

't:3

OZ'[- -

'b

?4

A

c:

o

0
I

I

I

1

I

I

I

I

I

I

I

I

I

I

I

I

I

I

data-snooping procedure is very effective, it is relatively slow and

therefore very expensive since five to six iterations are needed for the

editing process to finish. Alternative ways to speed up this process

have been developed by various investigators (Pope, 1976; Gruen, 1979).

Finally, we should keep in mind that constant and time dependent

biases, no matter how large, cannot be detected by the described

procedure simply because if such biases do exist they will be absorbed

by the recovered coefficients of the Chebychev polynomials.

3.6 GENERATION OF SIMULTANEOUS RANGES AND SIMULTANEOUS RANGE

DIFFERENCES

Since Lageos is a passive satellite, it is quite unlikely for the

coobserving stations to record strictly simultaneous observations even if

the same part of the the Lageos orbit is coobserved. This happens

because for each observing site the tracking starts at a different epoch,

each laser instrument has a different repetition rate and last but not

least there will always be synchronization errors among the coobserving

sites. Implementation of the geometric and the SRD methods requires

strict simultaneity, and therefore, an interpolation of the observed laser

ranges is necessary.

Simultaneous observations for the geometric solution are obtained by

first identifying passes continuously coobserved (i.e., data gaps smaller

than 60 seconds) by four or more stations (see Sections 3.6.2 and 3.6.3).

For each of those passes the station with the least number of

observations is identified. At its observing epochs simultaneous

observations for all the remaining stations are generated through an

interpolation.

Simultaneous Range-Differences for the SRD method are obtained by

dividing the observing stations into pairs with quasi-simultaneous

observations. For each of these pairs the station with the least number

of observations is identified. Subsequently, at its observing epochs

interpolated ranges for the alternate station are generated. The SRD

observables are finally obtained by subtracting the actually observed

75
"_ li

C:3

"" ¢C_

0')

r"- n"

p,. 0

•r u'1

>- 03

,-,o

I

0_" ! GO"iO

c;

C3

3C:

r"- O::

:3'

U'I

_t U.J

'V'I

,'-*0

I

O_"_ OO°iO

Oh"JO Noo'0'- ob'o'- oe'o'-

NI _! S7UNOIS3W

Oh °iOH

i

oo'o'- o_'o- oe'o'-

NI ![ s7unoIs3U

?6

s--

ira

4_

OZ'I--

'E-

-¢=

"b

-_=

_r'.-

-¢=.,,

0"3a

'bl,--

I--.

"1=

-_=

"o

0¢'I-o

=.,

=8o

e_

"0

.p-I

GO@

oi=.¢

"0

e_

_0

._,=( O0

00

"0

s., 0

o_.=¢ @

r.,

q) U

0

oN

o_,.i m

m

i

lIlIIII

i

IIIIIIIIII
IIiIIIIIIIIIIIIIIII

It

C:)

0

f_ ,'I"

,'T-

_" U')

)-- CO " .0

02 'I 08 "iO

0h "iO X

0

I!

i_" rrnr_ _j

P-_0

i

OZ'! O0 "iO

oh'o'- o,'o'-

SIUNQ IS31=I

oh-'o oo-o'- .'o'J o_-o'- W NI II sIUnQIS3H

9'?

"b

"b

,g

.5_rl-,-

f_

_,-,-

"_ ,,,

s-

'b

{]

'b¢

Oil'--

'b

Ix)

,b_

f,tz

"b i--

_,-,-

l,--

"b

"b

O_'t-_

A

"o

.,..p

f=oU
ranges of the station with the least number of observations from the

corresponding interpolated ranges of the alternate station.

Therefore, it is essential for the success of this study to select an

interpolation method that is capable of generating laser ranges with an

accuracy compatible to that of the observations.

3.6.1 Chebychev Polynomials and Spline Functions in the Context of

Global and Piecewise Interpolation

A survey of the interpolation methods shows that these methods may be

divided into two basic categories:

- the global interpolation methods, and

- the piecewise interpolation methods.

With the global interpolation methods, a function is approximated over

the entire interval of approximation by the same linear combination of a

selected set of base functions. The function being approximated may be

known either analytically or quantitatively at a small number of base

points. The latter case depicts the situation in the present study since

the ranges to the satellite, apart from measurement errors, are known

only at each of its observing epochs. With a piecewise interpolation

method, a specific function whose values are given at a specified set of

base points is approximated by dividing the base points into successive

subsets, each of which contains two, three or more base points. Within

each subset the function is approximated by a different linear

combination of (possibly) different base functions. Boundary conditions

are imposed on the common points of adjacent subsets to make the

interpolating function continuous with (possibly) continuous firstand/or second-order derivatives over the interval of approximation.

Thus, with a piecewise interpolation method one obtains a continuous

interpolating function consisting of pieces, each of which is composed

from a different linear combination of (possibly) different base

functions.

The approximating function of a global /nterpo/at/on method, on the

other hand, generates approximate values with a relatively strong

dependence on the values the approximated function assumes at each

78

I

I

!

!

II

I

I

II

I

I

II

II

I

It

I

I

l

I
I

I

I

I

I

I

I

I

I

I

;I

I

I

I

I

I

I

I

base point. This is a desired property for our study because we would

like to obtain interpolated ranges that reflect to the highest degree of

accuracy the overall information inherent in the actually observed

ranges. Furthermore, the effect of the gaps (i.e., distances between

successive base points) should be effectively controlled, and if possible,

kept down below the noise level of the observations. In global

interpolation methods these effects of unevenly distributed gaps are

uniformly distributed over the interval of approximation provided that

the gaps are not large enough to corrupt the effectiveness of the global

interpolating function (see Section 3.6.2). With global interpolators the

generated approximate values exhibit strong correlations because the

same linear combination of base functions is used to generate

approximate values over the entire interval of approximation. Moreover,

the closer the approximated values the stronger the correlations are.

Strong correlations are not welcomed in the present study since the

generated SR and SRD observables will be considered uncorrelated in

the final adjustment when the station coordinates and baselines will be

estimated (see Sections 3.6.2, 4.3 and 4.4). The effectiveness of the

global interpolation methods is largely dependent not only on the choice

of the base functions but also on their implementation to "best"

represent the given set of data points.

The base functions most often encountered in practice are the

monomials, the Chebychev polynomials, the Fourier series, the

exponentials, etc. (Carnahan et al., 1969; Davis, 1975). Linear

combinations of either monomials or Chebychev polynomials are by far

the most important and most popular approximating functions. These

polynomials are easily operated on by addition, multiplication,

integration, differentiation, scaling and shifting and more importantly

they are closed with respect to any of these operations. Other base

functions also possess some or all of the above properties, and therefore

the polynomials would not be so important if it were not for the

Weierstrass approximation and uniform approximation theorems (Davis,

1975; pp. 24 and 107).

79

I
The choice of the base functions depends on the behavior of the

function being approximated. For instance, functions with certain

periodicities can be best approximated with Fourier base functions

sin(k_), cos(kT), k = 1, ... n, while functions with exponential behavior

are best described with exponential base functions. A combination of

Fourier series and exponentials would also be effectively used if the

behavior of the approximated function exhibits such a pattern. In the

present study, since the laser range observations are affected by

periodic and secular perturbations caused mainly by the gravity field it

is only fair to choose for their representation periodic base functions

supplemented with monomials of degree zero, one, and possibly two.

Chebychev polynomials exhibit such a behavior not only because of their

cos(ne) origin (see Section 3.5.2) but also because the Chebychev

polynomials of degree zero and one coincide with the monomials of the

same degrees. In addition to this, the optimal properties of the

Chebychev polynomials (see Section 3.5.2) make them ideal base functions

for an effective global representation of the observed laser ranges.

Having chosen a set of base functions, their linear combination

should be determined to "best" approximate the function implied by the

given set of data (i.e., function of the observed laser ranges). The best

representation of this function is determined on the basis of a chosen

criterion. Such a criterion may be chosen to either reproduce the

function at its base points or to reproduce the function and its

derivatives at a given point. We may also choose to either minimize the

maximum error of the approximation (i.e., minimax principle) or to

minimize a weighted sum of the squares of the residuals at the base

points (i.e., least squares approximation).

Reproducing the function at n base points results in an

interpolating polynomial of degree (n-l). Thus approximation of the

laser ranges with an interpolating polynomial would result in a

polynomial of degree ranging anywhere between 200 to 12000. Use of

such a high-degree polynomial is completely out of the question.

Reproducing the function and its derivatives at a given point or

minimizing the maximum error of approximation are also ruled out

8O

II

I

I

I

I

I

I

I

I

I

II

I

I

II

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

because for the former criterion the derivatives needed are not available

while for the latter the base points should coincide with the roots of

the minimax polynomial. Such a choice is not feasible in our study

simply because we have no control over when the laser ranges will be

recorded. However, if such a choice would have been possible it would

result in a minimax polynomial of very high degree, and therefore it

would have been ruled out again. Minimizing a weighted sum of the

squares of the residuals (i.e.,least squares principle) constitutes a very

good alternative not only because we can statistically select a relatively

low degree for the approximating polynomial but also because the least

squares do not reproduce the observations. The former property

prevents instability problems usually associated with the polynomial

interpolation, while the latter is desired because the observed laser

ranges are always contaminated by measurement errors. Therefore, the

functional representation of the observed laser ranges with Chebychev

polynomials whose coefficients are estimated with a least squares

adjustment, constitutes an alternative having many of the desired

properties necessary for a successful interpolation (see Section 3.6.2).

The piecewise interpolation methods generate approximate values that

are sensitive to the values of the approximated function in the

neighborhood of the interpolating point but largely insensitive to the

values of this function a little farther away from that point.

Furthermore, the effects of the gaps in the piecewise interpolation

methods are not uniformly distributed over the interval of

approximation. Thus, approximation errors committed closer to big gaps

are substantially larger than those committed closer to smaller gaps (see

Section 3.6.2). Such behavior of the approximation errors cannot be

easily controlled because their functional dependence on the magnitude

and the distribution of the gaps is not known. Direct evaluation of this

functional dependence requires laser ranges between adjacent

observations which of course are not available. The inability to

effectively control the effect of the gaps in the piecewise interpolation

methods constitutes a major drawback when these methods are compared

to the global interpolation methods (see Section 3.6.2).

81
The generated approximate values via piecewise interpolation

methods are largely uncorrelated because a different linear combination

of a (possibly different) set of base functions is used to generate the

approximate values between different pairs of successive observations.

This property differentiates the piecewise interpolation methods from the

global interpolation methods and constitutes a desired property for the

present study because the generated SR and SRD observables will be

considered uncorrelated in the final least squares adjustment when

station coordinates and baselines are estimated (see Sections 4.3 and

4.4).

With the piecewise interpolation methods a relatively small number of

base functions (usually three to five) is needed to represent the data

between adjacent data points. The choice of the base function should

be made along the same lines discussed in the global interpolation

methods. The coefficients of the different linear combinations of the

(possibly) different set of base functions for each subset of adjacent

base points are determined by imposing boundary conditions to

reproduce the functional values (i.e., observed ranges) and to obtain a

continuous interpolating function with continuous first and/or second

derivatives. One may also choose not to reproduce the observed ranges

but rather to introduce a weight function that would reflect a desired

relation between the predicted and the actually observed laser ranges.

The introduced weight function may be derived according to the

measurement errors and to the distribution of the gaps in the

neighborhood of the base points (i.e., observing epochs). Use of weight

functions with cubic splines leads to the weighted cubic splines

interpolation (Spath, 1974). Estimating the weight function is a difficult

task, and this function may not be valid for different cases. Therefore

it was considered appropriate not to introduce any weight function but

rather to reproduce the functional values (i.e., observed ranges) at each

base point (i.e., observing epoch). This would enable the evaluation of

the overall performance of the piecewise methods in regard to their

ability to generate good approximate values when relatively large gaps

exhibit an uneven distribution over the interval of approximation (see

82

I

I

I

I

i

I

I

I

I

I

I

I

I

i

I

I

!

I
I

i

I

I

I

I

|

I

I

I

I

I

I

i

I

I

I

I

Section 3.6.2). For this evaluation we made use of the easy-to-handle

cubic splines.

A cubic spline s(t) where t is in the interval [t,, t.] is defined as a

set of third-degree polynomials, each defined in the interval [tk, tk+,]

(k = I, ...n-l). These polynomials are joined at the base points (tk) SO

that the resulting cubic spline is twice differentiable at each base point.

A cubic spline defined on an interval with n base points consists of

(n-l) third-degree polynomials. Thus, its unique determination requires

4(n-l) independent sets of equations. The requirement to have a twice

differentiable cubic spline in the interval [t,, t,] introduces a set of

2(n-2) independent equations which results from the continuity

conditions required for the existence of the first and second derivatives

at the base points t2 through t,_,. The requirement to reproduce the

functional values (i.e.,observed laser ranges) at the base points (i.e.,

observing epochs) introduces another set of [2(n-2) + 2] independent

equations bringing up to (4n-6) the total number of independent

equations. The two additional equations, necessary for the unique

determination of a cubic spline, are obtained from the boundary

conditions specified for both ends of the approximation interval (i.e.,at

the base points tx and t,). Since the observed laser ranges vary

slowly within a few seconds of time, we have adopted in the present

study the following conditions:

s''(t,) = s''(t=) (3-33)

s''(t.) = s''(t.-,) (3-34)

Such an arbitrary choice influences the results very slightly in the

neighborhood of the end points (Spath, 1974; Pavlis, 1982). The

independent set of equations used to determine the (4n-4) coefficients of

a cubic spline is given in (Spath, 1974; Pavlis, 1982). The next section

contains an evaluation of the relative performance of the least squares

Chebychev and cubic spline interpolators. This evaluation is based on

their ability to effectively interpolate the observed laser ranges.

83
3.6.2 Chebychev Polynomials vs. Cubic Spline Functions in the

Functional Representation of Laser Ranges

Least-squares approximation of the laser ranges with Chebychev

polynomials requires a priori knowledge for the degree of the resulting

algebraic polynomial

Ra(T) = 60 + _lTl(r) + • • + _kTk(T) (3-35)

The degree of the algebraic polynomial Ra(T) coincides with the degree

of the kth Chebychev base function Tk(_). This degree should be

chosen to represent the data "sufficiently" in the sense that an

extension of equation (3-35) to include higher-degree Chebychev

polynomials will represent the data with the same accuracy. There is a

limit as to what degree we can go to because after a certain degree has

been reached instability problems will deteriorate the solution

substantially. With Chebychev polynomials there exists a wide range of

degrees that can be used to represent the data with the same accuracy.

This is possible because Chebychev polynomials are not seriously

affected by instability problems (see Section 3.5.2). This property will

be very useful in our study because in the generation of the SR and

SRD observables it reduces the computing time substantially (see Section

3.6.3).

The lowest-degree Chebychev polynomial that sufficiently represents

the available data is conveniently determined through statistical testing.

To construct this test it was assumed that the coefficients of the

equations (3-35) were estimated together with the weighted sum of the

squares of the residuals

R k = (VTpV)k ~=2X_--k--I (3-36)

through a least-squares adjustment. At this point we would like to

statistically test if the coefficient (_k) is significantly different from

zero. To accomplish this, we perform another adjustment supplemented

with the constraint that the coefficient (_k) is zero. If in this

adjustment the computed weighted sum of squares of the residuals is

84

I

I

I

I

I

I

|

I

I

I

I

I

i

I

I

I

I

I

I
I

i

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

Rk_i, then the weighted sum

Rc = Rk-, - Rk (3-37)

caused by the constraint (_k : 0) is distributed as _2X_ [ - _2N2(0, 1)]

independently from Rk (Hamilton, 1964). Thus, the percentage change of

the weighted sum of the residuals

Rc Rk-, - Rk X_

" ~ - Fl,n-k-Z = tn2--k--Z P = Rk •(n-k-l) = Rk (n-k-l) X2_k_ *

(3-38)

caused by the constraint (_k : 0) is distributed as tn2_k_,. Having P

and its distribution it is possible to set up a hypothesis test based on

the percentage change of the weighted sum of the residuals

Ho: _k = 0 (=: not significant change in P)

H,: _k = 0 (=: significant change in P)

The critical region of the test is obtained through the following equation

P(Fz,n-k-* >Pcr) = c__ P(Itn--k--Zl>(Pcr) %} = _ (3-39)

where _ is the significance level (i.e., Type I error) of the test.

Specifying a value for _ of 0.01 (1%} we can determine Pcr for n > 200

with the help of a t table and equation (3-39} {DeGroot, 1975)

(Pcr) % = 2.576

resulting in the critical region

C = (P _ R: P > (2.576) 2 } (3-40)

Thus the hypothesis Ho is rejected if

P > (2.576) 2 (3-41)

thereby suggesting that a higher degree is needed at the 1%

significance level. This process continues up to the degree for which

85
Ho cannot be rejected. With this test at our disposal we have the

means to determine the least squares Chebychev polynomial that

"sufficiently" represents the observed laser ranges. Therefore, we can

proceed with the relative evaluation of the cubic spline and Chebychev

interpolators. This evaluation is based on the ability of these

interpolators to "adequately" represent the laser ranges.

The relative performance of these interpolators is based on a

quantitative analysis of the orbit residuals obtained with two identical

semidynamic orbit adjustments. The SRD observables input to those

adjustment were generated through the cubic spline and the Chebychev

interpolators using the edited observations of the station pair 7114 and

7115. These observations correspond to the same two passes whose

Chebychev residuals are shown in Fig. 3. The generated cubic spline

and Chebychev SRD's were processed with the help of

- Lageos initial state vectors predicted, through numerical

integration, for the starting epochs of both passes

- the shortened version of the DE/LE200 lunar planetary ephemeris

file covering the time span of the observations, and

-the coordinates of the pole with the variations in UT1 at each

observing epoch

through two separate orbit least squares adjustments implemented by

the GEOSPP software. This software was developed by Pavlis (1982) and

during the course of this study was not only modified to comply with

the MERIT standards but also was corrected to successfully operate in

the real data environment (see Chapter 2). The TRF and the CRF frames

for each of the above adjustments are realized through the implicit

constraints discussed in Chapter 2 supplemented with the following

weighted constraints (see Section 4.1):

-coordinates of station 7114 were held fixed (i.e., ax = ay = az =

0.0001 m)

-coordinates of station 7115 were moderately weighted (i.e., _x =

ay = az = 10.0 m)

- initial state vectors were moderately weighted (i.e., ax = ay = _z

: 20 m; _ : a2 : az : 0.02m/s)

_6

I

I

I

I

I

I

I

I

I

I

I

I

I

l

I

I

I

I
I

i

I

i

I

I

I

I

I

!

I

I

I

,t|

I

i

i

The relative evaluation of the cubic spline and Chebyehev interpolators

is based on the comparison of the orbit residuals obtained with the

cubic spline and Chebyehev SRD observables. The cubic spline and

Chebychev orbit residuals for the same two passes depicted in Fig. 3

are shown in Figs. 4 and 6. The generation of the cubic spline SRD's

and the Chebychev SRD's is based on the interpolation of the station

having the larger number of observations.

A comparison of Figs. 4 and 5 shows a strong correlation between

the cubic spline orbit residuals and the data gaps of station 7115.

However, this strong correlation exists only for the pass whose

successive data gaps reach a maximum value of 60 seconds (see left

plots of Figs. 4 and 5). For the other pass, whose gaps are not larger

than 10 seconds, the correlation of the cubic spline orbit residuals with

the successive data gaps is not strongly pronounced. This happens

because the small magnitude data gaps are evenly distributed over the

entire pass (see right plots of Figs. 4 and 5). On the contrary, the

correlation of the Chebychev orbit residuals with the successive data

gaps of station 7115 is very low (compare Figs. 6 and 5). Furthermore,

it is quite clear that the cubic spline residuals are noisier than their

Chebychev counterparts. In fact, the RMS of the cubic spline residuals

for the passes shown in Figs. 4 and 6 is 0.41 m and 0.23 m respectively,

while for the same two passes the RMS of the Chebyehev orbit residuals

is 0.12 m and 0.11 m respectively. The noisy behavior and the high

RMS of the cubic spline orbit residuals is traced to the pieeewise nature

of the cubic spline interpolation. This assessment is confirmed by

computing the differences between the cubic spline and the Chebyehev

SRD's. These differences are shown in Fig. 7, and they seem to exhibit

an almost identical behavior with the cubic spline orbit residuals (see

Fig. 4), thereby confirming the anticipated fact, namely, that these

residuals are primarily caused by the cubic spline interpolator. This

result is not surprising since the piecewise nature of the cubic spline

interpolator makes it sensitive to the magnitude and distribution of the

data gaps. On the contrary, the global nature of the Chebychev

interpolator makes it insensitive to the magnitude and the distribution

_7
I

(I"1

0

Cr_ ]:

P" rr"-

'(JO

.C)

per"

_" U')

>-

O0 "10 011 "10 W

oo'oI- o_'J- oe'oINI Sq_AQIS3_

0_." [--

IIIIIIII

C}

C3 r_ ¢¢:

OD

{z:

r-_ uJoi

0_"[ 09)0

O0

u)=

O0

..QUv

nO

dS

¢0 0

IIIIIIII

88
IIIIIIII

C_

crl

O

OD

(Z:

ii

00"09 O0"OS

III

O0"Oh 00"0£ 00"0;'

a3s NI

MO0"Ot O0"Oa

Sc11:19

89

..,9,

;,,,,

¢I

_>

oe,-_

O0

O00O

r_
(:3

N

L

0;'" [ OG JO

9O
U

,.c: 0

0_

0) uu©

gl
of the data gaps provided that these gaps are not larger than 60

seconds. On the other hand, the 0.11 m and 0.12 m RMS of the

Chebychev orbit residuals is consistent with the expected acc:uracy of

stations 7114 and 7115 respectively. Thus one can safely assume that

the errors caused in the generation of the Chebychev SRD's by the

up-to-60-second data gaps are not larger than the noise level of the

observations. Extensive experimentation with dense and sparse data

sets has revealed that gaps larger than 60 seconds tend to be several

minutes long, thereby implying that the corresponding station ceased to

observe due to calibration, weather problems, etc. (see Section 3.6.3).

Thus, from now on we say that a station observes continuously if and

only if successive data gaps are not larger than 60 seconds.

Furthermore, interpolation is performed only over time intervals with

"continuous" coverage of observations.

The above discussion demonstrates that interpolation of the laser

ranges with Chebychev polynomials is superior to that of the cubic

sp]ines. Thus, for the MMC data set all the SR and SRD observables

have been generated through a Chebychev least squares approximation.

The interpolation is performed only over observing intervals with gaps

smaller than 60 seconds (see next section}.

3.6.3 Data Selection for the Generation of the Simultaneous Ranges and

the Simultaneous Range Differences

In the present study a considerable amount of time and effort was

devoted to generate the Simultaneous Range (SR) and the Simultaneous

Range Difference (SRD) observables. These two observables constitute

the input to the geometric and to the SRD methods respectively. A

geometric solution, apart from special cases, is possible if at least six

satellite positions are being coobserved by at least four stations

generally distributed in space (i.e., not lying on the same plane). If the

ground stations either form a plane (or close to forming a plane), a

geometric solution is possible only if six or more stations are involved

and each satellite position is coobserved by at least four stations (see

92

I

I

I

I

I

I

I

l

I

I

I

I

I

l

I

I

I

I
I

I

I

I

I

I

I

I

I

I

l

I

i

I

I

I

I

I

Secl.ion 2.1.3). An SRD solution, apart from special cases, is possible if

at least one pair of stations is coobserving (see Section 2.2.6).

The only data set having the potential to best fulfil] the above

reqirements is the MMC SLR data set (see Section 3.4). This data set

contains the satellite laser range observations collected during the

MERIT Main Campaign. Presently, a large amount of simultaneous laser

range observations are recorded in the WEGENER/MEDLAS project.

The MMC data set was received upon request from the CDDB data

bank on nine-track magnetic tapes, each containing a month's worth of

observations collected either by all or by some of the MERIT stations.

In each tape the information relevant to a specifi( z observation is stored

in records in the Seasat Decimal (SSD) format (see Section 3.4}. Each

record contains the observed laser range, tile epoch of the observation,

systematic corrections and indicators for the corrections that have been

applied and for those that have yet to be applied (Schutz, 19_3b). Some

observations, however, were missing from the received tapes due mainly

to data problems and delays. Thus, it was necessary at first to ide,'_tify

for each of the received tapes the ID's of the observing stations, the

number of passes per station and the number of observations per [)ass.

This information was compared against the same information published in

the monthly SLR reports issued for the entire MERIT Main Campaign by

the Center of Space Research (CSR) at the University of Texas (IJTX).

These monthly reports contained the number of passes together _^,'it._the

total number of observations recorded by each of the stations involved.

Any observations missing from the received tapes were obtained Lhrou_h

a subsequent request issued only if the utilization of those obset'vat[orls

was considered critical for the successful implementation of either" the

geometric or the SRD methods. Upon this request, new tapes identified

with the release letters A, B, C, etc. were received from the CDDB data

bank. These tapes contained the missing observations together with

those received in previous releases.

When the software necessary to process the laser range

observations was ready, it was decided not to wait any longer but

rather to use whatever releases were available up to that date {Le.,

93
Spring 1985). These releases are listed in Table 3, and they were

finally employed in our investigation. For each of those releases the

passes containing arcs coobserved by two or more stations were

identified and isolated. In the context of tile present study a pass is

defined as a satellite passage whose blind spots are shorter than 2.2

hours. With blind spots we designate the periods during which the

satellite was not observed by any of the stations involved. An arc, on

the other hand, denotes an observing period during which the

successive data gaps are shorter than 60 seconds. An arc overlap of a

Lageos pass recorded by American stations is shown in Figure 8 (left

plot).

The abscissae in these two plots denote the observing epochs

relative to the starting epoch of the pass shown on the top of these

plots. The starting epoch (84082442132) shown on Lhc. t_ _>f I!_,_: l_:f'!

plot reads 84 (1984), 08 (August), 24 (-th day), ,I (hours), 21 (minuLes)

and 32 (seconds). The ordinates on the right-hand side designate !'l::

numeric station ID's (see Table 2 in Section 3.2.4). The ordinates on the

left-hand side of these plots designate in ascending order the observing

sequence of the stations involved. Tile solid horizontal lines indb=ate

that for the period they cover the stations whose numerical ID is shown

on the right-hand side of this line have data gaps shorter than 60

seconds. The dotted horizontal lines indicate that for the period they

cover the corresponding station was not observing. The number abcJv_

the starting point of each solid line denotes the number of observations

whose successive data gaps are shorter than 60 seconds. Pig. 8

confirms the assertion made in the previous section, namely, if the gaps

are larger than 60 seconds they tend to be several minutes long. It is

evident from the left plot of Fig. 8 that all of the stations involv_,d

except 7086 have a good continuous coverage, thereby making it possible

to effectively interpolate their ranges. There is no need to interpolatr_

the ranges of station 7086 since this station has recorded the least

number of observations (see Section 3.6). The observing patDern shown

in the left plot of Fig. 8 is indicative for the performance of the

American stations during the MERIT Main Campaign. On the contrary,

94

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
i ! i i

o

o. no

III!IIIIIIIIII

0

m

..1

le

R

p.

o i

I!°

?;

95

o

8_

°glLu

mw

A

0

u

L_o
the performance of the European stations is not quite as good. This is

seen in the right plot of Fig. 8. This plot confirms the extended

overlap expected among the visibility regions of the European stations.

In spite of this, we can hardly interpolate the laser ranges for any of

the coobserving stations simply because each of the arcs involved

contains a small number of observations. Thus the question as to how

many continuous observations (i.e., observations with successive data

gaps shorter than 60 seconds) are enough for an effective interpolation

is examined next. This question will be investigated by analyzing the

errors committed in the recovery of ground truth observations. The;

recovery errors are computed by taking out a subset of observations,

referred to as ground truth observations, in such a way that the d:_t_

gaps of the observations left remain shorter than 60 seconds. The

observations left. are then used to determine a least-squares Ch_:bychev

polynomial which is subsequently employed to recover the ground truth

observations, that is, the observations not considered in its

determination. Finally, by subtracting the recovered ground truth

observations from the actually observed ones we obtain the re_¢'_v-,:.'

errors mentioned above. A sample plot of such errors for statior_ 7210

is shown in Fig. 9 together with the distribution of the ground truth

points used in the computation of those errors (i.e., left.- and

right-hand-side plots respectively).

These errors were computed for 285 ground truth observati_,ns

obtained from a total of 1427 available observations by taking out every

fourth observation. The ones left {i.e., 1142 all together) formed the

basis to determine the least-squares Chebychev polynomial whh:h was :in

turn employed to recover the ground truth observations. The abscissae

for both of the plots shown in Fig. 9 designate the epochs of the

ground truth observations relative to 0 h UT of the day showrl or, t,he

top of the figure. The ordinates of the left plot designate the recovery

errors (i.e., Range - Rint.) while the ordinates of the right plot

designate the laser ranges after they have been scaled and shifted for

plotting purposes. The RMS of these errors is 0.03 m which is exactly

equal to the expected precision of station 7210.

96

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

tit)

C_JD

ae(_

¢nt4

u_ u

¢Cu'J

0 .r'uj

,'_ u.J

aD nO:

OD=It uJ

MJ

It . t_

..... (:p¢J

_cO

mC_

O.I"L_J

_¢:_.J

_O_. LiJ

_¢:_

_J

-r

u .(-J

n->--

C_i,,rn

>.. n_ "1-

"m

/ ._

• _ _

__ j _,.e SS'S- ( i'm° _9"6- Oi'[i

tt.'O" 61"Z- u.'_" [ W')_ N| "'l_d _3H3*O]H_

mr

i j,

9_

¢

'$.

p.=

I--

•$ :_

l--

_m

Oi "0 LO'O hO'O _0'0 _0"0- kO'O" LO "0° O| "0"

( 14 ) _'INl_ - 35Ni:IEI.l'

9T

4)

m

c_

4)

"0

oo

l=

on

0

t.,

"0

0

s.,

tO

0

.o

W

._=.i

"0

"0

m

0

4)

4)

_o

o

©

I
Using a different arc observed by the same station and in the same

month we again have computed the recovery errors for" 47 ground truth

observations obtained from a total of 283 observations by taking out,

every sixth observation. These errors, together with their distributions,

are shown in Fig. 10.

The RMS of these errors is about 23% worse than the e-:pected

precision of station 7210 (Analysis of Lageos Range Data, August 1985).

Extensive experimentation with American stations has indeed confirmed

that on the basis of less than 500 observations the interpolated ranges

are affected by errors which are 15% to 20% worse than the,, noise levr-:l

of the observations. On the contrary, using more than 500 observations

the, errors caused by the interpolation hardly eve,- reach the noise level

of the observations.

The distribution of the ground truth points shown in Pigs. 9 and !0

reflects for each arc the distrib,_tion of the available observati_ms

because the ground truth observations were obtained from t}>,, availab!_,

observations by taking out either every sixth or ever'}" fou_'th

observation. With this in mind a closer inspection of PiKe,. 9 and i0

clearly reveals that small recovery errors are distributed around the

denser parts of these two arcs while for the denser arc, shown ir_ T"!:4". 0

the recovery errors are considerably smaller. Although the .sm-_!_

recovery errors are distributed around the denser parts (,f t.he,_r, !w_,

arcs, there exist in those parts recovery err'ors having an absol,_,e

magnitude of 0.10 m which is about three times larger than l,h_', c,,xpect_,d

accuracy of station 7210. Thus, the question arises as to wheth_:r the:s_ _.

large recovery errors are caused by the procedure used to c,omp,ate

them or by errors affecting the actually observed ranges.

To examine this we have plotted in Fig. 11 the Cheby(:hev r,<,_id,,mls

for both arcs shown in Figs. 9 and 10. These Chebychev residuals have

been computed by using all of the 1142 and 283 available observation,_;

for both arcs. Comparison of the recovery errors shown in Figs. 9 and

10 with the corresponding Chebychev residuals shown in Fig. 1l clearly

shows their high correlation and the equality of their RMS values as

well. Moreover, the large recovery errors seem to have their

98

I

I

!

I

I

I

i

I

I

I

I

I

!

!

I

I

I

I

I
I

I _

,

,,

I ! I I

_,l "0 O! "0 LO'O bO'O ilO'O I0"0 ° bn'O - 90"0-

I ( H ] ('£Nll:l - ":IONUI'I) •

i -!

I _ _ _

_.-.. ._g

I .'Is .-dl ,o:6 ,2:, w's -,"_ ,8'_ Lo'o_

( H i N! "'fOol ^]H3XO3H3

I

I 99

ms

o

m

Qo

v

o0

0

;>

m

,.Q

0

,._

0

0

,,0

m

°I,,4

0

0
::=e

O

• (JO

p.. _p .-.

11 li ||

t_J

CZ_ "r" n"

Q.. rr C._

¢%1

¢0

• LJ

...-. 0 _-J

0

B II II

_rl UJ

_r_ CC) L_

Q_ rr (._

GO

=:p IaJ

<X) "1-

rr"

r_ UJ rn

>. rr_ ;=

,.-.E) r.3

I

00"0

!

£0 "lO 0 O'ol- £0" 01" 90 "O_-

N[ Sq_qO [S_IEI

..=,,.

"0

-I"1

gO

¢;

L.

;>

,,,¢

e,,3

A3

0

i00
I

I

I

!

I

I

I

I

i

!

i

I

!

I

I

I

I

I

counterparts in the Chebychev residuals thereby suggesting that the

large recovery errors have originated from errors affecting the actually

observed ranges. More importantly the equality of the RMS values

between the recovery errors and the Chebychev residuals suggests that

analysis of the Chebychev residuals will give a good qualitative measure

for the accuracy of the interpolated ranges. This important result used

in the present study to control the quality of the interpolated ranges

has also been confirmed from the analysis of the orbit residuals (see

Section 4.4). Thus, we can safely state that interpolated laser ranges

based on a relatively dense data set (i.e., with more than 500

observations) will be affected on the average by approximation errors

that are smaller than the noise level of the observation. Since the

Chebychev polynomials are nearly orthogonal, deterioration due to

numerical instability ceases to exist as the number of observations

increases 11,000, 12,000 or more. Thus, since the recovery errors are

smaller for denser data sets it is preferable to interpolate dense rather

than sparse data sets provided that at least 500 observations with gaps

shorter than 60 seconds are available.

Inspection of the left and right plots of Fig. 8 reveals that the

American stations fulfillthe above requirements (i.e.,to have more than

500 observations with gaps shorter than 60 seconds) while the European

stations hardly ever fulfill such a requirement. The American station

7086, on the other hand, seems to have experienced several problems

during the pass shown in the left plot of Fig. 8. The observing

pattern, however, of this station was more or less the same for the

entire MERIT campaign thereby making its observation inappropriate for

interpolation, Fortunately, out of all the American stations used in this

investigation only two of them, station 7086 (Texas) and 7112 (Colorado),

seem to consistently have so many interruptions within all of their

observed passes. More specifically, the interruptions of station 7086 are

traced to the need of calibrating the laser within any of its observed

passes. Furthermore, for most of the passes coobserved by four or

more American stations, there exist at least three stations with dense

enough observations to be effectively interpolated. This is the

I01
necessary requirement only for the geometric solution since the required

four-station events are obtained by interpolating the observations from

three stations at the observing epochs of the fourth station for which

the actual ranges are used. For an SRD solution, the observations of

only one station are interpolated at the observing epochs of the

alternate station.

The European stations, on the other hand, seem to have experienced

several interruptions for most of their observed passes. These

interruptions have forced most of them to record relatively sparse

observations. Therefore, it is difficult to find passes having arcs

coobserved by four or more stations out of which three have dense

enough observations that can be effectively interpolated. For instance,

for the arc overlap shown in the right plot of Fig. 8 the observations

from stations 7839, 7834 and possibly from 7840 can be effectively

interpolated while the observations from any of the remaining stations

are inappropriate for interpolation. Since this is the case for most of

the passes coobserved by the European stations, it was considered

appropriate to drop these stations from their implementation in the

geometric solution. However, simultaneous observations collected in

(WEGENER/MEDLAS) project have the potential to effectively implement

the SRD method.

In other parts of the world there may exist two or three stations

coobserving. Overlap, however, for four or more stations is quite

unlikely to occur anywhere else except over North America and Europe

because only in these two parts of the world there exists a large

number of operating laser ranging stations. This and Lageos' high

altitude orbit are the main reasons why most of the passes recorded

during the MERIT Main Campaign by American stations were coobseved

by at least two of them. This resulted in a large number of American

station pairs with quasi-simultaneous observations, and therefore it was

possible to generate a large enough number of SRD observables that

have led to a steady state response of the semidynamic solutions (see

Chapter 4). Since the purpose of the present study is not to compute

all possible baselines but rather to study the performance of the SRD

102

I

i

I

I

I

I

I

i

I

!

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

!

I

I

I

I

I

I

and geometric methods in baseline determinations, it was considered

appropriate to employ the laser ranges recorded only by the American

stations for both the SRD and the geometric methods. Fig. 12 shows the

locations for those stations.

Since for each station the observing pattern in terms of data gaps

and observational density is relatively homogeneous from pass to pass,

we've shown in Fig. 13 three additional arc overlaps involving stations

that are shown in Fig. 12 but not in Fig. 8. A comparison of Figs. 8

and 13 confirms our assertion that the observing patterns of stations

7109, 7110 and 7122 which are involved in more than one pass are

homogeneous from pass to pass. Furthermore, the RMS of the

Chebychev residuals for each of the stations involved shows very small

fluctuations from arc to arc. This, however, is not true for stations

such as 7122 that were upgraded during the MERIT Main Campaign.

Table 4 lists the monthly mean RMS values of the Chebychev residuals

obtained by interpolating the edited with the data-snooping-procedure

laser range observations. The RMS values listed in Table 4 are in close.

agreement with those computed by the CSR at UTX listed in Table 2 in

Section 3.2. This close agreement confirms once again the effectivenes_

of the data-Snooping procedure to edit the laser range observations.

For each observing station the quality of the interpolated rangc, s a_,_

reflected through Lhe qualitative pattern of the Chebychev residuals is

relatively homogeneous for all of the arcs recorded during t.he MERIT

Main Campaign. Thus, in order to get a ,ee,lng for the quality of the

interpolated ranges we have plotted in Figs. 14 and 15 the Ch_._bychev

residuals for those arcs of Figs. 8 and 13 which have been marked with

an asterisk next to their numeric station ID's. If a specific station has

observed more than one of the arcs shown in Figs. 8 and 13, then the

arc whose Chebychev residuals are shown in Fig. 14 or 15 is also

marked just above its starting point with an asterisk. Table 5 lists, for

the same arcs shown in Figs. 14 and 15, the RMS values of the

Chebychev residuals together with the number of the available

observations before and after the data-snooping procedure was applied.

The condition numbers also shown in this table refer to the normal

103
"-_210

_265

_08

\ _

s (

C:,

_- _07 (

I

I

80" i

70' i

BO'

50' I

40" i

30"

20' I

iO' i

o. II

-10'

-20" i

, i

S

_ -3o' I

-qO'

-50" i

. -60' i

- __ -7o"

=_--_--_L_,_ -BO'

I

Fig. 12 Location of the American stations used in the present

study.

ORIGINAL PAGE IS

104 OF, POOR QUALIT_
I

I

I

I

I

I

!

I

I

!

i

I

I

I

I

I

I

I

2327

Bq0125519 qq

7122

5q95

.... 7110

812 _

)

261

I0

........................................... 7265 _"

IS27

.......... 7109

I I I I i J

0 0"0 0'0 0 _ 0 C)"O_ _ 00"I'f52 s OCt2 6'L_.8 s 00_ 5"_t;' 00'_ _4 0 s 00_3_36 '

8qOq07430 I0

• 169_ 7t2!_

"1 .................................... 2_2

| 832 172

t ,o/ 2

m ............................................................... 72_0

iS6 _

!

2

q

q

5

...................................................... 7112

_060

................................. ?10g

o I I I I I I

O0"O_d 00"I L_WO_ OCf2B"IO s 00_I _ 00":J o--"2_ 0 s 0 I'I0"2_ 0 IbU.__O'

Fig. 13 Arc overlap of three Lageos passes (American

stations).

105
2?

831015103029

B13 _'

235

9

......... 7062

688

537 7220 _

208B 7122"

132

o _ I

o@o"oo_ oo'b6_io_

IBIS

qq

. 7110

i : i i J

Fig. 13 (toni'd)

equations matrix that was used to estimate the coefficients of the

least-squares Chebychev polynomial. The condition numbers missing

from this table correspond to the stations whose actual rangrs w(;c(,

used to generate both the SR and SRD observables. It is evident from

this table that for some stations rejection of less than 10% of their

observations reduces the RMS of the Chebychev residuals from

unacceptable levels down to the expected accuracy of the observalJons.

Furthermore, the large condition numbers associated with stations 7210

and 7121 are caused by the strong irregularities in the distribution of

the data gaps. These strong irregularities usually occur when

malfunction of the laser instrument has resulted in concentrated

erroneous observations which in turn are rejected by the data-snooping

procedure. For instance, the observations of stations 7210 and 7121 for"

the arcs shown in Fig. 15 are affected by four large gaps and one large

gap respectively. The word large is used here in the sense that

although the data gaps are usually shorter than 60 seconds they are

still large as compared to the other gaps affecting the remaining

observations.

106

I

I

I

I

i

I

I

I

I

!

I

I

I

I

I

I

!

I
IIIIIIIIIIIIIIIIIII

._., !

.._ F-

_=_

_-_ 4

e_

o

[_I_GINAE PACg I8

D_ POOR ,_,.nr*_:, ._T"

,:;dd

d

_44

_, ¢x]

__ _-_

v:4

4d _.

r--.

.Z4

t',-i

4 _5

107

4

I

4

%-

I

i=
0

,5

II

'srrr',J UD 0

:::P t._

>-

_0

('_

II

u_

.nPd

GO C_

C_

:_' t.O

GO

r"_ t.i.J

's- U')

_0 I

LO "0 LO "0-

o

"10

¢1

¢.o

oo

oo

L"-

r'-I

¢,-

CO

¢0

_c

O_0

o_

O_

L. om

¢s

"0

*p-t m

>..:. ©

,.¢¢

0 o,.._ _ m

108

ORIGINAL PA_,_, T_,_

_OE POOR ,,-,, _L j+:,J.: ._ f _-+,-_ "

IIIIIIIIIII

IIIIIIIII
IIIIIIIII

_j 0

_ JJ

c_

rr-

_P

CD_pE

c_l._

_cn

>-- r_o

ORfGINAE PAGI_ Ig

OF POOR QUALITY,

(.o0

c_

II

(,f)

am

OuJ

>-'_0

_0

o0

109
p...

C:)

II

t.C)

• -s""

or"

=e

O0 ¢,,.j

Or_

::_ ('%j

>.- cr_

_E::)

0

u

(c)

or O

0or_

(:_ In

::7 (.C)

CO

r'7 i._

"r" ¢r_

>-

6[ "lO I

£[ "0

SO "lO

"b

,0

ii0

A

oo

v

IIIIiIIIIII!IIIIIi
IIII

C:)

II

III ,__ W NI S]_AOI53_

IIIIIIII

II | _'

IE?0

-s-"

• ,"rf,,,..

C_ O C_

c'_ L,.I

|1

I

_.a"0 _{ ?0 £070 90"0- 91 "OI H NI S7UQOIS3_

111

9_'0-0

fi

"o

"t:::)

S_'O-

_u8

rr

nu

v

¢o

L',-

_X

r-4

L'.-

_S

L',-

_S

L"-

e-4

L',-

¢o

L'--

_X

184

L',-

nO=

.o

oo

co

no

4) •=

=

4) 0

m
O

II

',x0

=I' i._E

r-I i._

"_"" (,o

>-nm

6

II

er

a_

E

rTuj

em I

$I "0!

9h'O

4I

11!

14, --

/

Ol "tO hO'O ' I0"0 i- 90 .0'- _.I .0L

14 NIS 71:I1701S31d

@

"b

C_

"b

L|'O- c_

"b@

IIi

-_@

8h "0-o

112

A

=oo

v

°_

IIIIIIIIIIIIIIIIIII
IIIIIIIIIIIIIIIIIII

U_

II

U3

,, nL_

_rd

C_r_

_0

r'_ LLJ

_--(_

O_ Q

C_

II

rrC_L_

et-

_UJ

_-(_ | !

Lh'O 6h'O-- f 1 .I , . I

II .Il_ r_T.i.r._.lIl

-_'= i- I I .... 4" fl

l• | . q. 1F _ . lm • -

. :l-r l I | --_.• II

--: .i -L'-L __" .z": .1

6Z'O

[ 61 "lO O[ "_0 OO .i0 60",,- ,i 61I "ON[ _7uno [ S:ll=l

,ira

If,

__r_

II

" "ol--

-,I= L,.I

l,--

-,Ira

"b

8_, "0--

113

ORTGINAU PA-GN _.

_:3R QU.4d_IT_

A

a,

o

v

t_

°ml
ITMO 831015

OBSEB. 543 RM5 : 0.09

6 ̧

z_

9-

.=Jo

rr_

U'J

l,L

9-

I I I I

TIME AT STATION = 7062

Fig. 15 (cont'd)

The large condition numbers, caused by the strong irregularities in

the distribution of the gaps, might deteriorate the precision of the

interpolated ranges since the number of the significant digits lost in the,

inversion of the normals is approximately equal to the base I0 logarithm

of the condition number (Forsythe et al., ]967). Even if the numerical

instability is not bad enough to affect the quality of the interpolated

ranges there might still exist large approximation errors in the

neighborhood of the large gaps. Therefore, it was decided to either

reject those arcs or to break them down into severs] subarcs with even

distribution of gaps (i.e.,small condition number) provided that for each

subarc there are enough available observations to be effectively

interpolated.

Sufficient

observations

representation of a relatively large number of

requires a high-degree Chebychev polynomial.

Determination of its coefficients through a least-squares adjustment does

not cause any numerical instability, as might have been expected,

I

I

I

I

I

t

I

I

I

I

I

I

I

I

I

I

I

114

ORIGINAL PAGE IS

_)_ iK)OR QUALITY
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

Table 5. Precision of Chebychev Interpolation

(Before and After Data Snooping).

Station

IDs

7907

7105

7086

7109

7886

7110

7122

7220

7062

7265

7112

72]0

7121

RMS

before

D.S.*(m)

26.9

7.35

0.07

O. 04

0.07

0.23

0.16

0.09

0.70

0.09

23.18

0.07

]5.88

No. of obs.

before

D.S.

648

5314

297

9611

2720

6681

2088

686

613

812

136

185

1694

RMS

after

D.S. (m)

0.13

0.02

0.06

0.02

0.07

0.02

0.15

0.09

0.09

0.09

O.11

O. 05

O. 14

No. of obs.

after

D.S.

620

5116

Condition

Number

600

561

Reference

Figures

8, 14

8, 14

296

9474

2720

6570

2072

685

543

811

107

]67

1653

- 8, 14

31 8, 14

102 8, 14

34 8, 14

48 13, 15

106 13, 15

- ]3, 15

- 13, 15

- 13, 15

4249 ]3, 15

1058 13, 15

*Data snooping

because of the near orthogonality properties of the Chebychev

polynomials (see Section 3.5.2). For instance, to sufficiently represent

9474 observations of station 7109, a 22-degree Chebychev polynomial is

required. In this case the condition number of the normals has a very

small value (i.e., 31, from Table 5), and therefore numerical instability

and ill-conditioning associated with large condition numbers (ibid.) are

not a concern in spite of the fact that a high-degree Chebycbev

polynomial was used.

All the information acquired in regard to the interpolation of the

laser ranges was implemented in the locally developed software that was

subsequently used to generate the SR and SRD observables. With this

software the passes containing arcs coobserved by two or more stations

were identified and isolated together with the starting and ending

epochs of each of the arcs involved. More specifically, for the MERIT

Main Campaign there have been identified 536 such passes observed by

the American stations shown in Fig. 12. For each of these passes the

observations recorded within an arc by all of the stations involved were

edited by using the data-snooping procedure and by taking into

115
consideration all of the necessary precautions in regard to the

magnitude, the distribution of the gaps, and the resulting numerical

instability. However, effective editing of the laser ranges requires

knowledge of the degree of the Chebychev polynomial that would

sufficiently represent the laser ranges recorded within any of the arcs

involved. This is accomplished with the help of a statistical test which

is explicitly outlined in Section 3.6.2. According to this test, any

adjustment for which the hypothesis Ho is rejected leads to another

adjustment (see Section 3.6.2). To eliminate these additional

adjustments, and therefore to reduce the computing time, it was decided

to determine a priori the degree of Chebychev polynomials that

sufficiently represent the observations of any of the stations involved

in our study (see Fig. 12). Accordingly, if the number of available

observations recorded by a specific station falls into a certain window

then for that station an appropriate degree of the Chebychev polynomial

is assigned (see Table 6). These degrees are predetermined for

specified windows and for all of the stations involved on the basis of

the actually observed ranges. Table 6 lists for prespecified windows

and for all of the stations involved the degree of the Chebychev

polynomials that sufficiently represent the ranges falling into any of

those windows.

Next the software proceeds with the identification of the passes

containing arcs that have been coobserved by four or more stations.

For each of those passes the exact overlapping times of the arcs

involved are identified and subsequently Simultaneous Ranges {SR) are

generated along the guidelines described in Section 3.6. For instance,

stations 7110, 7122 7220 and 7062 have arcs with quasi-simultaneous

observations (see Fig. 16). It is obvious from Fig. 16 that

quasi-simultaneity for all four stations occurs only in the intervals (bc)

and (hi). The epochs b and c designate the starting and ending epochs

of the first arc of station 7062, while the epochs h and i designate the

starting epoch of the second arc and the ending epoch of the fourth

arc of stations 7220 and 7062 respectively. Since station 7062 recorded

the least number of observations for both arcs bc and hi respectively,

116

I

l

I

I

I

I

I

I

I

I

I

i

I

i

I

I

l

I
I

!

I

I

I

I

i

I

I

I

|

I

I

I

I

I

I

I

No. of

_ _obs.

Station_

IDs "x_

7907

Table 6. Degree of Chebychev Polynomials.

0 - 500 -

500 1000

16

7105 12

7086 12

7109 12

7886 12

7110 12

7122 14

7220 11

7062 16

7265 II

7112 14

7210 ii

7121 14

17

16

17

17

lO00 -

2000

17

16

17

17

2000 -

3000

18

18

18

18

3000 -

5000

19

19

19

19

5000 -

9000

21

21

21

21

> 9000

22

22

22

22

13 15 17 17 21 22

16 16 19 20 21 22

19 19 20 21 21 22

ii 14 14 16 21 22

18 18 19 19 21 22

13 15 18 19 21 22

19 19 20 21 21 22

19 19 20 21 21 22

19 19 20 21 21 22

the four station events (i.e., exact simultaneous observations from four

stations) for these two arcs were generated by interpolating all the

available observations from the corresponding arcs of stations 7110, 7122

and 7220 respectively. Although in Fig. 16 we have only four-station

arc overlaps, there may exist arc overlaps including five, six or even

seven stations. In such cases all possible seven-, six-, five- or

four-station events are generated, provided that any duplication that

may occur is to be avoided.

i17
27

831015103029

9

235

613

......... 7062

597

686

722O

2088 7122

132

1615

a b c d ef g h i j

7!!]

o I ( I L _

0 0_00_00s 0 @ 8_10s 0 0"i 6_2Os 0 @ _b Cs 0 O"_2_ O_ 0 0'_40%O: O0_S<J

Fig. 16 Arc overlap of a Lageos pass involving four

American stations.

Finally all the passes with quasi-simultaneous observations from two

or more stations were employed to generate the SRD observables. For

instance, the station pair 7122-7110 shown in Fig. 16 has

quasi-simultaneous observations for the intervals (ad), (ef) and (gj)

respectively (see Fig. 16). For the first interval (ad) the SRD

observables are generated by interpolating the observations of the first

arc of station 7110 at the observing epochs of station 7122, while for

the intervals (ef) and (g j) the SRD observables are generated by

interpolating all the available observations of station 7122. Thus, having

all possible four-, five-, six- and seven-station events together with all

possible SRD's, the baselines are estimated through the geometric and

SRD methods respectively. This is the subject of the next chapter.

118

" :v<iGINAL Pt\GE IS

<.,£ POOR qUALI_

I

I

I

II

I

I

i

I

I

I

m

I

I

I

il

I

l

I

li
Chapter 4

BASELINE ESTIMATION

4.1 BASELINE ESTIMABILITY

Baseline estimability in connection with satellite geodesy is very closely

related to the concept of estimable parameters. In fact, if optimum

geometric configurations are fulfilled and enough observations are

available, the baselines can be estimated with a precision compatible to

that of the observations. This is possible only because the baselines

form a set of estimable parameters. This is true for both the range

geometric mode and the SRD semidynamic/dynamic mode methods. In

general, baselines estimated through either geometric or dynamic mode

methods are estimable only if through the appropriate adoption of the

necessary constants and units the scale is implied by the measurement

system being employed. In such cases the baselines estimated with a

dynamic mode method are referred to as best estimable parameters. The

use of this term is justified since out of all possible estimable

parameters associated with either semidynamic or dynamic mode methods,

the baselines are recovered with substantially reduced a posteriori

variances (i.e., an order of magnitude) (Van Gelder, 1978).

4.2 STEADY STATE RESPONSE OF THE GEOMETRIC AND SRD METHODS

The nature of as well as the spatial and temporal distribution of the

available observations quantify the inherently present information for

each of the estimable parameters involved (Sections 2.2, 4.3 and 4._).

The process of recovering estimable quantities from any given set. of

observations is referred to as an inversion process. This process is

effectively realized via an estimation method such as a least squares

adjustment (Sections 2.1.2 and 2.2.6). The available set of observations

119
forms the input to the inversion process, while the recovered estimable

quantities constitute the response (i.e., the output) of this process to

the given set of observations.

Different inversion processes have different responses to a specific:

set of observations and adjusted quantities. Nevertheless, as far as

satellite geodesy is concerned, these responses cannot be meaningfully

differentiated in the light of the current observational distribution and

accuracy.

The steady state response of an estimation method has been reached

if extension of its input with additional observations does not contribute

any additional information to the estimable quantities being recovered.

In the present study, steady state responses of both the geometric and

the SRD methods is sought through the extension of their input with

more and more observations. This process continues up to that point

where additional observations do not affect the recovered baselines

beyond the level of accuracy implied by the sophistication of the models

employed and by the accuracy of the available observations. With the. _

assumption that the accuracy of the baselines, recovered from laser

range observations collected during the Main MERIT Campaign, cannot

exceed the 1 cm level, the steady state response of both the SRD and

the geometric methods is said to have been reached if additional

observations do not affect the recovered baselines at the 1 am accuracy

level. Furthermore, we refer to "steady state response of the SRD

method" if the input observables are the dynamic'_ally modeled SRD

observables, while we refer to "steady state response of the geometric

method" if the input observables are the

simultaneous ranges (see Chapter 2).

In many circumstances the steady state

geometric or the SRD methods may not be

geometrically modeled

response of either the

possible because Lheir

response either diverges or" oscillates. Such a response, however, can

be reached if the geometric and/or physical characteristics of either the

geometric or the SRD methods are changed accordingly.

120

I

I

I

I

I

I

i

I

I

I

I

I

I

I

I

1

!

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

In the geometric method, for instance, near singularity cases

(Section 2.1.3) result in an extremely ill-conditioned normal equation

matrix which in turn leads to a divergent response. Divergence is

reversed to convergence if the input of the geometric method is

extended to include more and more observations collected by stations

located well away from the critical surfaces and the critical curves

(Blaha, 1971). However, if the distribution and number of available

observations is not good enough to warrant a steady state response,

then such a response can be reached via the iInp]ementation of

appropriately chosen constraints. This is accomplished either by

constraining baselines in the geometric method (Section 4.3.2) or by

increasing the lengths of continuous integration in the SRD method

(Section 4.4.2).

Increasing the length of continuous integration in the hnplcmentation

of the SRD method results in a "faster" steady state response in regard

to baseline estimation. The term "faster" is used to indicate that the

steady state response of the SRD method is reached on the basis of a

substantially reduced number of observations compared to those needed

to reach this steady state response via the short arc solutions (Section

4.4). This is the result of the geometric strength implied by the long

integration periods and manifested in the reduced order of the normal

equation matrix and in the larger values of its corresponding diagonal

elements.

Long integration periods constitute a potential source of

"instability" of the normal equation matrix. The term "instability" is

used to express the existence of high correlations among the adjusted

parameters. High correlations (greater than .99) raise the warning flag

of ill-conditioning of the normal equation matrix and a possible

divergent response of the SRD method. Employing long integration

periods independently of what the baseline pass geometry is, the state

vectors of the passes that are days away from the corresponding

adjusted initial state vector are largely insensitive to that initial state

vector. Therefore, for those passes there exists a tendency for lir_ear

dependence among the corresponding columms of the state transition

121
matrix Y*(t) (equation 2-63). This tendency further contributes to the

instability of the normal equation matrix. The instability of the normal

equation matrix gets even worse when single long baselines (longer than

about 4000 km) are estimated via the SRD method. For those baselines

there is a tendency for the coobserved passes to concentrate in the

areas between the end points of the baseline (Section 4.4.1).

With the weather-dependent laser range observations it is very

likely that the stations consitituting the end points of the estimated

baseline have coobserved only one pass falling within a specific

integration period. Under these circumstances the length of continuous

integration will be reduced to that of the duration of one single pass

which lasts about one hour for the Lageos satellite. Such short

integration lengths implemented in the same solution with integration

lengths of up to one week (Section 4.4) constitute a potential source of

instability because of the resulting inhomogeneity in the structure of

the normal equation matrix. In the design matrix, this inhomogeneity

manifests itself by the many zero entries affecting the columns

associated with the initial state vectors of those short arcs. Therefore,

care should be exercised to avoid extreme circumstances that may result

in an algorithmically singular normal equation matrix (Section 4.4.2).

This can be avoided either by employing homogeneous integration

lengths and/or by incorporating in the same SRD solution observations

from many stations. This would improve the stability of the normals

because of the strength implied by the geometry of the additional

observations. This, however, is true only if four or more stations are

involved since as it was mentioned in Section 2.1.3, at least four stations

and six targets are generally necessary for a nonsingular range

geometric solution. Therefore, in the single baseline solutions the

stability of the normals is mainly controlled by the constraint imposed

on the satellite to move along a six-parameter orbit (Chapter 2).

In the present study only single baseline solutions via the SRD

method have been performed. This was decided upon because, for

moderate baseline lengths (<3500 kin), the single SRD baseline solutions

seem to be simple, fast and very accurate.

122

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

The maximum length of continuous integration in the implementation

of the SRD method should be chosen not only to assure the steady state

response of the SRD method but also to warranty minimum propagation

of the accumulated residual orbital biases into the recovered baselines.

These biases are the ones not eliminated by the nature of the SRD

observables which reduces substantially the accumulated orbital biases

with an almost total cancellation of the accumulated radial biases (Pavlis,

1982}. Consequently, the choice of the integration lengths is limited

only to those which result in an a posteriori standard deviation of unit

weight the value of which is close to unity. Thus, in seeking the

steady state response of the SRD method we start out with short arc

solutions. If with short arc solutions the available observations are not

enough to warranty a steady state response, the integration lengths are

steadily increased up to those resulting in a steady state response. If

the a posteriori standard deviation of unit weight approaches unity

before such a response has been reached, then the corresponding

baseline is not estimated. However, any integration lengths leading to a

steady state response will be adopted even if the a posteriori standard

deviation of unit weight is smaller than one. This simply means that the

baselines via the SRD method are estimated on the basis of the minimally

required integration lengths.

The effectiveness of the a posteriori standard deviation of unit

weight to control the maximum length of continuous integation in the

single SRD baseline estimation depends heavily on the assumption of

having weighted the observations properly. This assumption has been

elaborated on and it has also been well justified in Chapter 3.

4.3 BASELINE ESTIMATION VIA THE STEADY STATE RESPONSE OF THE

GEOMETRIC METHOD

In the absence of ill-conditioning that may result from nearly critical

configurations, the steady state response of the geometric method will

be reached if a large enough number of observations is available. This

is a consequence of the fact that the weight coefficient matrix (Q_) of

the adjusted parameter vector "improves" as additional observations are

123
incorporated into the geometric solution (Blaha, 1971). The term

"improves" is used to indicate that with the same adjusted parameter

vector the matrix (Q_ - Q_a) is a positive (semi) definite matrix where

(Q}) and (Q}a) denote the weight coefficient matrices of the adjusted

parameter vector before and after the inclusion of additional

observations. Since (Q_ _ Q_a) is a positive (semi)definite matrix, it

follows that (ibid.)

Tr (O_) < Tr (Q_) (4-])

indicating that the precision of the adjusted parameters improves only if

the additional observations are drawn from the same population, thereby

making it possible to assume that the a priori variance of unit weight is

the same for both solutions and it is equal to one. Since the

mathematical model employed in the minimum constraint geometri(:

adjustments is an almost error free model, the a posteriori variance; of

unit weight tends towards one at the steady state. If these a posteriori

variances of unit weight are very close to unity so they preserve the

relative structure of the weight coefficient matrices of eq. (_1-1), the; a

posteriori variance of the adjusted parameters decreases. In the

geometric mode adjustments, the adjusted parameters are the Cartesian

coordinates of the ground stations.

In the present study we are particularly interested in estimating

baseline lengths and therefore the a posteriori variances of the

Cartesian coordinates have been mapped back into the estimated

baselines. As the a posteriori variance of the estimated Cartesian

coordinates improves with the incorporation of additional observations,

the a posteriori variance of the estimated baselines will also improve

only if the nonlinearity of the model employed allow such an

improvement to Lake place. Therefore, by including additional

observations into the geometric solutions a point, will be reachc.'d where

the estimated baseline length will not change beyond the 1 cm level.

It turns out that before reaching the steady state, th¢: a pos_(.riori

standard deviations of the estimated baselines do not reflect the change

of their length as additional observations are incorpor'ated into the

124

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

i

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

solution. For instance, it is possible for an estimated baseline, the a

posteriori standard deviation of which is 5 cm, to change its length by

as much as 76 cm with the incorporation of additional observations

(Table 8, baseline 7110-7122). This means that before reaching the

steady state response, the a posteriori standard deviations of the

estimated baselines do not reflect their accuracy but rather they

indicate how far away the solution is from its steady state. If, however,

the steady state response has been reached, the corresponding a

posteriori standard deviations assume millimeter-level values (Table 8).

Therefore, in the present study the a posteriori standard deviation of

the estimated baselines are used only as indicators showing whether or

not the steady state response has been reached. Successful utilization

of these indicators assumes close to unity a posteriori variances of unit

weight so they preserve the relative structure of the weight coefficient

matrices of the recovered parameters as additional observations are

incorporated into the solution (eq. 4-1}.

If steady state has been reached, it is furthermore assumed that the

corresponding baseline has been determined at the 1 cm accuracy level.

This is well justified only if the following three assumptions are

fulfilled: (1) the motion of the observing stations has been properly

modeled for the time span of the observations, (2) the steady state

response implies strong geometry, and {3) the accuracy of the available

observations allows recovery of the baselines at the 1 cm level (see

Section 4.2). As of the geometric modeling itself, the only errors

affecting it arise from using the three-dimensional Euclidean space

formulation rather than the four-dimensional post-Newtonian formalism

(Moyer, 1971; Bjerhammar, 1985} which of course does not affect the

estimated baselines at the 1 cm level (Moritz, 1979}.

4.3.1 Geometric Strength of the Available Observations

In the geometric approach the observed satellite positions are treated as

auxiliary independent points in space (see Section 2.1}; therefore the

strength of any minimum constraint geometric solution depends entirely

on the geometric strength implied by both the amount and the

distribution of the available observations. Thus, any meaningful

125
interpretation of the results obtained in geometric mode adjustments

should be preceded by an analysis aiming to assess the strength of the

geometry involved in each of those adjustments.

The analysis of the geometry involved in the geometric adjustments,

presented later in this section, is based on the examination of Table 7.

The first column of this table lists the ID of the stations involved, the

second, third and fourth columns contain the number of observations

per station when all the simultaneous events from Sept. 1983 to May

1984, from Sept. 1983 to Aug. 1984, and from Sept. 1983 to Oct. 1984 are

considered respectively. The bottom part of the second and third

columns contain the total number of observations, the degrees of

freedom and the a posteriori variance of unit weight obtained from

minimum constraint solutions on the basis of the data shown in the

corresponding column. The bottom part of the third column contains the

same information for two types of solutions, one minimum constraint and

one overconstraint obtained using the data listed in that column.

The geometric strength in each of the above solutions is primarily

drawn from the stations having the most observations. In assessing

this strength we assume that stations 7886 and 7220 coincide with

stations 7109 and 7110 because they are 8 m and 15 m away from each

other. With this in mind it can be easily deduced from Table 7 that 86%

of the available observations has been recorded by stations 7105, 7109,

7110 and 7122. A geometric solution, the strength of which is primarily

drawn from four stations, tends to be sensitive to how close these

stations are from their best fitting plane (Blaha, 1971) because with

stations close to forming a plane six are needed for a nonsingular space

range network (see Section 2.1.3). Thus, when six or more stations are

involved the solutions are not sensitive to the closeness of those

stations from the best fitting plane.

In the geometric solutions shown in Table 7 the additional 11% of

the available observations has been recorded by stations 7112, 7086 and

7265. Therefore, 97% of the observations has been recorded by seven

stations. The existence of these three additional stations reduces the

sensitivity of the solutions to how close stations 7105, 7109, 7110 and

126

I

i

I

i

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

Table 7

Station

ID

Global Statistics of the Geometric Adjustments

No. of Observations(*)

Sep 83 - May 84 Sep 83 - Aug 84 Sep 83 - Oct 84

7105 7143 18990 19214

7109 10198 22784 23936

7112 2884 3467 3467

7122 10996 11284 12212

7220 1969 1969 1969

7110 11549 24200 25352

7062 841 841 841

7086 1412 4400 4624

7907 176 245 245

7082 299 299 299

7210 712 756 1684

7265 4395 4395 4395

7886 m 11859 11859

N = 52574 N = 105489

DF = 14519 DF = 29478

_o 2 = 1.18 ao 2 = 1.03

N = 110097

DF = 30630

ao 2 = 1.03

DF = 30631( 2 )

_o 2 = 1.05(2)

(*) Minimum constraint solution, coordinates fixed(3): X,Y,Z for

7109; X,Y for 7122; Z for 7105.

(2) Overconstraint solution, coordinates fixed(3): X,Y,Z for 7109;

X,Y,Z for 7122; Z for 7105.

(3) Coordinates fixed to those of (CSR)85LOI (Section 4.4.2)

N = total number of observations

DF = degrees of freedom

ao 2 = a posteriori variance of unit weight

7122 are from their best fitting plane. The remaining four stations,

7062, 7907, 7082 and 7210, have contributed only the last 3% of the

available observations. Therefore, their contribution to the geometric

strength is very minor as compared to that of the previous stations.

Next we assume that stations 7112, 7086 and 7265, which contribute

the 11% of the available observations, make the geometric adjustments

insensitive to how close stations 7105, 7109, 7110 and 7122 are from their

best fitting plane. Even if this were true, the strength of the

geometric solutions tends to still be weak because the stations that have

127
Table 8 Baseline Steady State Responseof the Geometric Solution

No. of Correlation Data Solution

Baseline Observ. _ 0.80 Length (m) Set I Type 2

7109-7110 9,186 None 883601.637 ± 0.02 A 1 I

21,772 " .608 ± 0.02 B 1

22,924 " .661 ± 0.02 C 1

" " 883601.661 ± 0.02 C 2

" " 883602.245 ± 0.009 C 3

7109-7265 3,363 " 627043.412 ± 0.02 A 1

" " .452 ± 0.01 B 1

" " .535 ± 0.01 C 1

" " .535 ± 0.01 C 2

" " .988 ± 0.005 C 3

7109-7886 11,859 None 7.746 * 0.002 B 1

" " .746 " C 1

" " .746 " C 2

" " .746 " C 3

7109-7122 8,644 None 2280712.335 ± 0.07 A ]

8,932 " 2.700 ± 0.05 B 1

9,860 " 3.188 ± 0.05 C l

" " 3.188 ± 0.05 C 2

" " 4.949 ± 0.0005 C 3

7110-7122 10,060 PZFZs=0.998 1437137.428 ± 0.05 A 1

10,348 " .780 * 0.04 B 1

11,276 " 8.187 * 0.03 C 1

" PXFXs=0.996 8.187 ± 0.03 C 2

" None 9.288 ± 0.009 C 3

7110-7220 1,576 PXFXs=0.801 15.225 ± 0.006 A ]

PYFYS=0.990

PZFZs=0.998

" PYFYs=0.987 .221 ± 0.005 B 1

PZFZs=0.998

" PYFYs=0.981 .218 ± 0.005 C 1

PZFZs=0.997

" PXFXs=0.997 .218 ± 0.005 C 2

PYFYs=0.975

PZFZs=0.800

" PYFYs=0.965 .208 ± 0.005 C 3

7110-7265 3,866 PYFYs=0.9] 274069.453 ± 0.01 A 1

PZFZs=0.996

" PYFYs=0.887 .383 ± 0.008 B 1

PZFZs=0.995

" PYFYs=0.850 .355 ± 0.008 C 1

PZFZs=O.994

" PXFXs=0.994 .355 ± 0.008 C 2

" None .474 ± 0.007 C 3

128

I

I

I

I

I

I

I

I

I

I

I

I

i

I

I
Table 8 (cont' d)

Baseline No. of Correlation Data Solution

Observ. _ 0.80 Length (m) Set _ Type =

7110-7886 11,859 None 883605.698 ± 0.02 B 1

" " .751 ± 0.02 C 1

" " .751 ± 0.02 C 2

" " 6.335 * 0.009 C 3

4,184 PZFZs=0.995 1663980.848 * 0.05 A 1

" =0. 994 i.161 ± O.04 B 1

" =0.993 1.555 ± O.04 C 1

" PXFXs=O. 993 1.555 " C 2

" None 2.823 ± 0.005 C 3

0 None 2280718.021 ± 0.05 B 1

0 " .509 ± 0.05 C 1

0 " 18.509 ± 0.05 C 2

0 " 20. 269 * 0.002 C 3

0 274066.158 * 0.010 A 1

I

I

I

I

I

I

I

I

I

7122-7265

7122-7886

7220-7265

7265-7886

PYFYs=0.90

PZFZs=0.994

0 PYFYS=0.874 .090 * 0.009 B 1

PZFZs=0.993

0 PYFYS=0.833 .064 ± 0.008 C 1

PZFZs=0.991

0 PXFXs=0.991 .064 ± 0.008 C 2

0 None .189 ± 0.007 C 3

0 None 627048.351 * 0.01 B 1

0 " .434 ± 0.01 C 1

0 " .434 ± 0.01 C 2

0 " .887 ± 0.006 C 3

I

I

I

I

I

I

I

' Data Sets: A Sep 83 - May 84

B Sep 83 - Aug 84

C Sep 83 - Oct 84

= Solution Type:

1 Minimum constraint solution, Cartesian coordinates fixed:

X,Y,Z for 7109; X,Y for 7122; Z for 7105

2 Hinimmn constraint solution, Cartesian coordinates fixed:

X,Y,Z for 7109; Y,Z for 7122; Z for 7105

3 Overconstraint solution, Cartesian coordinates fixed:

X,Y,Z for 7109 an__dd7122; Z for 7105

recorded the 91% of the observations (i.e., 7105, 7109, 7122, 7220, 7110,

7062, 7265 and 7886) are concentrated around two intersecting lines

defined by station 7109 with 7122 and 7122 with 7105 (see Fig. 8). Since

129
two intersecting lines belong in the family of second-order curves, it is

reasonable to expect that near singularity B (Section 2.1.3) tends to

weaken the strength of those solutions. If we furthermore assume that

the 9% of the observations recorded by stations 7082, 7086, 7112, 7907

and 7210 which are located well away from these two intersecting lines,

make the geometric solutions insensitive to the previously mentioned

near singularity B, it is still possible that near singlarity C might be

present. This is a consequence of employing a network of relatively

large extent. With such a network the simultaneously observed satellite

positions tend to concentrate on the area extended above the middle

part of the network and therefore to be closer to a plane. This in turn

would lead to near singularity C because off-plane targets are needed to

avoid this type of singularity (Section 2.1.3).

It is evident from the above discussion that the geometric strength

implied by the available observations is relatively weak. Therefore the

minimum constraint geometric solutions will be strongly influenced by

the distribution of the available observations thereby making it difficult

to recover the relative geometry of the ground stations (next section).

4.3.2 Baseline Results

On the basis of the data listed in the last three columns of Table 7, we

have performed five least squares geometric mode adjustments. Four of

these adjustments are based on minimum constraints, and the fifth is

based on constraining one additional Cartesian coordinate than those

required for a minimum constraint solution. More specifically, the first

three minimum constraint geometric adjustments have been obtained by

using all the events from Sept. 1983 to May 1984, from Sept. 1983 to

Aug. 1984 and from Sept. 1983 to Oct. 1984. An event occurs when four

or more stations observe the same satellite position. In these three

adjustments, the same mimimum constraints have been used as they are

implied by fixing the X,Y,Z Cartesian coordinates of station 7109, the X

and Y coordinates of station 7122, and the Z coordinate of station 7105.

The fourth adjustment was obtained by using all the available events

and by fixing the Y and Z rather than X and Y coordinates of station

130

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

7122. This adjustment was necessary to confirm further the weakness

of the geometry that seriously affects the minimum constraint geometric

solutions. Finally, the overconstraint solution performed on the basis of

the data listed in the fourth column of Table 7 was necessary to reach

the steady state response for longer baselines. The results for all of

the adjustments decribed above are listed in Tables 7 and 8. Table 8

contains for each baseline and for all of the solutions performed the

number of observations per baseline, the correlations that are greater

than or equal to 0.80 and the estimated baseline lengths followed by

their a posteriori standard deviations. Out of all of the estimated

baselines, only those are listed for which the a posteriori standard

deviations in the overconstraint solution have reached the millimeter

level. The symbol PzFzs, in Table 8, designates the correlation between

the Z coordinate of the first {F) station and the Z coordinate of the

second {S} station of the corresponding baseline.

A close inspection of the minimum constraint solutions clearly

reveals that the estimated baselines whose a posteriori standard

deviation is at the centimeter level change their length by as much as

70 cm to 80 cm when additional observations are incorporated into the

solutions (baselines 7110-7122 and 7122-7265). The number of additional

observations incorporated in each solution as compared to the previous

one is easily deduced from the information given in Table 7 for each one

of those solutions.

The baselines whose length is smaller than 1000 km change their

length to within -10 cm to 12 cm with a tendency of positive increase.

This positive increase, which is clearly pronounced for the longer

baselines (7110-7122 and 7122-7265) takes place toward the correct

length of the corresponding baselines. This is easily seen by comparing

these lengths with the ones obtained in the overconstraint adjustment

(i.e., Solution C 3). The fact that the lengths of the baselines, which

have been estimated at millimeter precision level via the overconstraint

adjustment are very close to their true lengths is elaborated on at the

end of this section, and it has also been confirmed by comparing their

lengths with those obtained via the SRD method as well as with those

131
computed by other computational centers such as CSR and ZIPE via

dynamic long-arc solutions {see next section}.

The geometry implied by the number and the distribution of the

available observations manifests itself in the minimum constraint

solutions shown in Table 8. This geometry is not strong enough to

warrant a steady state response since all the baselines but 7110-7220

and 7109-7886 change their length by several centimeters with the

incorporation of additional observations. The existence of weak

geometry is the result of the expected near singularities described in

the previous section. This weak geometry is also confirmed by the high

correlations prevailing among the station coordinates that were not

constrained in the implementation of the minimum constraint solutions.

This weakness is further confirmed by using all the available events

from Sept. 1983 to Oct. 1984 and by changing the minimum constraints

from (1) to (2} {Table 8). In solution (2) the correlations prevail among

the X coordinates of the stations for which the correlations in solution

(1) were high among their Z coordinates. It is interesting to note that

the correlations among the X coordinates in solution (2) are almost the

same as those among the corresponding Z coordinates of solution (1).

This is the result of fixing the X and Y coordinates of station 7122 in

solution (1) and the Y and Z coordinates of the same station in solution

{2). The reduction of the correlations among the Y coordinates of the

corresponding stations from (1) to (2} simply reveals that the minimum

constraints (2) have better stability characteristics than those of

minimum constraint (1).

The correlations among the coordinates of station 7110 and 7220 are

also high because these stations are very close to each other, about 15

m apart. This, in light of weak geometry, makes their separation

difficult thereby bringing the corresponding correlations to high levels.

The correlations among the coordinates of station 7886 and the

coordinates of all of the other stations involved are small because 7886

is only about 8 m away from 7109 which was held fixed in the

implementation of both solutions (1) and {2}. The estimated baseline

lengths and their a posteriori standard deviations remained of course

132

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

!

!

I

the same for both solutions because the baselines are invariant

Quantities under any minimum consraints and as such their estimated

lengths together with their variances should not change from one

minimum constraint solution to the next.

The above discussion evidently reveals that the geometry implied by

all of the simultaneously observed satellite positions is weak thereby

making the minimum constraint solutions susceptible to the distribution

of the available observations. For reasons mentioned in the previous

section, the simultaneously observed satellite positions tend to

concentrate on the area extended above the middle part of the employed

network (see Fig. 12). This observational coverage together with the

weak geometry will imply, via the minimum constraint solutionsp a range

space network with a tendency to shrink towards its center and more

specifically towards the area where most of the available observations

are concentrated. This simply means that the scale of the recovered

range space network is not properly implied by the geometric strength

of the available observations. This fact is also confirmed by the

positive increase of the longer baselines when additional observations

are incorporated into the minimum constraint solutions. Implementation

of the scale in the geometric solution has been attempted by fixing the

third coordinate of station 7122 (solution (3)) in addition to the

coordinates fixed in the implementation of the minimum constraint

solutions. By fixing the third coordinate of this station we implicitly fix

the length of baseline 7109-7122.

As mentioned earlier in this section, Table 8 contains the results

only for those baselines for which the a posteriori standard deviations

in the overconstrained solution are smaller than 1 cm. Although all of

the remaining baselines have been estimated in all solutions shown in

this table, they are not included there because the steady state

response for those baselines is assumed not to have been reached. The

reason for claiming this is based on the fact that baselines 7110-7220

and 7109-7886 are the only ones that do not change their lengths at the

centimeter level either when additional observations are incorporated

into the solution or when additional constraints are applied.

133
Furthermore, only these two baselines have been estimated in the

minimum constraint solutions with an a posteriori standard deviation at

the milimeter level thereby suggesting that a steady state response in

the geometric mode environment assumes millimeter level a posteriori

standard deviations of the estimated baselines.

Application of an additional constraint brings the a posteriori

standard deviation of all the remaining baselines shown in Table 8 to

the millimeter level thereby implying that for those baselines the steady

state response has been reached via the implementation of the scale in

the resulting range space network.

This was expected because in the geometric solutions the baselines

constitute a set of estimable quantities and as such constraining one of

them (i.e., 7109-7122) will lead, apart from nonlinear terms, to reduced

unscaled standard deviations for the remaining baselines. Their scaled

standard deviations (i.e., a posteriori standard deviations) will also be

reduced if the scaling factors (i.e., the a posteriori standard deviations

of unit weight) are close enough to unity so they preserve the relative

structure of the corresponding unscaled standard deviations. In fact

the a posteriori standard deviation of unit weight in the overconstrained

solutions is 1.05, just slightly larger than that of the corresponding

minimum constraint which is 1.03, making it possible to preserve in

those solutions the same relative structure of the scaled variancecovariance matrices as that of their corresponding weight coefficient

matrices (Table 8). This simply means that the additional constraint did

not distort the geometry but rather made it stronger, and therefore it

was possible to reach through the solution (3) steady state response for

longer baselines.

Since the geometric solutions are not affected by the erroneous

modeling of either the satellite motion or the motion of reference frames,

we have taken solution (3) to constitute the standards of comparison in

assessing the effectiveness of the SRD method versus the dynamic mode

methods. Solution (3) was chosen over the corresponding minimum

constrained solutions because, of all of the available events, the steady

state response for more than the two very short baselines {i.e.,

134

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

l

I

7110-7220 and 7109-7886) was not possible via the minimum constrained

solutions (Table 8).

The only errors affecting solution (3) result from the ineffectiveness

to implement the ill-defined scale by fixing the third coordinate of

station 7122 to that of the (CSR)85L01 solution (see Section 4.4.2). This,

however, is difficult to assess, but on the basis of the results presented

in the next section it can be safely stated that the scale along a band

in the direction of baseline 7109-7122 has been properly defined by

constraining this baseline. In fact the value of the scale factor varies

from solution C2 to that of C3, from 6.61x10 -7 to 7.72x10 -7 for baselines

longer than 627 km while for the two baselines of about 274 km the

value of the scale factor drops down to 4.5x10 -7. This dependency of

the scale factor on the baseline lengths is expected since for shorter

baselines the implied geometry is stronger.

4.4 BASELINE ESTIMATION VIA THE STEADY STATE RESPONSE OF THE

SRD METHOD

The incorporation of additional observations into the SRD solutions will

lead to a steady state response for the same reasons mentioned in the

previous section. This response is even faster if observations of

"improved quality,, are available. The term improved quality indicates

improvement of the weight coefficient matrix. For instance, from two

sets of observations (a) and (b) having identical sizes, the sc_'t {b) is of

improved quality if the difference {Pb - Pa) of their weight matrices is

a positive (semi)definite matrix. Based on this, it is a trivial exercise to

prove that the matrix (Q_a _ Q_b) is positive (semi)definite if the design

matrix is the same for both sets of observations (a) and (b) and of full

rank. Q_a and Q_b denote the weight coefficient matrices of the same

adjusted parameter vector obtained on the basis of sets (a) and {b).

Since (Q_a _ Q_b) is a positive (semi)definite matrix, it follows that

Wr(Q_ b) _ Tr(Q_ a) (4-2)

Although the two sets (a) and (b) are drawn from different populations,

consistency between the accuracy of the models employed and the

135
accuracy of the observations leads to a close-to-unity a posteriori

variance of unit weight, thereby reducing, as it follows from equation

(4-2}, the a posteriori variances of the adjusted parameter vector. This

vector contains the earth-fixed Cartesian coordinates of the ground

stations and the inertial initial state vectors which are treated in the

present study as nuisance parameters. The adjusted station coordinates

and their statistics are used to estimate the baseline lengths and their a

posteriori standard deviations. Reduced variances in station

coordinates, obtained from improved quality observations, yield baselines

the variances of which are also reduced only if the non-linearity of the

employed models allows such a reduction to take place. These

reductions, taking place for most of the SRD solutions presented in the

next section, imply a faster steady state response.

Constraining of estimable quantities in any SRD solution results in

baselines with reduced variances thereby implying again a faster steady

state response. Thus, incorporation of additional observations, improved

quality observations, and constraining of estimable quantities constitute

the three major factors leading to a steady state response. In the SRD

solutions, this response is achieved by balancing the contribution of the

first and the third factors because we have no control over the second

factor since the quality of already recorded observations does not

change. Furthermore, the goal has been to restrict as much as possible

the contribution of the third factor because the steady state response

achieved on the basis of constraining estimable quantities may be

affected severely by the errors affecting those quantities. These

quantities are in error because their recovery has also been based on

erroneous observations. If the propagation of the errors affecting the

constrained estimable quantities is not well controlled, it will lead to

erroneous recovery of the adjusted parameter vector and therefore to

erroneous baseline estimates. Even worse, the a posteriori variances of

those estimates might also be reduced making it more difficult to

reliably assess their accuracy.

136

I

l

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

l

I

4.4.1 Paes Baseline Geometry and Its Manifestation in the Design Matrix

of the SRD Observable

The underlying characteristic inherently present in any factor leading

to a steady state response is its ability to strengthen the geometry

implied in the SRD solution under question. Thus, with stronger

geometry, the steady state response will be achieved on the basis of a

reduced number of observations. In other words, observations

containing more information for the recovered estimable quantities will

lead to a faster steady state response.

Fig. 17 shows the horizontal plane of the starting point (I) of the

baseline (12). The ending point (2) of the baseline and the

simultaneously observed satellite positions have been projected on this

plane. This simplification is aimed at revealing the information the SRD

observables contain about the estimated baseline for two characteristic

geometric configurations, namely, one when the subsatellite tracks are

parallel to the estimated baseline (Fig. 17a) and one when the

subsatellite tracks are perpendicular to that baseline (Fig. 17b and c).

It is a matter of trivial trigonometric manipulations to confirm that

with subsatellite tracks parallel to the estimated baseline the magnitude

of the SRD observable tends towards the length of the estimated

baseline as the subsatellite point moves away from either end of the

baseline along the subsatellite track denoted by (5) in Fig. 17a. Thus,

with this geometry the SRD observable directly relates to the length of

the baseline. However, as the subsatellite point moves towards the

perpendicular bisector of baseline (12"), the SRD observable tends

towards zero.

With the subsateUite tracks perpendicular to the estimated baseline,

the SRD observable tends towards zero as the subsatellite point moves

away from the projected baseline along the line {z) in Fig. 17b and c.

This observable also tends towards zero as the subsatellite track (_)

tends toward the perpendicular bisector of the projected baseline (12").

However, the SRD observable tends towards the baseline length as the

subsatellite point moves toward point (I) and as this point moves in the

direction of the baseline and away from it (Fig. 17c). Thus, for a faster

137
2 !

I

I

I

I

j' I. _ |

I

I

I

I

I

I

(a) (b) (c)

z subsatellite track

j" subsatellite position

12" projected baseline length

p," projected observed range from 1

pa" projected observed range from 2

6p projected SRD observable

Fig. 17 Subsatellite track-baseline geometry.

138
l

i

I

I

I

I

steady state response, parallel passes, the coobserved parts of which

are extended well away from the endpoints of the baseline, should be

preferred. When perpendicular passes are available, the ones located

well away from the perpendicular bisector plane of the estimated

baseline should be preferred.

The previously described favorable geometry manifests itself in the

strong independence prevailing among the columns of the resulting

design matrix. The term "strong independence" indicates the sharp

variation characterizing the ratio of the corresponding entries between

any two columns of the design matrix. The sharp variation resulting

from favorable geometry is seen in the partial derivatives of the SRD

observable taken with respect to the coordinates of one of the ground

stations and to those of the inertial initial state vector (Appendix A)

in

a_2 pj 2

aXo Pj 2 Pj, axo

I

I

I

I

I

I

I

I

<4-3)

(4-4)

where X j, Yl, Y2 are the earth-fixed position vectors of the satellite at

epoch j and stations (1) and (2) respectively.

Since the partial derivatives of equation (4-3) have a common

denominator, the variation of their ratios is controlled by the variations

prevailing among the coordinates of the satellite as it moves along the

coobserved part of the pass referred to from now on as "common part."

The longer the common parts, the sharper the variations are for the

ratios of the partial derivatives of equation (4-3). Common parts tend

to be longer as the baseline decreases thereby resulting in sharper

variations among the ratio of those partial derivatives for shorter

baselines. However, these variations become sharper for both short and

long baselines as more and more passes of different spatial distribution

are incorporated into the solution. For long baselines (i.e.,> 3500 kin),

the common part tends to be short if the pass is parallel to the baseline

thereby reducing the variation of the ratio of those partial derivatives.

139
For perpendicular passes, the common parts tend to be a little longer

for long baselines, allowing therefore for sharper variations in the ratio

of the partial derivatives of equation (4-3}.

The variations of the ratios between the partial derivatives of

equation (4-4) as well as between the partial derivatives of equations

{4-3) and (4-4} tend for perpendicular passes to be substantially

reduced. The variations of the ratios among the partial derivatives of

equation (4-4) are controlled by the changes in the observed ranges p j2

and P j l, the changes in the satellite coordinates and the changes in the

columns of the state transition matrix. When the SRD observables tend

to become shorter (Fig. 17b and c}, the observed ranges tend to be

equal a,d therefore the coordinates of the coobserved satellite positions

tend to cancel out. With such cancellations taking place and with the

tendency to have equal denominators in those partial derivatives, the

variations of their ratios are primarily controlled by those among the

transformed state transition matrix ((SNP)'aRj/aXo, ...), the columns of

which, apart from the common denominator, tend in this case to be

multiplied by the components of the estimated baseline. Thus, for

shorter baselines the variations among these partial derivatives tend to

be sharper than those for longer baselines (Appendix B). The

unfavorable geometry associated with almost perpendicular passes

manifests itself with an almost constant numerator in the partial

derivative taken with respect to the station coordinate measured along

the direction of the estimated baseline. This happens because in this

direction the satellite coordinates change very little as the satellite

moves along a pass almost perpendicular to the direction of the baseline.

This geometry also leads to small variations of the numerators in the

partial derivatives of equation (4-4); therefore, in the design matrix the

ratio between the entries of the columns corresponding to the station

coordinate, recorded along the baseline direction, and to the coordinates

of the initial state vectors will exhibit very little variation because the

denominators of the corresponding entries are nearly equal.

Reduced relative variations among the columns of the design matrix

lead to high correlations among the corresponding adjusted parameters.

140

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

l

l

I

I

I

I

i

I

I

I

I

i

l

I

i

I

I

These correlations tend to increase or decrease as the homogeneity and

the strength of the implied geometry decreases or increases

respectively, a fact also manifested in the resulting poor or strong

conditioning of the normal matrix.

The above discussion, the aim of which is to help in understanding

the results presented in the next section, evidently reveals that in

single SRD baseline solutions, the strength of the geometry fades away

as the length of the estimated baseline increases. For longer baselines

(i.e., > 4000 km), the steady state response cannot be reached because

poor geometry results in an ill-conditioned normal equation matrix.

The weak geometry associated with long single baseline solutions is

substantially improved through a network solution. Such a solution has

not been performed in the present study because of the large number

of observations and the limiting computer capabilities. If SRD normal

points were available a network solution would have been possible

which, however, is beyond the scope of this study.

4.4.2 Baseline Results

This section describes the process of achieving steady state response

for those baselines only for which such a response, through single SRD

baseline solutions, was possible on the basis of the data collected during

the main MERIT campaign.

The adjusted parameter vector in single SRD baseline solutions

contains the coordinates of the baseline end-points and the components

of the initial state vectors for all of the arcs involved. In these

solutions, obtained via a weighted least squares adjustment, the

coordinates of one baseline end are held fixed while the coordinates of

the other end and the components of the initial state vectors are

allowed to adjust by assigning, through their weight matrix, the

following standard deviations

#x = #U = #z = 0.0001 m

o"x = crg = o"z = 20 m

O'x o = o'g ° = O'Zo = 50 m

_Xo = _Yo = _Zo = 5 cm/s

(fixed baseline end)

("free" baseline end)

(initial position)

(initial velocity)

141
The term "free" designates the adjusting baseline end.

The standard deviation of 20 m assigned on the coordinates of the

free baseline end reflects an accuracy estimate of their approximate

values. These values have been obtained by altering from 20 m to

100 m the coordinates of the (CSR)S5L01 solution (Tapley et al., 1985b).

The standard deviations of 50 m and 5 cm/s assigned to the approximate

values of the components of the initial state vectors, reflect the

accuracy of those values. They have been obtained by integrating, over

a period of a month, Lageos' equations of motion using as initial state

vectors those computed at EG&G for the entire main MERIT campaign

(Pavlis, 1986; private communication). At the beginning of each month

the EG&G estimates differ from those obtained by integrating Lageos'

equation of motion over the period of the previous month because a

simplified orbital modeling was used in the integration. These

differences ranging from 50 m to 150 m for the position and from 5 cm/s

to 13 cm/s for the velocity reflect the assigned accuracies to the

components of the initial state vector.

Constraining one baseline end, by applying large weights to its

coordinates, results in implicit constraining of its latitude and geocentric

distance which both constitute, in the dynamic/semidynamic environment,

a set of estimable quantities if the scale has been defined through the

appropriate adoption of units and constants (i.e., units of time, velocity

of light, etc.) and if a solution is possible through appropriately

adopting physical constraints (Van Gelder, 1978}. Weighting moderately

the components of the initial state vectors results also in moderate

implicit weighting of the estimable quantities associated with the

geometric characteristics of the orbit. Adjusting only one baseline end

and the initial state vectors, in either dynamic or semidynamic single

baseline solutions, results in implicit and/or explicit constraining of a

large number of estimable quantities. For instance constraining the Cnm

and Snm potential coefficients together with the geocentric gravitational

constant (GME) and the mean radius of the earth (aE) results in explicit

constraining of the Anm and Bnm coefficients which in the dynamic/

semidynamic environment constitute a set of estimable quantities {ibid.)

142

I

I

i

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

Erroneous constraining of estimable quantities distorts the geometric

and/or physical characteristics of the dynamic environment, thus if

enough time is allowed the accumulated errors can exceed the noise level

of the observations to the extent that the adjusted estimable quantities

will be affected by errors larger than those implied by the accuracy of

the observations. Even worse with a priori variance of unit weight

equal to one, the a posteriori standard deviations of these quantities

will also be reduced only if the a posteriori variance of unit weight is

close enough to unity so it preserves the relative structure of the

weight coefficient matrices (ibid.)

Under those conditions baselines estimated on the basis of Cartesian

coordinates will have their a posteriori variances reduced only if the

non-linearity of the models employed allow such reductions to take

place. In single SRD solutions the errors affecting the baselines,

computed from the erroneously adjusted Cartesian coordinates of the

free baseline end, result from the errors originated by erroneously

constraining the parameters entering in the computation of the

discrepancy vector L of equation (2-60). Erroneous baseline lengths

with reduced variances constitute a potential source for misleading

inferences in regard to the accuracy of the estimated baselines.

Therefore, the a posteriori standard deviations of the estimated

baselines will only be used, as in the geometric solutions, to indicate

how far away the solution is from its steady state. The actual accuracy

of the estimated baselines will be inferred on the basis of the

comparison with the baselines estimated through geometric method

(present study) and range dynamic methods (CSR at UTX and ZIPE at

Potsdam) (next section).

Baselines estimated through Cartesian coordinates will also be

affected in the dynamic environment by the errors commited in the

implementation of the terrestrial reference frame, as it happens, for

instance, when the barycenters of the observations of the baseline end

stations are well apart in time (Pavlis et al., 1983). Strict simultaneity,

however, implied by the nature of the SRD observable eliminates such

inconsistencies in the present study.

143
For various baseline lengths ranging from 8 m to 4000 km, Tables 9

through 17 list the final results and those obtained at intermediate

stages of the process leading to a steady state response. The a priori

variance of unit weight is taken equal to one for all of the baseline

solutions presented in this section. Fig. 18 shows the locations of the

stations involved in the SRD solutions (lower part) together with a

typical LAGEOS groundtrack for some of their coobserved passes (upper

part).

The lengths of these common parts increase or decrease as the

length of the baseline being estimated decreases or increases

respectively.

Tables 9 and 10 contain the results of the steady state response

reached for two very short baselines via both short and long arc

modes. The first column lists the number of passes coobserved by the

baseline end points while the second lists the number and the duration

of the arcs, the position and orientation of which were adjusted to fit

"best" the available observations. For instance, 6 (lh) in the first row

of Table 9 indicates that six arcs were adjusted, each of one hour long,

and 1 (4h) and 1 (2d) in the fifth row indicate that observations from

eight passes (first column of this row) were adjusted using two arcs,

one four hours long, and one two days long. The thirdp fourth and

fifth column list the total number of the SRD observables, the a

posteriori variance of unit weight and the root mean square of all of the

SRD residuals respectively. The last column contains the estimated

baseline length together with its a posteriori standard deviation. All of

the remaining tables through 17 have the same format except for some

obvious very minor differences.

For baselines 7109 - 7886 and 7110 - 7220 steady state response has

been reached via both short arc and long arc solutions on the basis of

10 passes (Table 9, rows 3-7) and 17 passes (Table 10, rows 2-5)

respectively. As in the geometric solutions, the steady state response is

also associated with a posteriori standard deviations below the

centimeter level. After the steady state, the estimated lengths of these

two baselines do not change at the centimeter level either through a

144

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I _fi_IGINA[7 PPfGI_ TS

80 ° ' ' ' ' ' ' ' _, 80"

I T"- ......

' 1 7 30

I . l le_V,',/;_ I_.: I I,o.

| ! i 1 - " "

, T

?_,oo- u_ qo° I _o. x.c.o _7o_e -- ( so'

i ,_:.

20" '_,k [ 20"

I0" 10"

I o. ., __ o.

( . ",_

:,_ / .

i -_o._ _ _ o _ o _ _ o _ k_ _, _ -o

Fig. 18 Station location and Lageos groundtracks

145

LI
Table 9

No. of

Steady State Response of Baseline 7109 - 7886, Parallel

Passes

All length units in meters.

Integration No. of !2) P_S of All SRD Baseline

Passes Lengths (:) Observ. _ Residuals 7.

I

I

I

Short-Arc Mode I

6 6 (1 h) 8,414 0.94 0.070 .719 ± 0.002

8 8 (1 h) 12,807 0.90 0.069 .731 ± 0.002

10 10 (1 h) 18,589 0.86 0.067 .737 * 0.001

13 13 (1 h) 25,865 0.85 0.067 .738 ± 0.001

£ong-ArcMode

8 1 (4 h) 13,527 0.90 0.069 .733 ± 0.002

1 (2 d)

l0 1 (4 h) 20,216 0.86 0.067 .741 ± 0.001

i (4 d)

13 1 (4 h) 27,697 0.85 0.067 .738 * 0.001

1 (7 d)

(,) k (_ h) = k arcs of _ hours

I

I

!

I

i

k (_ d) = k arcs of $ days

(2) a posteriori variance of unit weight

short arc or a long arc solution and furthermore these lengths are the

same for both short arc and long arc solutions, thereby indicating that

the accumulated orbital biases for so short baselines cancel out

completely, a plausible property of the SRD observable. The correlations

among the components of the adjusted parameter vector are for so short

baselines substantially reduced because the lengths of the coobserved

parts of the passes tend to be long and their orientation is not

important since the subsatellite tracks for these short baselines are

located well away from their end points (Section 4.4.l).

The 34 cm a posteriori standard deviation of the baseline solution,

listed in the first row of Table 11, shows that although 17 passes of

about one hour long are available, steady state response, with a short

arc solution, is not possible for baseline 7110 - 7265 because the

pass-baseline geometry has deteriorated as the length of the estimated

146

I

!

!

I

!

I

I

I

I
I

!

I

I

!

I

I

!

!

I

!

!

I

I

I

!

I

I

Table 10 Steady State Response of Baseline 7110 - 7220, Parallel

Passes

All length units in meters.

No. of Integration No. of ¢_(2) P_4S of All SRD Baseline

Passes Lengths (*) Observ. Residuals 15.

Short-ArcMode

ii ii (i h) 12,040 0.95 0.084 .233 ± 0.003

17 17 (1 h) 17,971 1.09 0.091 .236 ± 0.002

Long-ArcMode

17 1 (4 h) 17,982 1.07 0.090 .240 ± 0.002

1 (1 d)

1 (2 d)

3 (3 d)

17 1 (I h) 17,983 1.08 0.090 .250 ± 0.002

1 (3 d)

1 (4 d)

2 (5 d)

17 1 (1 h) 17,984 1.08 0.090 .238 ± 0.002

1 (3 d)

1 (S d)

1 (7 d)

(*) k ($ h) = k arcs of _ hours

k (_ d) = k arcs of _ days

(2) a posteriori variance of unit weight

baseline has increased from 8 m (Tables 9) to about 274 km (Table 11).

This deterioration is even worse since the geometry of the orbit is such

that passes parallel to this baseline do not exist and only passes

intersecting it at about ± 30 to ± 50 degrees are available creating,

therefore, a geometry which is worse than that of the parallel passes

and better than that of the perpendicular ones. Steady state response,

however, for this solution can be reached by strengthening its

geometric characteristics on the basis of constraints imposed on some

additional estimable quantities. In the dynamic environment this is

accomplished by increasing the maximum length of continuous integration

(maximum arc length) which intensifies, through the implied geometric

strength, the effect of the coefficients Anm and Bnm on the resulting

long orbital arcs. Care should be excercised not to destroy the

147
Table ii Steady State Response of Baseline 7110 - 7265, Passes

Within ±30" - *50"

All length units in meters.

Passes

No. of Integration No. of _=) P_MS of All SRD Baseline

Lengths (*) Observ. Residuals 2740

I

I

I

Short-Arc Mode I

17 17 (1 h) 21,767 0.77 0.0734

Long-Arc Node

17 1 (1 d) 21,780 0.80 0.075

3 (2 d)

17 2 (2 d) 21,781 0.80 0.075

1 (4 d)

17 1 (2 d) 21,781 0.80 0.075

1 (3 d)

1 (7d)

70.453 ± 0.335

69.391 ± 0.012

69.482 ± 0.009

69.494 ± 0.009

(,) k (_ h) = k arcs of $ hours

I

I

!

I

k (_ d) = k arcs of _ days

(=) a posteriori variance of unit weight

accuracy of the available observations by the errors accumulated, over

these long arcs, due to errors affecting the Anm and Bnm potential

coefficients.

With this in mind and on the basis of all observations, three long

arc solutions have been performed by allowing maximum arc lengths of

up to three, five and seven days respectively. The results of these

three solutions are listed in 2nd through 4th row of Table 11. In the

first long arc solution, employing one arc of one day long and three

arcs of two days long, the length of the estimated baseline changed by

106 cm as compared to that of the short arc solution. Steady state

response, however, for this solution has not yet been reached since the

length of this baseline changes by about 10 cm when arcs up to five

and seven days are allowed {Table 11, rows 3-4). However, increasing

the maximum arc length from four to seven days (rows 3-4) results in a

change of the estimated baseline length of about 1 cm, thereby implying

that steady state response is being reached and that the baseline length

of the solution listed in 3rd row of this table is the least affected by

148

!

I

!

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

l

I

I

I

Table 12 Steady State Response of Baseline 7109 - 7110, Parallel

Passes

All length units in meters.

No. of Integration No. of (2) _MS of All SRD Baseline

Passes Lengths (*) Observ. _ Residuals 883602.

Short-Arc Mode

4 4 (1 h) 11,512 0.5 0.029 .676 ± 0.30

8 8 (1 h) 36,236 0.36 0.026 .309 ± 0.04

12 12 (1 h) 59,011 0.35 0.026 .249 ± 0.03

16 16 (I h) 69,083 0.43 0.026 .251 ± 0.03

Long-Arc_e

12 1 (1 h) 59,020 0.36 0.026 .220 ± 0.001

i (2 d)

1 (3 d)

12 1 (4 h) 59,021 0.36 0.026 .224 ± O.OO1

1 (3 d)

12 i (7 d) 59,022 0.86 0.042 .217 ± 0.002

16 2 (2 d) 69,81.6 0.37 0.027 .225 • 0.001

i (3 d)

16 1 (2 d) 69,817 0.38 0.027 .226 ± O.OO1

1 (4 d)

16 1 (1 d) 69,817 0.79 0.057 .219 ± 0.001

1 (6 d)

(*) k (_ h) = k arcs of t hours

k (_ d) = k arcs of _ days

(2) a posteriori variance of unit weight

any accumulated orbital errors, because at steady state this solution

employs the shorter long arcs and the a posteriori variance of unit

weight is smaller than one. The correlations among the components of

the adjusted parameter vector do not cause, for this baseline, any

instability since almost all of them are less than 0.90 with very few just

exceeding this value.

Tables 12 and 13 contain the results for the steady state response

of two baselines with the same length, the same orientation, and one

common end point occupied by station 7110. The observations recorded

149
Table 13 Steady State Response of Baseline 7110 - 7886, Parallel

and Perpendicular Passes

All length units in meters.

No. of Integration No. of (2) }{MS of All SRD Baseline

Passes Lengths (*) Observ. ¢_ Residuals 883606.

I

I

I

Short-Arc Mode I

33 33 (1 h) 58,261 0.74 0.064 .467 * 0.056

parallel

passes

33 39 (I h) 61,037 0.75 0.064 .459 ± 0.056

parallel

+6

perpendicular

£ong-Jh'cMode

33 4 (1 h) 58,248 0.75 0.065 .347 ± 0.004

parallel 4 (4 h)

passes 2 (1 d)

1 (1.5 d)

3 (2d)

33 2 (1 h) 58,252 0.76 0.065 .342 ± 0.003

parallel 1 (4 h)

passes i (i d)

3 (2 d)

3 (3 d)

33 1 (4 h) 58,256 0.77 0.065 .335 ± 0.003

parallel 1 (2 d)

passes 1 (3 d)

1 (4 d)

i (4.5 d)

1 (7 d)

I

I

I

I

I

i

I

I

(*) k (_ h) = k arcs of _ hours

k ($ d) = k arcs of $ days

(2) a posteriori variance of unit weight

by stations 7109 and 7886, occupying the other ends of these two

baselines, have an accuracy of 0.028 m and 0.070 m respectively,

therefore, examination of the results presented in these two tables will

show how the quality of the observations affects the speed of the

steady state response (Section 4.4). Steady state response of baseline

7109 - 7110 has been reached via short arc solution on the basis of 12

150

I

!

I

I

!

I
I

I

I

I

l

I

I

l

I

I

I

I

I

I

I

i

I

I

passes having approximately 59,000 observations, as is confirmed by

comparing the third and fourth rows of Table 12. The a posteriori

standard deviation of the recovered baseline is at the 3 cm level

because, in short arc solutions, the components of the initial state

vectors are recovered with a relatively low precision (20 m to 50 m).

To examine the effects of using long arc solutions to achieve steady

state response, we have performed six long arc solutions, three of

which are based on observations from 12 passes (Table 12, rows 5-7)

and three on observations from 16 passes (rows 8-10} when maximum arc

lengths up to three, five and seven days are allowed respectively. The

baseline lengths, estimated through all of these long arc solutions, do

not change their length at the centimeter level, when either the number

of observations or the maximum arc length increases and furthermore,

the a posteriori standard deviation of the estimated baselines assumes

millimeter level values. However, the baseline length changed from the

short arc solution to the long arc solution by about 2.5 cm due to

accumulation of orbital errors. The RMS of all of the SRD residuals for

the long arc solutions (Table 12, rows 5-6,and 8-9) is about the same as

that of the short arc ones (rows 1-4),but it is twice as much when arcs

up to seven and six days are allowed in those solutions {rows 7 and 10).

Steady state reached through a short arc solution assumes a

posteriori standard deviation of 3 cm (Table 12), thus, for baseline

7110-7886 shown in Table 13 such a response has not yet been reached,

because a short arc solution, on the basis of all of the available

observations, results in a 6 cm a posteriori standard deviation of the

estimated baseline length (Table 13, row 1). This claim is also confirmed

by the 10 cm change of the estimated baseline length when long arc

solutions are performed (rows 3-5). The perpendicular passes

incorporated in the short arc solution (2nd row}, tend, in the light of

weak geometry, to decrease the estimated baseline length because for

those passes the SRD observable tends to become very short (Section

4.4.1). Thus, the 8 mm decrease (2nd row) on the basis of just 2,776

additional observations is another indicator confirming that steady state

151
response for this baseline through a short arc solution has not been

reached. Consequently, such a response is sought through long arc

solutions performed on the basis of arcs up to two, three, and seven

days long {Table 13, rows 3-5}. Since in these three long arc solutions

the estimated baseline length did not change at the cm level and since

its a posteriori standard deviation is below the cm level, it is assumed

that steady state has been reached in those solutions.

For this baseline, although of same length and direction as baseline

7109 - 7110_ steady state response was not possible on the basis of 33

passes with ._8,261 observations processed via a short arc solution,

because the accuracy of the observations of station 7886 is twice as

large as those of station 7109. This confirms the claim, made in Section

4.4, that steady state response on the basis of improved quality

observations is faster.

As of the correlations among the components of the adjusted

parameter vector, they are well below the 0.90 value for all short arc

solutions shown in Tables 12 and 13, except for very few of them,

among the components of some of the recovered state vectors, that

tend to be just a little higher than this value without, nevertheless,

affecting the conditioning of the normals to the extent that could result

in an algorithmically singular normal equations matrix. This behavior of

the correlations is not surprising since in estimating these two baselines

only parallel passes, the coobserved parts of which are extended well

beyond the baseline end points, were employed, thereby, leading

according to the discussion in Section 4.4.1 to a favorable geometry

which is manifested in the reduced correlations among the components of

the adjusted parameter vector. On the contrary, when in some of the

long arc solutions shown in Tables 12 and 13, arcs of one or four hours

long are employed together with arcs of two, three or seven days long,

the relative geometry of the hours long arcs is very weak as compared

to that implied by the days long arcs, thereby leading to high

correlations among the components of initial state vectors of the hours

long arcs. This inhomogeneity together with the potential to have many

passes with weak geometry, when long baselines of about 4000 km are

152

I

I

I

I

I

I

|

I

I

!

!

!

I

I

I

I

I

I
I

I

I

I

l

I

I

l

I

I

I

I

I

I

I

I

I

I

estimated, results in an algorithmically singular normal equations matrix

(see below).

Baseline 7110 - 7086, the steady state response of which is shown in

Table 14_ has a relatively large number of coobserved passes each

having a substantially reduced number of observations as opposed to

those recorded in the passes used in the estimation of the baselines

described in Tables 9 through 13. Examination of the steady state

response for this baseline will reveal the effect of the geometry, implied

by the many passes, as opposed to the number of the available

observations. The 41 cm a posteriori standard deviation obtained with

the short arc solution, shown in the first row of this table, indicates

that steady state has not yet been reached. Thus, six long arc

solutions have been performed, the three of which use observations from

29 passes (Table 14, rows 2-4) while the other three use observations

from 40 passes (rows 5-7) and arc lengths of up to three, five and

seven days respectively. On the basis of 29 passes, processed in the

long arc mode, steady state response has not been reached because

changing the maximum arc length from five days to seven days the

length of the estimated baseline changes as much as 18 cm (Table 14,

rows 3-4). However, when 11 more passes with about 4,000 additional

observations are included, the lengths of the estimated baseline change

only at the centimeter level (rows 4-7) and their a posteriori standard

deviations has dropped below the centimeter level, thereby suggesting

that steady state response for this baseline has been reached. Since

the a posteriori variance of unit weight is smaller than one, the

estimated lengths are assumed not to have been influenced by any

errors accumulated over these long periods and therefore, the length

estimated with the smaller standard deviation is assumed to be the

closest to that of the steady state (last row of this table).

In the short arc solution (lst row of Table 14), the correlations

among the components of the adjusted parameter vector are almost all

of them less than 0.90 with few just exceeding this value. These

correlations (i.e., >0.90) exist among some components of the initial state

vectors corresponding to arcs, the geometry of which is not favorable

153
Table 14 Steady State Response of Baseline 7110 - 7086, Passes

Within ±20" - ±60"

All length units in meters.

No. of Integration No. of ¢_(2) RMS of All SHD Baseline

Passes Lengths (*) Observ. Residuals 119829

I

I

I

Short-Arc Mode I

40 40 (i h) 16,752 0.76 0.078

Zong-Arc

29 8(lh), l(4h) 12,395 0.81 0.081

3(ld), 2(2d)

3(3d)

29 5(lh), 2(4h) 12,414 0.82 0.081

2(ld), 2(2d)

2(4d), 2(5d)

29 4(lh), l(ld) 12,400 0.84 0.083

l(2d), 2(4d)

2(5d), 2(7d)

40 lO(lh), l(4h) 16,771 0.82 0.081

3(ld), 2(2d)

5(3d)

40 7(lh), l(4h) 16,790 0.80 0.080

3(Id), 3(2d)

l(3d), 2(4d)

2(5d)

40 7(lh), 1(ld) 16,776 0.87 0.084

1(2d), 1(4d)

2(5d), 2(6d)

2(7d)

1.110 ± 0.411

0.814 ± 0.042

0.821 ± 0.036

0.998 ± 0.015

0.987 ± 0.009

1.024 ± 0.009

1.005 ± 0.008

I

I

I

I

i

!

I

I

(*) k (_ h) = k arcs of _ hours

k (_ d) = k arcs of _ days

(a) a posteriori variance of unit weight

with respect to the estimated baseline (Section 4.4.1), that is, when the

passes are close to being perpendicular to the baseline and/or when

their coobserved part is short and close the perpendicular bisector

plane of the estimated baseline.

Increasing the maximum arc length results in arcs of different

lengths ranging from one hour to either three, five, or seven days long,

all of which are employed in the same long arc solution. This in turn

154

I

I

I

l

I

I
I

I

I

I

!

I

I

l

I

I

I

I

I

I

I

I

l

l

weakens the geometric strength of the short arcs because in those

solutions the long arcs are the ones dominating the implied geometric

strength. The weaker the geometry of those short arcs, in a short arc

solution, the higher the correlations among their components adjusted in

long arc solutions, in which these short arcs could not be matched with

any other arcs. Furthermore, if only passes of weak geometry are

matched together they result in a long arc also of weak geometry,

thereby leading again to high correlations among the components of its

adjusted initial state vector. These correlations become larger as the

maximum arc length increases and the long arcs of weak geometry

cannot be matched with any other ones of strong geometry. This

pattern, which is present in the long arc solutions of Table 14, does

not lead to an algorithmicaUy singular normal equation matrix, simply

because the length of the baseline is not long enough to result in a

geometry so weak that could lead to near singularities, although passes

parallel to this baseline do not exist (Fig. 18).

Table 15 contains the results for the steady state response of

baseline 7110 - 7122_ the geometry of which is such that passes

parallel to this baseline do exist and in fact all of the passes employed

in the solutions of Table 15 are parallel to this baseline. Since steady

state response for this baseline was not possible through a short arc

solution, long arc solutions have been performed on the basis of 15

passes {rows 1-2) and 22 passes (rows 3-5) with maximum arc lengths

up to three, five and and seven days respectively. The temporal

distribution of the observed passes is such that when 15 of them are

used the solution allowing maximum length of three days coincides with

that allowing maximum arc length of five days, and therefore only one of

them is shown in Table 15 {lst row}. These passes contain observations

from station 7122 before and after it was upgraded (Table 2).

The 3 mm and 2 mm level of the a posteriori standard deviations

associated with the baselines obtained through all of the long arc

solutions, shown in Table 15, indicate that a steady state response for

this baseline has been reached. This is also confirmed by the fact that

the estimated baseline length changes at just the centimeter level when

155
Table 15 Steady State Responseof Baseline 7110 - 7122, Long-Arc

Mode, Parallel Passes (all length units in meters)

No. of Integration No. of

Passes Lengths (*) Observ. a_(2) P_MS of All SRD Baseline

Residuals 1437139. I

15 1 (1 h) 42,328 0.65 0.044 .307 ± 0.003

5 (1 d)

1 (2 d)

15 1 (1 h) 42,329 0.66 0.045 .309 ± 0.002

4 (1 d)

i (6 d)

22 2 (1 h) 69,803 0.52 0.035 .305 ± 0.002

4 (i d)

3 (2 d)

1 (3 d)

22 1 (1 h) 69,804 0.54 0.037 .302 ± 0.002

5 (i d)

2 (2 d)

i (3 d)

22 2 (l d) 69,807 0.81 0.055 .293 ± 0.002

1 (2 d)

I (6 d)

2 (7 d)

I

I

I

I

I

I

(*) k ($ h) = k arcs of _ hours

k (_ d) = k arcs of $ days

(2) a posteriori variance of unit weight

on the basis of 22 passes the maximum arc lengths change from 3 to 7

days (Table 15, rows 4-5). The correlations follow the same pattern as

that described for baseline 7110 - 7086 but since the passes are now

parallel to the estimated baseline, weak geometry is implied only by

those the coobserved part of which is extended only in between the

baseline endpoints (Section 4.4.1). However, the geometry implied by

those passes is not so weak to cause any algorithmic singularity in the

normal equations matrix. Such weakness in the geometry causes

problems when the lengths of the estimated baselines increases to

3500 km and 3700 km (Tables 16 and 17).

156

I

I

I

I

I

I

I

I

I
I

I

I

Table 16 Steady State Response of Baseline 7109

Mode, Passes Within *20" - *30"

All length units in meters.

- 7105, Long-Arc

(2) NMS of All SRD Baseline No. of Integration No. of ¢_ Passes Lengths (1) Observ. Residuals 3703351.

I

I

I

l

l

I

l

I

I

I

I

I

I

I

l

Naxim.m Arc Iength= 2 days

42 10(lh), 5(4h) 78,801

2(7h), 6(ld)

54 12(lh), 5(4h) 106,770

2(7h), 8(ld)

62 14(lh), 5(4h) 129,147

2(7h), 9(ld)

72 18(lh), 5(4h) 151,901

2(7h),ll(ld)

Maximm Arc I,en&,'_h= 3 days

42 5(lh), 2(4h) 72,783

l(7h), 4(ld)

2(2d), 3(3d)

54 5(lh), 2(4h)

l(7h), 4(ld)

4(2d), 4(3d)

62 5(lh), 2(4h)

l(7h), 5(ld)

5(2d), 4(3d)

72 6(lh), 2(4h)

l(7h), 7(ld)

6(2d), 5(3d)

0.54 0.032 .628

0.52 0.032 .686

0.50 0.031 .702

0.49 0.031 .693

0.57 0.033 .677

106,989 0.58 0.033 .689

129,155 0.54 0.032 .699

151,910 0.53 0.032 .695

157

* 0.003

* 0.003

* 0.002

* 0.002

* 0.003

* 0.002

* 0.002

* 0.002
Table 16 (cont'd)

No. of Integration No. of (2) RMS of All SRD Baseline

Passes Lengths (*) Observ. _ Residuals 3703351. I

MaxiJ_mArcZe_gth = 7 days

42 4(lh), 1(4h) 78,809 1.11 0.046 .706 ± 0.004

I(7h), l(ld)

l(2d), l(3d)

l(5d), 4(6d)

l(7d)

54 4(15), 1(4h) 106,780 3.33 0.081 .815 ± 0.005

1(75), l(ld)

l(2d), l(3d)

l(5d), 5(6d)

2(7d)

62 4(lh), l(4h) 129,159 2.78 0.073 .809 ± 0.005

1(7h), l(ld)

l(2d), l(3d)

l(5d), 6(6d)

2(7d)

72 4(15), 1(45) 151,916 2.45 0.069 .801 ± 0.004

1(7h), 2(ld)

l(2d), l(3d)

l(4d), l(5d)

7(6d), 2(7d)

I

I

I

I

I

l

I

(*) k (_ h) = k arcs of _ hours

k (_ d) = k arcs of t days

(2) a posteriori variance of unit weight

For both of these baselines the pass-baseline geometry is weak and

that of baseline 7110 - 7105 is even worse since passes intersecting this

baseline at ±20 to ±60 degrees are only possible, and many of them are

located close to the perpendicular bisector plane of this baseline.

For both of these baselines steady state response was not possible

through a short arc solution, therefore, long arc solutions were

performed on the basis of 42, 54, 62 and 72 passes and with maximum

arc lengths up to two, three, and seven days respectively.

Examination of the results listed in rows 2-4 of Table 16 reveals

that with maximum arc length up to two days steady state response has

158

I

I

I

I

I

I

I

I
I

I

I

I

i

I

I

l

I

I

I

I

I

I

I

I

I

I

been reached when 54 passes are available and the a posteriori standard

deviation of the estimated baseline is 3 mm. Increasing the maximum arc

length up to three days leads, as expected, to a faster steady state

response since such a response is close to being reached with 42 rather

than with 54 passes (Table 16, rows 2 and 5). When maximum arc length

of 7 days is allowed, we see that steady state response has clearly been

reached on the basis of 42 passes (Table 16 cont'd, row 1). However,

the a posteriori standard deviation is greater than one, thereby

suggesting that either accumulated orbital errors and/or possible

ill-conditioning may be on their way up to corrupt the solution.

Ill-conditioning, however, seems to be at work since many components of

the recovered initialstate vectors exhibit correlation at the .9999 level

resulting mainly from the weak geometry and the inhomogeneity of the

arc lengths employed in the solution shown in the 1st row of Table 16

cont'd. This is also confirmed, when on the basis of 54 passes (2nd

row), the estimated baseline length jumps by II cm and the a

posteriori standard deviation of unit weight jumps to 3.33. These

behavior is caused by the high correlations existing among the

components of the initial state vectors of those short or long arcs, the

geometry of which became worse by the additional long arcs when 54

rather than 42 passes were employed (Table 16 cont*d, rows 1-2).

Inclusion of additional observations, however, seems to lead the response

of this solution to the right direction, as it is seen from the 2nd

through 4th row of this table, but at a very slow pace, since the

geometry is stilltoo weak for a steady state response to take place.

The divergent response resulting from weak geometry manifests itself

when baseline 7110 - 7105 is estimated on the basis of 56, 60 and 65

passes and with maximum arc lengths up to three, five and seven days

respectively. For many of these passes the geometry with respect to

baseline 7110 - 7105 is very weak (Section 4.4.1). This in turn leads to

high correlations among the components of the corresponding initial

state vectors. With arc lengths up to three days, steady state response

seems to have been reached because with the incorporation of additional

observations the length of the estimated baseline just changes at the

159
Table 17 Steady State Response of Baseline 7110 - 7105, Long-Arc

Mode, Passes Within *20" - *60"

All length units in meters.

No. of Integration No. of _(2) RMS of All SRD Baseline

I

I

I

Passes Lengths (i) Observ. Residuals 3559743.

Maxi_u_ Arc I,ength = 3 days I

56 12(lh), l(5h) 120,111 0.58 0.036 .594 ± 0.003

l(8h), 5(ld)

3(2d), 4(3d)

60 12(lh), l(5h) 136,356 0.56 0.036 .615 ± 0.003

l(8h), 4(ld)

3(2d), 5(3d)

65 15(lh), l(5h) 146,956 0.55 0.035 .619 ± 0.003

l(8h), 5(ld)

3(2d), 5(3d)

_i_ArcI, e_g_h = 5 days

56 ll(lh), l(5h) 120,114 0.64 0.038 .542 ± 0.003

l(8h), 3(ld)

2(2d), l(3d)

3(4d), l(5d)

60 ll(lh), l(5h) 136,359 0.63 0.038 .554 z 0.003

l(8h), 2(ld)

2(2d), 2(3d)

3(4d), l(5d)

_ei_ArcZe.ng_ = 7 days

56 8(lh), l(ld) 120,118 1.25 0.053 .637 ± 0.004

l(2d), 2(4d)

5(6d), 2(7d)

60 8(lh), l(ld) 136,362 1.14 0.051 .615 ± 0.004

2(2d), 2(4d)

5(6d), 2(7d)

I

I

I

I

I

I

I

I

I

(1) k (_ h) = k arcs of _ hours

k (_ d) = k arcs of $ days

(2) a posteriori variance of unit weight

centimeter level and the associated standard deviations assume 3 mm

values (Table 17, rows 1-3}. However, correlations at the 0.9999 level

exist among the components of the initial state vectors of those arcs

that are primarily composed of passes crossing this baseline close to its

160

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

perpendicular bisector plane. These correlations become larger when

maximum arc lengths up five and up to seven days are allowed and

when additional observations result in sharp inhomogeneity in regard to

the arc lengths incorporated in the same solution (Table 17, rows 4-5,

and 6-7). Thus, on the basis of 65 passes and when arcs up to five or

seven days are allowed, the existence of high correlations leads to a

singular normal equation matrix, thereby implying divergent response.

Because of this, solutions were not possible with those arc lengths and

therefore they are not shown in Table 17.

Thus, it is questionable whether for long baselines (>3500 kin)

steady state response through single SRD baseline solutions can be

reached. For shorter baselines, however, steady state response can be

reached even if passes parallel to the baseline do not exist.

The ill-conditioning affecting longer baselines could be prevented in

a network solution if there exist observations from baselines that are

parallel to the passes responsible for the ill-conditioning. Such

solutions, however, for the reasons mentioned in the previous section

could not be performed in the course of this study.

4.5 Baseline Comparison

The baseline lengths estimated with the geometric method are

independent of any orbital errors and inconsistencies affecting the

implementation of the terresrial reference frames; therefore, these

lengths (Section 4.3.2) will constitute the standards of comparison when

assessment is made about the accuracy of the baselines estimated via

the SRD method. Both the geometric and the SRD baseline estimates will

also be compared with those obtained via long arc range dynamic

methods, on the basis of the MERIT data set, by the Central Institute

for Physics of the Earth (ZIPE) in East Germany and by the Center for

Space Research (CSR) at The University of Texas.

Both of these centers used in their solutions the gravitational

constant proposed by MERIT standards (i.e., GM =3.98600448 * 10 _"

m_s -2) which is different from that employed in the single SRD baseline

solutions (Section 2.2.5.1). A change in the gravitational constant will

161
not only affect the metric scale (the semimajor axis a) but also the

dynamic scale (mean anomalistic motion n) according to Kepler's modified

third law :

n2a s =GM(I+_) (4-9)

where _ is a small parameter depending on the satellite orbit. It is

uncertain how the effect of this change will be divided between the

metric and the dynamic scales. For regional baseline estimation

however, the metric scale is primarily controlled by the velocity of light

which is the same for all the SRD, ZIPI_ and (CSR)85L01 solutions. For

this reason, changing the value of the gravitational constant (GM) in the

SRD solutions presented in the previous section in accordance with

MERIT Standards will affect at the centimeter level the lengths of only

the two very long baselines {i.e., 7109-7105 and 7110-7105). The

remaining baselines are affected either at the few millimeter or

submillimeter level depending on the length of the estimated baseline

and on the duration of the arcs employed in the corresponding solution.

Nevertheless, for the sake of consistency the baseline differences shown

in Table 18 have been adjusted to the same value of the gravitational

constant by reestimating the SRD baselines on the basis of the (GM)

value proposed by the MERIT Standards.

Table 18 lists the differences of those baselines only for which

steady state response has been reached either through the SRD or the

geometric solutions. The first column of this table lists the baseline

ID's, the second one contains the length of the baseline rounded to the

nearest meter, the third and fourth columns contain the differences of

the baselines estimated through the SRD method and through the range

dynamic method by both ZIPE and CSR, the fifth and sixth columns list

the differences of the geometric and dynamic ZIPE and CSR baselines.

The last column contains the differences of the baselines obtained

through the SRD and the geometric methods. The baseline differences

not listed in this table were not computed because steady state response

for the SRD and/or the geometric method was not possible or because

one of the stations constituting a baseline endpoint was not included in

162

i

I

I

I

I

I

l

I

I

I

I

I

I

I

I

I

I

I
1

! Table 18 Baseline Differences

Length (m) ZIPE (CSR)85L01 ZIPE (CSR)85L01 SRD

Baseline to nearest -SRD -SRD -GEOM -GEOM -GEOM

meter

II

II

II

I

II

II

7109-7110 883602. -0.03 -0.03 -0.03 -0.03 0.00

7109-7265 627044. -- m 0.00 0.04

7109-7886 8. -O.O1 0.00 -0.02 -O.O1 -O.O1

7110-7122 1437139. -0.03 -0.03 -0.02 -0.02 O.O1

7110-7220 15. m -0.04 m -O.O1 0.03

7110-7265 274069. -0.04 -0.07 -0.03 -0.06 O.O1

7122-7265 1663983. -- _ -0.02 -0.06

7122-7886 2280720. _ _ 0.01 0.02

7220-7265 274066. m m _ 0.00

7265-7886 627049. _ _ 0.01 0.06

7110-7886 883606. -O.O1 0.00 -O.O1 0.00 0.00

7110-7086 1198291. 0.01 0.02 _ m

7109-7105 3703352. 0.06 0.05 _ _

I

II

II

I

I

I

I

I

II

the corresponding dynamic method. For instance, the difference (ZIPE -

SRD) for baseline 7220-7265 could not be computed because steady state

response for this baseline was not possible through a single SRD

baseline solution; also for the same baseline the difference (ZIPE -

GEOM) is not included in Table 18 because station 7220 was not listed in

the reported ZIPE solutions (Montag et al., 1985).

The baseline differences between the SRD, ZIPE and CSR solutions

shown in the third and fourth columns are negative for north-south

baselines (rows i-II) and positive for the east-west ones (rows 12-13)

(Fig. 18). Since some of these differences are larger in magnitude for

shorter baselines, scale difference between SRD and CSR or ZIPE

solutions would account for part of these differences. The remaining

differences, at the 2 cm or 3 cm level, must be caused by orbital errors

affecting mostly the dynamic solutions, because the SRD solutions

(column 7) are clearly closer to the geometric solution than both ZIPE

and CSR solutions (Table 18, columns 5-6). The large differences of 6

163
cm between (CSR)85L01and geometric solutions (column 6) are associated

with station 7265 (MOHAVE). This station during the MERIT Main

Campaign was equipped with a TLRS-1 laser instrument which

experienced many problems, as was confirmed in the present study when

editing the data with the data snooping procedure (Section 3.5).

Therefore, it is very likely that erroneous observations from this station

are still present in the (CSR)85L01 solution. This in turn would explain

these relatively large baseline differences.

Although the geometric solution (Table 8, Type C3) used in Table 18

has been overconstrained to the (CSR)85L01 solution, SRD baseline

estimates are on the average the "best" ones as compared to those of

either the ZIPE or the (CSR)85L01 solutions. The term "best" indicates

that a solution is the one closest to a geometric solution at its steady

state. However, ZIPE baseline estimates when compared to those

obtained in the (CSR)85L01 solution, are closer to the geometric solutions

and therefore more accurate.

The root mean square of the differences between the geometric

baseline estimates and those of the SRD, ZIPE and CSR are l cm, 2 cm,

and 4 cm respectively.

Since these differences are based on baselines up to 1500 kin, it is

fair to state that SRD baseline estimates are at least as accurate as

those obtained through the range dynamic mode methods.

4.6 RESPONSE OF THE SRD METHOD TO THE SIMPLIFICATIONS OF THE

ORBITAL MODEL

It was mentioned in Section 2.2.5 that the aim of the present study is

not to estimate Lageos' orbit with the highest degree of accuracy but

rather to employ models as simple as possible and yet be able to

recover baselines with an accuracy compatible to that of the

observations.

Since the temporal variations of the baseline endpoints have been

accounted for to the required degree of accuracy (Section 3.3), and

since any inconsistencies in the implementation of the Terrestrial

Reference Frame do not affect the SRD observables (Section 4.4.2), the

164

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

l

I

I

l

I

I

I

I

I

I

I

I

I

I

I

I

I

errors affecting the baselines estimated via the steady state response of

the SRD method {Section 4.4.2) result only from those orbital errors

accumulated over the integration periods and not cancelled out in the

computation of the SRD observable. Thus, the question that should be

addressed and investigated is twofold:

1) Is the sophistication of the orbital model, employed in the present

study (Section 2.2.5) sufficient to result in baselines the accuracy of

which is compatible to that of the observations?

2) If the answer is yes, how much could the sophistication of the orbital

model be reduced without affecting the accuracy of the estimated

baselines? If the answer is no, how much should the sophistication

of the orbital model be enhanced so the estimated baselines have an

accuracy compatible to that of the observations?

The errors affecting the estimated baselines are propagated from

those affecting the corresponding Cartesian coordinates. These errors,

which originated by the erroneous constraints imposed on a large

number of estimable quantities entering in the computation of the SRD

observable (eq. 2-60), tend to accumulate as the employed integration

periods become longer.

For baselines of moderate length {<2000 kin) accumulated radial

errors are cancelled out almost totally in the computation of the SRD

observable. These errors are propagated almost unaltered into the

computed range observable (Pavlis, 1982). Depending, however, on the

location of the observed satellite positions, accumulated latitudinal and

longitudinal errors may affect the computed value of the SRD observable

worse than they affect the range observable {ibid.} For shorter

baselines (<200 kin) the computed SRD observable is less affected by all

these three errors {i.e., radial, latitudinal and longitudinal).

Consequently, the answers to the questions posed in the beginning of

this section depend on the length of the estimated baseline, on the

magnitude of the accumulated orbital errors, and on whether these

errors are radial, latitudinal or longitudinal. If these errors are

primarily latitudinal and longitudinal, these answers depend also on the

165
relative orientation of the estimated baselines and its observed passes.

This is the result of the anisotropy characterizing the latitudinal and

longitudinal error surfaces. The orientation of the estimated baseline is

not important for the propagation of the radial orbital errors because

the radial error surfaces are isotropic. Since all of the baselines of

Table 19 are extended in both north-south and east-west directions (Fig.

18), they will be affected by all three orbital errors (radial, latitudinal

and longitudinal}.

The estimable quantities, entering the satellite perturbations and

being neglected in the simplification process, implicitly assume zero

values in the resulting orbital model. If these constraints result in

radial errors, then depending on the length of the estimated baselines,

these errors should be relatively large in order to affect those baselines

beyond the centimeter level. This in turn would allow longer lengths of

continuous integration and, therefore, a reduced number of observations

would be necessary to achieve steady state response of the SRD method

(Section 4.4.2}. However, if these constraints result, in addition to

radial errors in both latitudinal and longitudinal errors, then even

smaller lengths of continuous integration may affect baselines of

moderate length (<2000 kin) beyond the centimeter level. To set up the

guidelines as to what simplifications of the orbital model can be applied

without affecting the accuracy of the estimated baseline beyond the

centimeter level, several tests have been performed the results of which

are shown in Table 19.

Table 19 contains the baseline differences obtained as the orbital

model was simplified from one containing a 12x12 gravity field, the

direct point mass (PM) effects of Sun and Moon, the tidal(TD) effects

due to Sun and Moon, the solar radiation (SR) pressure effects and the

along-track (AT) acceleration effects, to that containing only a 2x2

gravity field and the direct (PM) effects of the Sun and the Moon. The

first column of Table 19 lists the orbital models employed to estimate the

baselines which subsequently were differenced from those estimated on

the basis of orbital model shown on the title of this table {i.e.,

12x12+(I)).

166

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

l

l

Table 19 Baseline Differences (in meters) With Respect to Those

Computed Using an Orbital Model Including a 12x12

Gravity Field + (I)

Force Model Gravity 7110-7265 7109-7110 7110-7122

Gravity Field + ( ) (274069.50) (883602.25) (1437139.30)

I

I

I

I

I

12x12 + (2)

12x12 + (3)

12x12 + (4)

12x12

lOxlO + (2)

lOxlO + (4)

8x8 + (2)

8x8 + (4)

6x6 + (2)

6.6 + (4)

4.4 + (2)

4x4 + (4)

3x3 + (4)

2x2 + (4)

0.01 0.00 0.00

0.02 0.00 0.04

-0.08 0.00 0.10

5.14 -0.02 -3.15

0.01 0.00 0.01

-- 0.00 --

0.03 0.00 0.01

-- 0.00 --

1.02 0.00 2.02

-- 0.01 --

0.90 0.00 3.37

-- -0.01 --

-- -0.02 --

-- -0.07 --

I

I

I

I

I

I

I

I

I

I

(I) (PM) + (TD) + (SR) + (AT)

(2) (PM) + (TD) + (SR)

(3) (PM) + (TD)

(4) (I'M)

PM = point mass effects of sun & moon

TD = tidal effects due to sun & moon

SR = solar radiation pressure effects

AT = along-track acceleration effects

The resulting differences are shown in the corresponding rows for only

three baselines estimated on the basis of arcs up to seven days (Table

II, row 4), up to one hour (Table 12, row 4) and up to three days

(Table 15, row 4).

Elimination of the AT acceleration, the SR pressure, the TD

acceleration and the direct PM effects of the Sun and the Moon,

introduces into the resulting orbital errors secular, long-period and

short-period terms. Short-period terms directly related to Lageos' mean

anomaly are common in all of these errors.

Since the eccentricity of Lageos' orbit is very small (~0.0039),

elimination of the AT acceleration results primarily in latitudinal and

longitudinal orbital errors (Obenson, 1970). The long periods of the

resulting errors range from 66 days to 1100 days (Smith et al., 1985).

Elimination of the SR pressure acceleration results also in radial,

latitudinal and longitudinal orbital errors. Besides the additional

167
short-period terms related to the "ramp-like" behavior of the solar

radiation pressure, the long periods of the resulting orbital errors are

directly related to the motion of the Earth around the Sun and to the

motion of the ecliptic in space (Musen, 1960; Kozai, 1961; El'Yasberg,

1967; Blitzer, 1970). In addition to the short-period perturbations

related to the daily motion of the Earth, the nonstationary disturbances

of the Earth's potential due to the attraction of the Sun and the Moon

perturb the satellite periods with periods greater than a week (Kozai,

1973; Goad, 1977). As a result radial, latitudinal and longitudinal orbital

errors, having those periods will be introduced from the elimination of

the TD acceleration. Finally, elimination of the direct PM effects of the

Sun and the Moon introduces again the three types of orbital errors

having, besides the short periods associated with Lageos' mean anomaly,

intermediate periods associated with the motion of the Earth around the

Sun( ~ multiple of 180 days) and the motion of the Moon around the

Earth( ~ multiple of 14 days), and long periods associated with the space

motions of the ecliptic and the orbital plane of the Moon (Kozai, 1959;

El'Yasberg, 1967; Fisher, 1971; Blitzer, 1970).

Therefore, elimination of each of the mentioned accelerations

introduces radial, latitudinal and longitudinal orbital errors. Since

periods longer than a week are present, these errors will accumulate as

the integration periods increase from one hour to seven days. The

effects of these errors on the estimated baselines are shown in Table 19

(rows 1-4).

The perturbations caused in the motion of the satellites by the

ocean tides can reach as much as 20% of the perturbations caused by

the tides of solid Earth (Musen, 1973). Inspection of Table 19, row 3,

reveals that elimination of ocean tidal effects from the orbital model, as

is the case in this study (Section 2.2.5), will hardly affect the estimated

baselines at the centimeter level.

Therefore, centimeter level accuracy for baselines up to 1500 kin,

estimated via long-arc solutions, can be achieved if in addition to

gravitational effects of the Earth, the SR pressure, the TD and the

direct PM effects of Sun and Moon are included in the orbital model.

168

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

l

I

I

I

l

l

I

However, if steady state response is possible via a short arc solution,

centimeter level accuracy for baselines up to 1000 km can be achieved if

in addition to the gravitational effects of the Earth only the direct PM

effects of the Sun and the Moon are included (Table 19, column 3).

The question as to how much the gravity field of the Earth can be

reduced without affecting at the centimeter level the estimated baselines

is investigated by performing several solutions, the results of which are

shown Table 19 (rows 5-14). In these solutions the gravity field is

being reduced by two degrees at a time while the SR pressure, the TD

and the direct PM effects of the Sun and the Moon are included in the

orbital model (Table 19, rows 5, 7, 9 and 11). For baseline 7109-7]]0,

estimated via short-arc solutions, an additional solution has been

performed with each reduced gravity field and with only the direct PM

effects of the Sun and the Moon (Table 19, column 3, rows 6, 8, 10, 12,

13 and 14}. Eliminating two degrees at the time introduces secular

orbital errors due to elimination of the even zonal harmonics,

long-period errors due to the elimination of the zonals, and short

periods due to the elimination of all of the harmonics included in those

degrees. Thus, as the length of continuous integration increases radial,

latitudinal and longitudinal orbital errors tend to increase to the extent

that the estimated baselines will be affected beyond the centimeter level.

Careful study of the results shown in Table 19 reveals that

baselines of up to 1500 km estimated via the SRD method will be affected

at just the centimeter level if the orbital model includes :

1) in short arc solutions:

gravity field 4x4, and

the direct point mass effects of the Sun and the Moon.

2) in long arc solutions with arcs up to three days:

gravity field (8x8},

the PM effects of the Sun and Moon,

TD effects, and

SR pressure effects.

169
3) in long-arc solutions with arcs up to seven days:

gravity field (10x10),

the PM effects of the Sun and Moon,

TD effects, and

SR pressure effects.

Therefore, the sophistication of the orbital model employed in the

present study results in baselines having centimeter-level accuracy.

This accuracy is compatible with the accuracy of the laser observations

employed in this study (Section, 4.2). This constitutes the answer to

question (1); question (2) has already been answered by the conclusions

stated above.

170

I

I

I

I

I

I

I

I

i

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

l

I

l

I

I

I

I

I

I

Chapter 5

CONCLUSIONS AND RECOMMENDATIONS

5.1 CONCLUSIONS

The severe requirements of the geometric solutions to have simultaneous

observations from four or more stations (Section 2.1.3}, results in

rejection of a large number of nonsimultaneous observations. The

rejection is even greater with weather dependent satellite observations

not specifically designed for simultaneous tracking as happens with the

laser range observations.

Although during the MERIT Main Campaign many stations committed

themselves to collecting simultaneous observations {Section 3.4), these

observations were not enough to yield strong geometric solutions

because of their sensitivity to configurations being close to those

leading to singularities B and C (Sections 2.1.3 and 4.3.1). Because of

these configurations, a steady state response in minimum constraint

geometric solutions was possible only for two very short baselines.

However, stronger implementation of the scale through an overconstraint

solution (Table 8, solution C3) resulted in a steady state response for

baselines of up to 2280 km (Section 4.3.2}. These baselines formed the

standards of comparison in assessing the accuracy of those obtained via

both the SRD semidynamic and range dynamic mode methods (Section

4.5). Although the scale in this overconstraint solution was partly

implied by that of the {CSR)85L01 (Section 4.3.2), it is feasible in the

forseable future to incorporate the VLBI scale on the basis of only one

or two baselines. This in turn could lead through a geometric solution

to the estimation of a large number of baselines of compatible accuracy.

These baselines could be effectively used to assess the accuracy of

those obtained via the range dynamic methods. This practice might be of

171
great importance with the upcoming millimeter accuracy laser range

systems.

As a result, regular baseline estimation through geometric solutions

using laser range observations is impossible and therefore, since

geometric solutions are free of errors affecting the orbit and the TRF

frame they should be used, whenever possible, as standards of

comparison.

In contrast to the geometric method, the SRD semidynamic mode

method can be effectively used for regional baseline estimation,

especially with laser range observations to Lageos, since the altitude of

its orbit makes it possible for stations as far as 3703 km apart to collect

enough simultaneous observations to result in a steady state response of

the SRD method. The number of observations required for such a

response to take place is a function of

o

the baseline length,

o

the geometry and the length of the arcs employed,

g

and the accuracy of the available observations.

Table 20 shows the number of passes used in the present study to

achieve steady state response for different baseline lengths. In

general, as the length of the estimated baseline increases, more arcs,

stronger geometry, and more observations are required for a steady

state response to take place. The speed of the steady state response

will primarily depend on the number, distribution and accuracy of the

available observations, and on the length of the arcs employed in t.he

SRD solutions.

For instance, when the geometry is favorable (i.e., parallel passes

are available) and the observations of the baseline end stations have an

accuracy of about 3 cm, baselines of 883 km can be estimated via both

short arc and long arc solutions (Table 12). This is possible on the

basis of only 12 passes collected within a period of a week. However,

as the accuracy of the observations decreases, even in the light of

favorable geometry, the steady state response of a baseline having

172

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

l

I

I

I

I

I

I

I

I

Table 20 Steady State Response of the SRD Method

No. of Passes Approx. Occupation Time Steady State

Needed in This Dataset Up to __ km

I0 - 15 1 week I000

20 - 25 3 months 1500

25 - 30 3 months 2500

50 - 55 8 months 3500

about the same length as that in Table 12 (see Table 13) was not

possible on the basis of approximately the same number of observations

and through a short arc solution because the accuracy of its end

stations was worse than that of Table 12 (see Table 2). Nevertheless, a

steady state response for this baseline was possible through a solution

employing arcs of up to two days long. For longer baselines of up to

1437 km, the end stations of which recorded observations of 3 cm and 9

cm accuracy respectively, a steady state response in the light of

favorable geometry was possible on the basis of 22 passes and only

through a long arc solution employing arcs up to three days long (Table

15).

When the baseline-pass geometry is not favorable, (i.e.,passes *30"

to *50") and the accuracy of the observations of the baseline end

stations is approximately 3 cm and 8 cm respectively, a steady state

response for a 274 km baseline was possible on the basis of 17 passes

and only through a long arc solution employing arcs up to four days

long (Table 11). For another baseline of about 1198 km having

unfavorable pass-baseline geometry, a steady state response was

possible on the basis of 40 passes and through a solution employing

arcs up to seven days long (Table 14). The number of observations

required for the steady state response of this baseline is substantially

reduced as compared to those required for all the other baselines shown

in Tables 9 through 17. This is not surprising because the stronger

geometry implied by the 40 passes compensated for the reduced number

of observations. However, if more observations could have been

collected during each pass, fewer passes would be required for a steady

173
state response.

For very long baselines of about 3703 km, although passes of

favorable geometry are not available, a steady state response was still

possible on the basis of 54 passes and through a solution employing

arcs up to one day long (Table 16). This response, however, started

diverging when for the determination of this baseline arcs up to seven

days long were employed (Table 16). Furthermore, for a baseline of

about 3559 kin, the pass-baseline geometry of which is very weak (i.e.,

passes within t20" to ±60°), a steady state response was possible only

with arcs up to three days (Table 17). Therefore, for so long baselines

and when the geometry is fairly weak care should be exercised as to

whether or not a steady state response is possible.

From the above discussion it is evident that for baselines up to

1500 kin, steady state response can be reached without any problems

even if the geometry is not favorable. However, tracking passes of

favorable geometry with accurate (i.e., 3 cm) and high repeatability laser

instruments increases the resolution of baseline recovery to weekly

estimates. Although this was shown only for baselines of up to 883 kin,

baselines of up to 1437 km can also be recovered on a weekly basis,

because steady state response for these baselines was reached on the

basis of about 69,000 observations (Table 15).

In addition to the increased temporal resolution the comparison of

the SRD method with the geometric solution (Table 8, solution C3)

evidently demonstrated an accuracy of lcm (Table 18), accomplished on

the basis of a simple modeling and a limited orbit adjustment (Section

4.6). This demonstrates the potential of the SRD method for accurate

differential positioning, and the insensitivity of the SRD observable not

only to the errors affecting the orbit but also to those affecting the

implementation of the TRF frame. The simple modeling of the orbit and

its limited adjustment result in relaxing requirements in regard to

availability of sophisticated software and extensive computer facilities

which in turn makes it feasible to implement this method in the PC

environment.

174

I

I

I

I

I

I

I

i

I

I

i

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

The stability characteristics of the SRD single baseline semidynamic

solutions make it possible to estimate baselines of up to 1437 km without

any need for a local support tracking network which is required in

regional baseline estimation via the range semidynamic mode method

(Christodoulidis et al., 1981).

The steady state response of baseline 7109-7105 (Table 16}, reached

on the basis of weak pass-baseline geometry strongly suggests that with

parallel passes single SRD baseline solutions may very well lead to

accuracy at the 2 to 3 cm level for baselines of up to 3500 km or 4000

km.

The stability characteristics and the high accuracy of the SRD

method make it ideal for regional baseline estimation. This in turn

makes it possible to obtain baselines up to 1500 km on a weekly basis

free of the fluctuations usually affecting the monthly estimates of the

range dynamic method. The stability characteristics could be greatly

enhanced with a network solution which on the basis of SRD normal

points would also substantially reduce the bulk of the computations

without affecting the speed of the steady state response since the

almost noise free SRD normal points will lead to a faster steady state

(Section 4.4). However, with a network solution there might exist.

baselines for which simultaneous observations will not be available or"

the SRD observables will contain very little information for those

baselines. This happens when the intersection of the visibility regions

of the baseline end points is either an empty set or one having very

few points. If a steady state response for those baselines has been

reached it will entirely depend on the strength of the network which in

the semidynamic environment tends to depend on the geometric strength

of the observations which may not lead to a steady state response for

longer baselines as it was the case with the geometric solutions. Steady

state response of those baselines would require increasing the

integration period. If this increase leads to a steady state, the

estimated lengths for those baselines would primarily depend on the

orbital strength, and therefore they would be susceptible not only to

accumulated orbital errors but also to any inconsistencies affecting the

175
implementation of the TRF frame. However, employing the same orbital

model as in the range dynamic solutions, the SRD baseline estimates will

be less affected by the accumulated orbital errors.

In a network solution care should be exercised to account for the

correlations resulting from the redundancy of having pairs of stations

with common observations.

In the implementation of the SRD method a great deal of effort was

spent to edit the data and create the SRD observables. Editing the

observed laser ranges before generating those observables is crucial

because generating them on the basis of erroneous laser ranges leads to

the distribution of these errors among all of the generated SRD

observables and therefore it will be difficult, if not impossible, to

effectively edit the SRD observables in the final orbit adjustment.

The effects of gaps should also be well controlled and should be

kept below the noise level of the observations. This was the reason

why so much care was taken in the generation of SRD observables

(Chapter 3). This in turn increased considerably the bulk of the

computations required for the implementation of the SRD method.

These two factors do not by any means limit the potential of the

SRD method because through proper arrangements either full rate SRD

observables or preferably SRD normal points could be very easily

generated with some slight modifications of the software employed by

the computational centers responsible for the generation of the range

normal points. With SRD normal points at our disposal a network SRD

solution could be performed even in a PC environment, and since normal

points are almost noise free they would lead to a faster steady state

response of the SRD method.

5.2 RECOMMENDATIONS

Since the SRD method, based on a relatively simple orbital modeling and

a limited orbit adjustment, is for regional baseline estimations at least as

accurate as the dynamic mode methods, it is recommended that for such

applications serious consideration should be given to regularly implement

176

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

i

I

I

I

I

I

I

I

this method for accurate differential positioning. If such a step is to

be taken it is also suggested that:

° Further effort should be made by the observing stations to

achieve simultaneous tracking.

° The direction of the baselines being estimated should be chosen, if

this is feasible, to closely resemble the two main Lageos groundtrack

directions.

SRD normal points should be generated by the computational center

responsible for the generation of range normal points.

Research the response of the method, on the basis of normal points,

when different network configurations are employed.

Projects specifically designed for studies of regional plate tectonic

motions, determined on the basis of laser range observations to Lageos,

are ideal for implementation of the SRD method, such a project is

currently under way in the area around the Mediterranean Sea

(MEDLAS/WEGENER project), and therefore for this project the SRD

method could offer an accurate and inexpensive alternative to study

regional plate tectonic motions.

177
Appendix A

PARTIAL DERIVATIVES OF THE SRD OBSERVABLE

Expressing the SRD observable relative to an earth-fixed frame takes

the following form:

(A-l)

with

Ej2 : Xj - :2 and Ej, = Xj - :, (A-2)

where X j, Y2 and Y, are the earth-fixed position vectors of the satellite

at the epoch j and stations (2) and (1) respectively.

Differentiating 6pj with respect to Y2 one obtains:

a6pj a6pj aEj2 1 Ej2 T

_--=-:= - _j_T.(-I)-

aY2 Ej2 aY2 Pj2 Pj2

Thus

a6pj (Xj-Y2) T

a:2 pj2

(A-3)

where pj2 is the distance from station (2) to satellite position J.

Suppose that

Ro = (Xo, Yo, Zo) (A-4)

and

-y

Ro : (Xo, 40, _o) T (A-5)

are the initial position and velocity vectors with respect to an inertial

frame. Differentiating equation (A-l) with respect to the initial state

7"

vector (Ro, Ro) one obtains:

178

!

!

I

I

I

I

I

I

I

i

I

I

I

I

I

I

I
I

I

i

I

I

I

I

I

I

I

I

I

I

I

I

i

I

I

a6pj a6pj aEja aXj a6pj aEj_

a(_o,ho) a_.j_a_j a(Ro,_o) aEj, af(j (A-6)

From equ. (A-l) and (A-2) it follows that:

a6pj 1 _ XjT-y2T

a'_j2 = PJ2 EJ2T = PJ2

(A-7)

CA-B)

_(_o,_o)-,_o'""'_To,: , _o '"" _ ,

(A-9)

where p j2 is the distance from the station (2) to the satellite position j;

Rj is the inertial position vector of the satellite at the epoch j; and S,

N, P are the earth rotation, nutation and precession matrices

respectively. Similarly to eq. (A-7), the partial a6pj/aEjl takes the

following form:

.a6pj XjT-y1T

aEjl pjz

(A-IO)

where p j, is the distance from station (1) to satellite position j.

Substituting eq. (A-7) - (A-10) into eq. (A-6) the following formula

results:

a6pj {(Xj-Y2) T

a(Ro,Ro) - pj2 - - } [a_j a_j]

(Xj-Yz) T

pjl "l_o ' "'" _-_oj (A-II)

or

a6pj {(Xj-Y2) T

a(_o,_o)- pj_ - - ] [ a_j ahjl (A-12)

Thus, from eq. {A-11) and {A-12) it is easily deduced that the partial

derivative of the SRD observable with respect to the first coordinate

7"

(i.e., Xo) of the initial state vector (Ro,Ro) takes the following form:

179
_-K-.- ,_ ;7:" _×--7 <A-14)

Similarly one can obtain the partial derivative of the SRD observable

with respect to all of the remaining componenets of the inertial initial

state vector.

180

I

I

i

I

I

I

I

I

I

I

I

I

I

I

I

i

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

Appendix B

SENSITIVITY OF THE PARTIAL DERIVATIVES WITH RESPECT TO THE

INITIAL STATE VECTORS

With the tendency of the observed ranges pj2 and pit from stations (1)

and (2) to become equal (Section 4.4.1).

pj= = pit = pj

equation (A-11) becomes

a(_o,fio) pj

(B-l)

(B-2)

This equation can also be expressed as follows

a6pj 1

a<_o,fio) pJ

where (.) and (5_o , .... 5_oj) denote the dot product and the partial

derivative vectors taXo' "'" a_o" at the epoch 3 respectively.

If _xJo,, ..., _zJo designate the _ngles between the baseline vector

and the partial derivative vectors respectively, equation (B-3) takes the

following form

a6pj 1

a(Ro,_o) PJ IY_-Y=l(IDxjl.cos(ox_), ..., ID_lcoso_j) (B-4)

Therefore, from epoch to epoch the variations among the entries of

the columns, corresponding to the components of an initial state vector,

are controlled by the variations of the projection of the partial derivative vectors -Dx_, ..., DzJo in the direction of the estimated baseline.

These variations tend to be reduced as the observed part of the pass

181
becomes short and/or perpendicular to the estimated baseline. This

geometry is closer to reality for longer baselines because as the

baseline length increases the intersection of the visibility regions of the

baseline endpoints tends to become smaller and closer to the

perpendicular bisector plane of the baseline, thereby leading to a

geometry close to also fulfilling assumption (B-l).

182

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

!

i

I

REFERENCES

Abshire, J.B., C.S. Gardner (1985), "Atmospheric Refractivity Corrections

in Satellite Laser Ranging," IEEE Transactions on Geoscience and Remote

Sensing, Vol. GE-23, No. 4.

Abshire, J.B., T.W. Zagwodzki, J.F. McGarry, J.J. Degnan (1984), "An

Experimental Large Aperture Satellite Laser Ranging Station at GSFC,"

Proc. of Fifth International Workshop on Laser Ranging Instrumentation,

Herstmonceux Castle, U.K., Sept. 10-14; J. Gaignebet, ed., Groupe de

Recherches de Geodesie Spatiale.

Afonso, G., F. Barlier, C. Berger, F. Mignard, J.J. Walch (1985),

"Reassessment of the Charge and Neutral Drag of Lageos and Its

Geophysical Implications," J. of Geophys. Res., Vol. 90, No. Bll,

9381-9398.

"Analysis of Lageos Laser Range Data, Sept. 1983 - Oct. 1984," Center

for Space Research, Univ. of Texas at Austin.

"Analysis of Lageos Laser Range Data, Aug. 1985," Center for Space

Research, Univ. of Texas at Austin.

The Astronom/ca] Almanac, 1984, U.S. Govt. Printing Office, Washington,

D.C.

Baarda, W. (1968), "A Testing Procedure for Use in Geodetic Networks,"

Publications on Geodesy, New Series, Vol. 2, No. 5, Netherlands Geodetic

Commission, Delft.

Bjerhammar, A. (1985), "On a Relativistic Geodesy," Bull. Geodesique, Vol.

59, No. 3, 207-220.

Blahs, G. (1971 ), "Investigations of Critical Configurations for

Fundamental Range Networks," Dept. of Geodetic Science Rep. 150, Ohio

State Univ., Columbus.

Blitzer, L. (1970), Handbook of Orbital Perturbations, TRW Systems

Group, Redondo Beach, Calif.

Boucher, C., M. Feissel (1984), "Realization of the BIH Terrestrial

System," Proc. of the International Symp. on Space Techniques for

Geodynamics, Research Inst. of the Hungarian Academy of Sciences,

Sopron, Hungary.

183
Boucher, C., Z. Altamimi (1985), "Towards an Improved Realization of the

BIH Terrestrial Frame," Proc. of the International Conf. on Earth

Rotation and the Terrestrial Reference Frame, Vol. 2, Dept. of Geodetic

Science and Surveying, Ohio State Univ., Columbus.

Boucher, C., Z. Altamimi (1986), "Status of the Realization of the BIH

Terrestrial System," Proc. of IAU Symp. 128, The Earthts Rotation and

Reference Frames for Geodesy and Geodynamics, Reidel Publ.

Bureau International de l'Heure, Annual Reports, 1983, 1984.

Cappellari, J.O., C.E. Velez, A.J. Fuchs, eds. (1976), "Mathematical Theory

of the Goddard Trajectory Determination System," NASA/Goddard Space

Flight Center, NASA X-582-76-77, NTIS N76-24291 - N76-24302.

Carnahan, B., H.A. Luther, J.O. Wilkes (1969), Applied Numerical Methods,

Wiley and Sons, New York.

Christodoulidis, D.C., D.E. Smith (1981), "Prospects for TLRS Baseline

Accuracies in the Western USA Using Lageos," NASA Technical

Memorandum 82133, Goddard Space Flight Center, Greenbelt, Maryland.

Christodoulidis, D.C., D.E. Smith, S.M. Klosko, M.H. Torrence, P.J. Dunn

(1985), "A GSFC Alternative to the SLR MERIT Standards," Proc. of the

International Conf. on Earth Rotation and the Terrestrial Reference

Frame, Vol. 2, Dept. of Geodetic Science and Surveying, Ohio State Univ.,

Columbus.

Coates, R.J., H. Frey, G.D. Mead, J.M. Bosworth (1985), "Space-Age

Geodesy, The NASA Crustal Dynamics Project," IEEE Transactions on

Geoscience and Remote Sensing, Vol. GE-23, No. 4, 360-367.

Cohen, S.C., D.E. Smith (1985), "Lageos Scientific Results: Introduction,"

J. of Geophys. Res., Vol. 90, No. BII.

Davis, P.J. (1975), Interpolation and ApproximaLion, Dover Publications,

New York.

Degnan, J.J. (1985), "Satellite Laser Ranging: Current Status and Future

Prospects," IEEE Transactions on Geoscience and Remote Sensing, Vol.

GE-23, No. 4, 398-413.

DeGroot, M.H. (1975), Probability and Statistics, Addison-Wesley Publ.,

Reading, Massachusetts.

Diamante, J., R. Williamson (1972), "Error Models for Solid Earth and

Ocean Tidal Effects in Satellite Systems Analysis," Wolf Research and

Development Corp., Greenbelt, Maryland

Drain, L.E. (1980), The Laser Doppler Technique, Wiley and Sons, New

York.

184

I

I

I

l

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

Eanes R.J., B. Schutz, B. Tapley {1983), "Earth and Ocean Tide Effects

on Lageos and Starlette," Proc. of the Ninth International Syrup. on

Earth Tides, E. Schweizerbart'sche Verlagsbuchhandlung, Stuttgart, FRG.

El'Yasberg, P.E. (1967), "Introduction to the Theory of Flight of Artificial

Earth Satellites," translated from Russian by Z. Lerman, Israel Program

for Scientific Translations, Jerusalem. Avail. Clearinghouse for Federal

Scientific and Technical Information, Springfield, VA 22151.

Fisher, D. (1971), "Lunisolar Perturbations of the Motion of Artificial

Satellites," Mission Trajectory Determination Branch, Goddard Space

Flight Center, Greenbelt, MD, Preprint X-552-71-60.

Fitzmaurice, M.W., P.O. Minott, J.B. Abshire, H.E. Rowe {1977), "Prelaunch

Testing of the Laser Geodynamic Satellite (LAGEOS)," NASA Tech. Paper

1062, NASA/Goddard Space Flight Center, Greenbelt, Maryland.

Forsythe, G., C.B. Moler (1967), Computer Solution of Linear Algebraic

Systems, Prentice-Hall, Englewood Cliffs, New Jersey.

Fukushima, T. (1986), "Coordinate System in General Relativity," Proc.

of IAU Syrup. 128, The F_rth's Rotation and Reference Frames for

Geodesy and Geodynamics, Reidel Publ.

Goad C.C. (1977), "Application of Digital Filtering to Satellite Geodesy,"

NOAA Technical Rept. NOS 71 NGS 6, Rockville, Maryland.

Grafarend, E., B. Schaffrin (1974), "Unbiased Free Net Adjustment,"

Survey Review, Yol. XXII, No. 171, 200-218.

Gruen, A. (1979), "Gross Error Detection in Bundle Adjustment," pres. at

International Syrup. on Current Trends in the Development, Execution

and Use of Aerial Triangulation, Brisbane, Australia, Oct. 15-17.

Halliday, D., R. Resnick (1962), Physics, Part TT, Wiley and Sons, New

York.

Hamilton, W.C. (1964), S_atistics in Physical Science, Ronald Press, New

York.

Heck, B. (1982), "Report of the FIG Working Group on the Analysis of

Deformation Measurements," FIG Syrup., Budapest.

Heiskanen, W.A., H. Moritz (1967), Physical Geodesy, W.H. Freeman, San

Francisco.

Himwich, W.E., E.J. Harder (1986), "Direct Estimation of Nutation

Coefficients from VLBI Data," Proc. of IAU Syrup. 128, The Earthts

Rotation and Reference Frames for Geodesy and Geodynamics, Reidel

Publ.

185
Johnson, T.S., J.J. Degnan, T.E. McGunigal (1978), "200 Picosecond Laser

Development. A Status Report," Workshop on Laser Ranging

Instrumentation, III, NASA TM-79643, NASA/Goddard Space Flight Center,

Greenbelt, MD.

Kaula, W.M. (1961), "Analysis of Gravitational and Geometric Aspects of

Geodetic Utilization of Satellites," Geophysical J. of R. Astron. Soc., Vol.

5, No. 2.

Kaula, W.M. (1966), Theory of Satellite Geodesy, Applications of Satellites

to Geodesy, Blaisdell Publ., Waltham, Massachusetts.

Kovalevsky, J., I.I. Mueller (1981), "Comments on Conventional Terrestrial

and Quasi-Inertial Reference Systems," Reference Coordinate Systems for

Earth Dynamics, E.M. Gaposchkin and B. Kolaczek, eds., Reidel Publ.

Kozai, Y. (1959), "On the Effects of the Sun and the Moon Upon the

Motion of a Close Earth Satellite," Smithsonian Institution Astrophysical

Observ., Research in Space Science Special Rep. 22.

Kozai, Y. (1961), "Effects of Solar Radiation Pressure on the Motion of an

Artificial Satellite," Smithsonian Institution Astrophysical Observ.,

Research in Space Science Special Rep. 56.

Kozai, Y. (1973), "A New Method to Compute Lunisolar Perturbations in

Satellite Motions," Smithsonian Astrophysical Observ., Special Rep. 349.

Krakiwsky, E.J., A.J. Pope (1967), "Least Squares Adjustment of Satellite

Observations for Simultaneous Directions or Ranges," Dept. of Geodetic

Science Rep. 86, Ohio State Univ., Columbus.

Krakiwsky, E.J., B. Wanless, B. Buffett, K.P. Schwarz, M. Nakiboglu

(1985), "GPS Orbit Improvement and Precise Positioning," Proc. of the

First International Syrup. on Precise Positioning with the Global

Positioning System, Rockville, Maryland.

Krogh, F.T. (1969a), "A Variable Step Variable Order Multistep Method for

the Numerical Solution of Ordinary Differential Equations," Information

Processing, 68, North-Holland Publ., Co., Amsterdam.

Krogh, F.T. (1969b), "VODQ/SVDQ/DVDQ, Variable Order Integrators for

the Numerical Solutions of Ordinary Differential Equations," Technology

Utilization Document No. CP-2308, Jet Propulsion Laboratory, Pasadena,

California.

Krogh, F.T. (1973a), "On Testing a Subroutine for the Numerical

Integration of Ordinary Differential Equations," J. of the Assoc. for

Computing Machinery, Vol. 20, No. 4, 545-562.

Krogh, F.T. (1973b), "Algorithms for Changing the Step Size," SIAM J.

Numer. Anal., Vol. 10, No. 5, 949-965.

186

I

I

I

I

I

I

I

I

I

I

I

!

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

Krogh, F.T. (1974), "Changing Stepsize in the Integration of Differential

Equations Using Modified Divided Differences," Proc. of the Conf. on the

Numerical Solution of Ordinary Differential Equations, Lecture Notes in

Mathematics, Vol. 362, Springer-Verlag, New York, 22-71.

Landau, L.D., E.M. Lifshitz (1960), Mechanics, transl, from Russian by

J.B. Sykes and J.S. Bell, Pergamon Press.

Latimer, J.H., E.M. Gaposchkin (1977), "Scalar Translocation Using Laser

Range Data," Smithsonian Astrophysical Observ., Center for Astrophysics

Preprint Series No. 750, Cambridge, Massachusetts.

Lerch, F.J., S.M. Klosko, G.B. Patel, C.A. Wagner (1985), "A Gravity Model

for Crustal Dynamics (GEM-L2}," J. of Geophys. Res., Vol. 90, No. Bll,

9301-9311.

Lieske, J.H. (1979), "Precession Matrix Based on IAU (1976} System of

Astronomical Constants," Astronomy" and Astrophysics, 73, 282-284.

Marini, J.W., C.W. Murray, Jr. (1973), "Correction of Laser Range

Tracking Data for Atmospheric Refraction at Elevation Angles Above l0

Degrees," NASA Technical Rep. X-591-73-351, NASA/Goddard Space Flight

Center, Greenbelt, Maryland.

Melbourne W., R. Anderle, M. Feissel, R. King, D. McCarthy, D. Smith, B.

Tapley, R. Vicente (1983), "Project MERIT Standards," U.S. Naval

Observatory Circular No. 167, Washington, D.C. 20390.

Montag, H., G. Gendt, R. Dietrich, K. Kurth (1985), "Investigation of Polar

Motion and the Length of Day by Means of SLR Data of the MERIT

Campaign," Proc. of the International Conf. on Earth Rotation and the

Terrestrial Reference Frame, Dept. of Geodetic Science and Surveying,

Ohio State Univ., Columbus.

Moritz, H. (1979), "Concepts in Geodetic Reference Frames," Dept. of

Geodetic Science Rep. 294, Ohio State Univ., Columbus.

Moritz, H., I.I. Mueller (1987)9 Earth Rotation, Theory and Observation,

Ungar Publ., New York.

Moyer, T.D. (1971), "Mathematical Formulation of the Double-Precision

Orbit Determination Program (DPODP)," TR-32-1527j Jet Propusion

Laboratory, Pasadena, California.

Mueller, I.I. (1964a), "The Geodetic Applications of Satellites," Institute of

Geodesy, Photogrammetry and Cartography Rep. 34, Ohio State Univ.,

Columbus.

Mueller, I.I. (1964b), Introduc_on to Sate//ite Geodesy, Ungar Publ., New

York.

187
Mueller, I.I. (1969), Spherical and Practical Astronomy As Applied to

Geodesy, Ungar Publ., New York.

Mueller, I.I., B.H.W. van Gelder, M. Kumar (1975), "Error Analysis for the

Proposed Close Grid Geodynamic Satellite Measurement System

(CLOGEOS)," Dept. of Geodetic Science Rep. 230, Ohio State Univ.,

Columbus.

Mueller, I.I. (1981), "Reference Coordinate Systems for Earth Dynamics, A

Preview," Reference Coordinate Systems for Earth Dynamics, E.M.

Gaposchkin and B. Kolaczek, eds., Reidel Publ.

Mueller, I.I. (1985), "Reference Coordinate Systems and Frames:

Concepts and Realization," Bull. Geodesique, Vol. 59, No. 2, 181-188.

Mueller, I.I., G.A. Wilkins (1986), "Earth Rotation and the Terrestrial

Reference Frame, Adv. Space Res., Vol. 6, No. 9, 5-11.

Obenson, G. (1970), "Direct Evaluation of the Earth's Gravity Anomaly

Field from Orbital Analysis of Artificial Satellites," Dept. of Geodetic

Science Rep. 129, Ohio State Univ., Columbus.

Musen, P. (1960), "The Influence of the Solar Radiation Pressure on the

Motion of an Artificial Satellite," J. of Geophys. Res., Vol. 65, No. 5,

1391-1400.

Musen, P. (1973), "A Semi-Analytical Method of Computation of Oceanic

Tidal Perturbations in the Motion of Artificial Satellites, X-590-73-190,

NASA/Goddard Space Flight Center, Greenbelt, MD.

Pavlis, E.C. (1982), "On the Geodetic Applications of Simultaneous

Range-Differencing to LAGEOS," Dept. of Geodetic Science and Surveying

Rep. 338, Ohio State Univ., Columbus.

Pavlis, E.C., I.I. Mueller (1983), "The Effect of Earth Orientation Errors

in Baseline Determination," Bull. Geodesique, Vol. 57, No. 3, 273-282.

Pope, A.J. (1976), "The Statistics of Residuals and the Detection of

Outliers," NOAA Technical Rep. NOS 65 NGS1, National Geodetic Survey,

Rockville, Maryland.

Rao, C.R. (1973), Linear SLatistical Inference and Its Applications, 2nd

ed., Wiley and Sons, New York.

Rinner, K. (1966), "Systematic Investigations of Geodetic Networks in

Space," U.S. Army Research and Development Group (Europe), May.

Rizos, C., A. Stolz (1985), "Force Modeling for GPS Satellite Orbits," Proc.

of First International Syrup. on Precise Positioning with the Global

Positioning System, Rockville, Maryland.

188

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I
I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

l

Rubincam, D.P. (1980), "Atmospheric Drag As the Cause of the Secular

Decrease in the Semimajor Axis of LAGEOS' Orbit," Geophys. Res. Lett., 7,

468-470.

Rubincam, D.P., N.R. Weiss (1985), "The Orbit of LAGEOS and Solar

Eclipses," J. of Geophys. Res., Vol. 90, No. Bll, 9399-9402.

Schutz, B.E. (1983a), "Satellite Laser Ranging Participants During Project

MERIT," Center for Space Research, Univ. of Texas at Austin.

Schutz, B.E. (1983b), "Satellite Laser Ranging Procedures Guide for

Project MERIT," Center for Space Research, Univ. of Texas at Austin.

Schutz, B.E. (1985), "Report of the Satellite Laser Ranging Coordinator,"

Report on the MERIT-COTES Campaign on Earth Rotation and Reference

Systems, Part I, G.A. Wilkins, ed., Royal Greenwich Observ.,

Herstmonceux, U.K.

Schwarz, C.R. (1969), "The Use of Short Arc Orbital Constraints in the

Adjustment of Geodetic Satellite Data," Dept. of Geodetic Science Rep.

118, Ohio State Univ., Columbus.

Shawe, M.E., A.G. Adelman (1985), "Precision Laser Tracking for Global

and Polar Motion," IEEE Transactions on Geoscience and Remote Sensin$,

Vol. GE-23, No. 4.

Shelus, P.J. (1983) "Transportable Laser Ranging Station TLRS-I,"

International Coordination of Space Techniques for Geodesy and

Geodynamics/' CSTG Bulletin No. 5, 112-116, Dept. of Geodetic Science

and Surveying, Ohio State Univ., Columbus.

Silverberg, E.C. (1982), "The Feedback Calibration of the TLRS Ranging

System," Proc. of the 4th International Workshop on Laser Ranging

Instrumentation, Geodetic Inst., Univ. of Bonn.

Smith, D.E., D.C. Christodoulidis, R. Kolenkiewicz, P.J. Dunn, S.M. Klosko,

M.H. Torrence, S. Fricke, S. Blackwell (1985), "A Global Geodetic

Reference Frame from Lageos Ranging (SL5.1AP)," J. of Geoph:ys. Res.,

Vol. 90, No. Bll, 9221-9233.

Spath, H. (1974), Splines Algorithms for Curves and Surfaces, transl.

from German by W.D. Hoskins and H.W. Sager, Utilitas Mathematica Publ.,

Winnipeg, Canada.

Standish, E.M. (1981), "Letters," Astron. Astrophys., 101, L17.

Tapley, B.D., B.E. Schutz, R.J. Eanes (1985a), "Station Coordinates,

Baselines, and Earth Rotation from LAGEOS Laser Ranging: 1976-1984," J.

of Geophys. Res., Vol. 90, No. BII.

189
Tapley, B.D., R.J. Eanes. B.E. Schutz (1985b), "UT/CSR Analysis of Earth

Rotation from Lageos SLR Data," Proc. of the International Conf. on

Earth Rotation and the Terrestrial Reference Frame, Vol. 1, Dept. of

Geodetic Science and Surveying, Ohio State Univ., Columbus.

"Transportable Laser Ranging System No. 2 (TLRS-2), Tech. Manual,"

NASA/GSFC Publ. MS-428-10, Goddard Space Flight Center, Greenbelt,

Maryland.

Tsimis, E. (1972), "On the Geometric Analysis and Adjustment of Optical

Satellite Observations," Dept. of Geodetic Science Rep. 185, Ohio State

Univ., Columbus.

Tsimis, E. {1973}, "Critical Configurations (Determinantal Loci) for Range."

and Range-Difference Satellite Networks," Dept. of Geodetic Science Rep.

191, Ohio State Univ., Columbus.

Uotila, U.A. {1967}, "Introduction to Adjustment Computations with

Matrices," lecture notes, Dept. of Geodetic Science, Ohio State Univ.,

Columbus.

Uotila, U.A. (1986}, "Notes on Adjustment Computations, Part I," Dept. of

Geodetic Science and Surveying, Ohio State Univ., Columbus.

Van Gelder, B.H.W. (1978), "Estimability and Simple Dynamical Analysis of

Range (Range-Rate and Range-Difference) Observations to Artificial

Satellites," Dept. of Geodetic Science Rep. 284, Ohio State Univ.,

Columbus.

Veis, G. (1960), "Geodetic Uses of Artificial Satellites," Smithsonian

Contrib. _o Astrophys., Vol. 3, No. 9.

Wahr, J.M. {1981a), "The Forced Nutations of an Elliptical, Rotating,

Elastic and Oceanless Earth," Geophys. J. R. astr. Sot., 64, 705-727.

Wahr, J.M. (1981b}, "Body Tides on an Elliptical, Rotating, Elastic

Oceanless Earth," C_ophys. J.R. astr. Sot., 64, 677-703.

Wilkins, G.A., ed. (1980), "A Review of the Techniques to Be Used During

Project MERIT to Monitor the Rotation of the Earth," publ. jointly by

Royal Greenwich Obs., Herstmonceux, U.K., and Inst. f. Angewandte

Geod_isie, Frankfurt, FRG.

Wilkins, G.A., I.I. Mueller (1986), "On the Rotation of the Earth and the

Terrestrial Reference System," Bull. Geodesique, Vol. 60, No. 1.

190

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I