MAXIMUM LIKELIHOOD ESTIMATION The maximum

European Geosciences Union General Assembly 2013 Vienna | Austria | 07 – 12 April 2013

NOISE CHARACTERISTIC IN GPS SUB-DIURNAL COORDINATES: MAXIMUM LIKELIHOOD ESTIMATION APPROACH 1)

1)

1)

Anna Kłos , Janusz Bogusz , Mariusz Figurski , Wiesław Kosek

CGS

1)

Military University of Technology, Centre of Applied Geomatics, Warsaw, Poland 2) Space Research Centre, Polish Academy of Sciences 3) Environmental Engineering and Land Surveying, University of Krakow e-mail: [email protected]

ABSTRACT

CATS SOFTWARE WITH MLE METHOD

The aim of this study was to evaluate the noises appearing in the changes of coordinates of Polish permanent stations belonging to EPN (EUREF Permanent Network). The analysis was made on the topocentric coordinates, but calculated not in the standard daily, but 3-hour solutions to determine changes in the sub-diurnal frequencies. The data covered the period of nearly 2 years from 2008 to 2010. The time series are characterized by the trend, but seasonal, cyclic and random variations as well. In the analyzed data the systematic component, which interests us the most and the random noise, which is a kind of disturbance that does not bring any significant information were isolated using CATS (Create and Analyze Time Series) software. It allows fitting the multi-parameter models in time series and analyzing the noise component with the use of the maximum likelihood estimation. The procedure which is used ensures the simultaneous solutions' searching for all of parameters. As the result of the calculations, the chosen software provides the noise parameters as well as their errors. The noise parameters include: the spectrum index value of coloured noise, the amplitudes of both white and coloured noises that appear in data, the amplitudes of sinuses and co-sinuses corresponding to determined frequencies and the plot's slope. The noise model occurring in the geodetic time series was presented for each of the topocentric coordinates: North, East and Up on each of considered stations and compared with the white noise model. The presentation also deals with the character of the noise existing in the sub-daily GPS solutions.

CATS (Create and Analyze Time Series) is the software that enables fitting of multiparameter models into considered time series and analyzing the noise component with the use of maximum likelihood estimation. Programme allows the analyzing of white and power-law noise parameters including flicker and random walk noise models. It also provides the option for combining the above models into one. CATS uses the least squares method to fit the mutliparameter model into time series and estimates the type and the amplitude of existing noise. The estimation process is based on determination of λ parameter that describes the probability distribution of xi data. When the distribution is known, we can define the probability of occurrence of individual xi. The values of parameters λ are chosen as to maximize the probability of occurrence of any xi. The estimation of noise components is realized by maximization of likelihood function of their occurrence l for the given data set x.

NOISES Noise in GPS time series is defined as power-law process with the power spectrum of the form: P0 and f0 are normalizing constants and f is temporal frequency.

æ f ö Px ( f ) = P0 çç ÷÷ è f0 ø

2),3)

RESULTS:

l ( x, C ) =

BOR1 station N component

stations BOGI BOGO BOR1 GWWL KRAW LAMA REDZ USDL WROC ZYWI

1 T -1 ˆ exp( 0 . 5 v C vˆ) 1/ 2 N /2 (2p ) (det C )

MLE function -32214,20 -26368,50 -23457,49 -23911,68 -25769,43 -23897,03 -25117,90 -36167,52 -21944,37 -28649,18

Intercept -19,50 -20,19 -19,41 -17,43 -18,04 -19,19 -19,37 -19,69 -20,13 -21,74

N component Slope White amplitude 13,39 x 13,72 x 13,38 x 12,43 x 12,77 x 13,18 x 13,23 x 12,36 3,83 13,61 x 14,70 x

E component

κ -0,44 -0,80 -0,59 -0,97 -0,77 -0,59 -0,68 -1,52 -0,90 -0,80

PL amplitude 5,53 7,34 3,53 8,63 6,50 3,65 5,00 65,44 6,15 9,11


MLE function -36395,19 -30034,52 -26877,66 -27766,07 -28642,97 -27908,88 -31252,60 -50610,56 -24801,07 -32093,45

Intercept -27,70 -28,55 -27,39 -26,33 -28,06 -27,21 -25,96 -27,85 -28,21 -28,86

E component Slope White amplitude 19,30 x 19,88 x 19,02 x 17,83 x 19,60 x 18,79 x 17,81 x 25,56 17,33 19,30 x 20,00 x

U component

κ -0,39 -0,37 -0,27 -0,90 -0,36 -0,29 -0,80 -1,71 -0,53 -0,59

PL amplitude 7,28 3,88 2,35 10,49 3,31 2,65 11,40 281,33 3,51 7,72


MLE function -41981,32 -39041,75 -33729,27 -35695,44 -37324,15 -35855,60 -38334,46 -41099,51 -34456,17 -40036,30

Intercept -2,99 -0,70 -3,99 7,87 0,28 -1,90 -4,48 -11,03 -3,02 -0,15

U component Slope White amplitude 1,81 x 0,20 x 2,50 x -5,96 x 0,01 x 0,88 x 3,39 x 4,81 x 1,80 x -0,32 x

κ -0,68 -0,71 -0,82 -0,98 -0,86 -0,72 -0,91 -1,43 -1,08 -0,92

PL amplitude 23,06 19,06 15,15 25,96 22,94 14,73 28,15 118,59 29,23 33,72

where C is the covariance matrix that represents the possible noise, N is the amount of epochs, v stands for the residua vector of linear function obtained by least squares method for the same covariance matrix C. In practice, we maximize the logarithm of likelihood function:

k

where κ is the spectral index,

The spectral index which characterizes the obtained noise model can range in: -3 < κ< -1 for „fractional Brownian” motion – the nonstationary processes, -1 < κ< 1 for „fractional Gaussian” processes – the stationary or independent of time processes, κ=0 for classic uncorrelated white noise with flat power spectrum, κ=-1 for „flicker” noise, κ=-2 for classical „random walk” or Brownian motion.

„flicker „ran d

om

” noise walk

Spectral index: κ = -0.27 0.25 Noise amplitude: A = 2.35 mm/yr

1 ln[l ( x, C )] = - [ln(det C ) + vˆT C -1vˆ + N ln(2p )] 2 CATS provides the following results: the values of log-likelihood function for defined model or the best spectral index, the intersection of the graph with the vertical axis (estimated on the beginning of the year that first data come from), the graph slope, the amplitudes of sine and co-sine functions corresponding to defined periods, the amplitudes of white and power-law noises as well as spectral index. Beyond above, programme 1/4 ensures the one sigma uncertainties for all mentioned values. The noise amplitudes are given in mm for white noise, in mm/year for „flicker 1/2 noise” and in mm/year for „random walk”.

white noise

”

Spectral index: κ = -0.59 Noise amplitude: A = 3.53 mm/yr0.25


Noise recognized by spectral index for which the log-likelihood function was maximized is stated to be the one occurring in the time series.

amplitude

The classical white noise κ=0 has a flat power spectrum. It appears in nearly all natural phenomena but does not bring any significant information. It is a strictly random process. It can be reduced by increasing the number of data.

amplitude

DATA

The power spectral density of „flicker noise” κ=-1

amplitude

1 is inversely proportional to frequency: f

frequency [Hz]

amplitude

time [s]

Data taken to this research were the topocentric coordinates in ITRF2005 Reference Frame of 10 selected Polish permanent stations belonging to EPN (EUREF Permanent Network) and covered the period of one year (from June 2008 till June 2009). The time series are characterized by trend as well as seasonal, cyclic and random variations. In the analyzed data the systematic component, which interests us the most and the random noise, which is a kind of disturbance that does not bring any significant information were isolated using CATS (Create and Analyze Time Series) software. The coordinates were calculated not in the standard daily, but 1-hour solutions to determine changes in the sub-diurnal frequencies. This kind of data processing allowed authors to estimate the amplitudes of frequencies equal to 24- and 12-hours, the slopes and intercepts of functions, spectral indexes of power-law noise and amplitudes of both white and coloured noise. No a-priori assumption on the spectral index was made.

GWWL station

Spectral index: κ = -0.97 Noise amplitude: A = 8.63 mm/yr0.25

time [s]

frequency [Hz]

amplitude

The power spectral density of “random walk” κ=-2 is inversely 1 proportional to the square of frequency: f2

amplitude

A combination of flicker plus white noise is considered to describe the noise characteristic in Global Positioning System (GPS) the best.


time [s]

frequency [Hz]

The geodetic monument instability has two major compounds: a seasonal one and a random walk one. This local monument motion is considered to follow the random walk process which corresponds to a strictly nonstationary process and affects the estimations of geodetic uncertainties.

xi = x0 + rti + e x (ti )

Every single station coordinate component can be estimated by equation: where εx(ti) is the error term in time ti. It is a combination of uncorrelated unit-variance random variables α and a coloured noise β:

e x (ti ) = aa (ti ) + bk ¹ 0 b (ti )

where a and b are the amplitudes of white (κ=0) and coloured (κ≠0) noise respectively. Due to above, the covariance matrix Cx can be written as:

Average κ=-0.80

MLE RESULTS

2 bold f -1 Pold = 2p 2

2 k

where I is the identity matrix and Jκ is the covariance matrix for the coloured noise of spectral index κ.

n

L = Õ f ( xi , l )

p

2 bnew f -1 f s × 24 × 60 × 60 × 365 .25

1 L=Õ e i =1 s i 2p

On the assumption that the probability distribution is normal, the likelihood function will take the form:

n

- ( xi - l ) 2 2s i2

USDL station


14˚

16˚

18˚

22˚

20˚

14˚

24˚

16˚

18˚

22˚

20˚

14˚

24˚

16˚

18˚

22˚

20˚

24˚

i =1

( xi - l ) - å ln( 2p s i ) l = -å 2 2s i i =1 i =1

After taking into account the normal probability distribution, the logarithmic function will take the form:

118,59

n

i =1

n

281,33

Tables present the results obtained in CATS software. Programme provided the values of log-likelihood function maxima, the intercept and slope for each of functions. Moreover, the amplitudes for white noise were also estimated but only for USDL station. „x” in table indicates that the amplitude was lower than 1 μm or programme treated κ parameter between integer values as the sum of coloured and white noise.

L = ln L = ln Õ f ( xi , l ) = å ln f ( xi , l )

Usually the estimated parameters are determined by calculation of logarithm of likelihood function:

65,44

2p

i =1

n

2

n

5.00 REDZ

54˚

8.63 GWWL

54˚ 54˚

5.53

10.49 GWWL

7.34

54˚ 54˚

7.28

3.88

52˚ 52˚

BOGI BOGO

3.51 WROC

WROC

9.11

Spectral index

In maximum likelihood estimation approach, parameters that maximize the likelihood function are stated to be the real. The second derivatives of logarithm of likelihood function express the dispersion of likelihood function and can be used in determination of approximate uncertainties.

-1

65.44 USDL

Spectral index

20.17

+1 0

Noise amplitude

-1

16˚

18˚

20˚

22˚

24˚

19.06

52˚ BOGI BOGO

The spectral indexes for topocentric coordinates indicate the existence of power-law noise. κ parameters for N component with the average value of -0.80 correspond to well-known combination of flicker plus white noise that characterize the GPS system. Similar situation can be noticed for E component. Spectral index κ has the average of -0.62 for East. For U component few of stations have κ close to „random walk” and therefore correspond to monument instability in vertical coordinates.

29.23

50˚ 50˚ 281.33 USDL

Spectral index

20.17 Noise amplitude

-1

20.17

16˚

18˚

20˚

22˚

24˚

48˚ 14˚

16˚

18˚

Worth indicating is USDL station which is highly affected by „random walk” with the amplitudes greater than on any other station.

ACKNOWLEDGMENTS

Noise amplitude

48˚ 14˚

20˚


118.59 USDL

ZYWI

+1 0

50˚

22.94 KRAW

33.72

U component

48˚ 14˚

The slopes of functions are quite obvious for horizontal components and correspond to well-known trend from plate motion (data are in ITRF2005 reference frame) whereas for vertical ones can indicate some geophysical or artificial effects related to the monument type or antennae stabilization. The fact if slope for vertical component can be interpreted as real trend needs further analysis.

WROC

3.31 KRAW 7.72 ZYWI

E component

ZYWI

+1 0

50˚ 50˚

6.50 KRAW

Maximum likelihood estimation implemented in CATS software provides estimation of the reliability of GPS stations using two parameters: spectral index of noise κ and its amplitude.

15.15 BOR1

6.15

N component

23.06

GWWL

BOGI BOGO

50˚

54 ˚ 14.73 LAMA

25.96 2.35 BOR1

52˚ 52˚

n

DISCUSSION

28.15 REDZ

2.65 LAMA

3.53 BOR1

52˚

xi åi =1 s 2 i = n 1 åi =1 s 2 i

11.40 REDZ

3.65 LAMA

To determine the function's maxima the derivative of logarithmic function of likelihood for λ parameter should be equal to zero: n n n xi ( xi - l ) 2 dl 1 =0=å = ® = ® lNW 0 l å å 2 2 2 dl si i =1 s i i =1 s i i =1

Pnew =


where: Pnew, Pold are the new and old power spectra, respectively, bnew, bold are the amplitudes of new and old coloured noise, respectively, fs is the sampling interval. ( f s × 24 × 60 × 60 × 365.25)1 / 4 For the diurnal sampling interval [in Hz] the new amplitude bnew will be equal to: bnew = × bold = 1.7440 × bold

MAXIMUM LIKELIHOOD ESTIMATION The maximum likelihood estimation (MLE) maximizes the function of likelihood, which is determined as a product of a-posteriori probability for n samples available:

Average κ=-0.91

As the result of this research, the mentioned parameters from CATS software were obtained for topocentric coordinates with 1-hour sampling interval. For nine of ten stations the „fractional Gaussian” process with κ between -1 and 1 was obtained. USDL station is characterized by fractional Brownian motion and therefore is stated as nonstationary process. Standard CATS input is 24-hour sampling interval data. To process 1-hour data it has to be previously re-scaled. The amplitude of coloured noise is estimated as follows:

Cx = a I + b Jk 2

Average κ=-0.62

22˚

24˚

This research is financed by the Ministry of Science and Higher Education, grant No. 2011/01/B/ST10/05384.
