Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.gao.spb.ru/english/personal/malkin/publ/zm13m.pdf Äàòà èçìåíåíèÿ: Fri Nov 8 18:01:00 2013 Äàòà èíäåêñèðîâàíèÿ: Thu Feb 27 23:21:32 2014 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: ìåòåîðèòèêà

Using Modified Allan Variance for Time Series Analysis
Z. Malkin

39

Abstract

Allan Variance (AVAR) was introduced more than 40 years ago as an estimator of the stability of frequency standards. Now it is also used for investigations of time series in astronomy and geodesy. However, there are several issues with this method that need special consideration. First, unlike frequency measurements, astronomical and geodetic time series usually consist of data points with unequal uncertainties. Thus one needs to apply data weighting during statistical analysis. Second, some sets of scalar time series naturally form multidimensional vector series. For example, Cartesian station coordinates form the 3D station position vector. The original AVAR definition does not allow one to process unevenly weighted and/or multidimensional data. To overcome these deficiencies, AVAR modifications were proposed in Malkin (2008. On the accuracy assessment of celestial reference frame realizations. J Geodesy 82: 325­329). In this paper, we give some examples of processing geodetic and astrometric time series using the classical and the modified AVAR approaches, and compare the results.
Keywords

Allan variance Time series analysis

1

Introduction

The scatter of a geodetic time series provides a good measure of the series quality (and its more explicit statistical characteristics). One of the most effective approaches to analyze the noise component of a measured signal is the Allan Variance (AVAR) originally developed for the evaluation of the stability of time and frequency standards (Allan 1966). In particular, its advantage is the weak dependence of the noise parameter estimate on low-frequency components of the signal under study. AVAR has already had a rather long history in geodesy and astrometry. It was used to study station coordinates time series (Malkin and Voinov 2001; Le Bail and Feissel-

Z. Malkin (*) Pulkovo Observatory, St. Petersburg, 196140, Russia e-mail: malkin@gao.spb.ru

Vernier 2003; Le Bail 2006; Feissel-Vernier et al. 2007), radio source position stability (Feissel et al. 2000; Gontier et al. 2001; Feissel-Vernier 2003), the Earth orientation parameters (EOP) series (Gambis 2002), geocenter motion (Feissel-Vernier et al. 2006). AVAR can be also used to analyze the spectral characteristics of the signal (Feissel et al. 2000; Feissel-Vernier 2003; Feissel-Vernier et al. 2006, 2007) and Hurst parameter (Bregni and Primerano 2005). A detailed review and history of using AVAR in astrometry and geodesy can be found in Malkin (2011). However, application of AVAR to analysis of the time series of geodetic and astrometric measurements yields sometimes unsatisfactory results because the analyzed series consists of the data points with unequal uncertainties. This necessitates proper weighting not provided by the original AVAR definition. Besides, we often deal with multidimensional values, e.g., terrestrial coordinates and/or velocities (3D or 6D), celestial coordinates and/or proper motions (2D or 4D). To provide adequate analysis of such kind of data,
271

Z. Altamimi and X. Collilieux (eds.), Reference Frames for Applications in Geosciences, International Association of Geodesy Symposia 138, DOI 10.1007/978-3-642-32998-2_39, # Springer-Verlag Berlin Heidelberg 2013


272

Z. Malkin

AVAR modifications WAVAR, MAVAR, and WMAVAR have been proposed by Malkin (2008a). They are explained in Sect. 2. In Sect. 3, several practical examples are given to show how proposed AVAR modifications can be applied to investigation of real geodetic and astrometric data. These examples also allow us to better understand some practical features of various types of AVAR estimates. We do not consider here the applications of AVAR to spectral or persistency analysis. Such applications require computation of modified AVAR on different averaging intervals, which should be straightforward.

Strictly speaking, the weights pi should be computed in accordance with the error propagation law. ! k i2 j 2 j 2 !' þ1 X & h j j pi ¼ yi þ yi×1 = d i si × si×1 :
j ¼1

(39.4) However, the expression Eq. 39.4 has a singular point in the case of di ¼ 0, i.e. of equal adjacent measurements. Hence, this case require special treatment, e.g. assigning unit weight for close (near-)zero value of di. To avoid the singularity, the following simplified formula for pi is recommended in Malkin (2008a) for practical use: !þ1 k X j 2 j 2 ! pi ¼ si × si×1 : (39.5)
j¼1

2

AVAR and Its Modifications

In this section, we give definitions of AVAR modifications proposed in Malkin (2008a) to analyze real time series. The classical AVAR is defined as s2 ¼
n þ1 X 1 Ïy þ y 2Ïn þ 1÷ i¼1 i

i ×1

÷2 :

(39.1)

where yi are measurements, i ¼ 1,...,n. This definition is not satisfactory for the analysis of real measurements, which generally have different uncertainties and thus should be unevenly weighted during analysis to obtain realistic and statistically meaningful estimates for the investigated parameters. For such kind of data, the weighted AVAR, WAVAR, is introduced in Malkin (2008a). It can be defined as s2 ¼ 1 p¼
n þ1 1X p Ïy þ yi 2p i ¼ 1 i i

No significant difference between results computed with weights given by Eqs. 39.4 and 39.5 was found during real data processing. In the case of pi ¼ 1, we have unweighted multidimensional AVAR (MAVAR). It is easy to see that WMAVAR is a universal definition which includes all others as special cases. It can be also noted that for homogeneous time series with close si, MAVAR (WMAVAR) values should be equal to AVAR (MAVAR) ones computed for one of k vector dimensions and multiplied by k. In noise component analysis, one often uses Allan Deviation ADEV computed as the square root of AVAR. Correspondingly, ADEV modifications, WADEV, MADEV, and WMADEV can be used in most practical applications, including this study.

×1

÷2 ; (39.2)
þ1 2 i×1 ÷

n þ1 X i¼1

3

Practical Examples

pi ;

pi ¼ Ï

s2 i

×s

: In this section, we give several examples of using modified ADEV. Four examples will be considered: comparison of station displacement time series, investigation of source position stability, comparison of celestial pole offset (CPO) time series, and quality assessment of celestial reference frame (CRF) realization (radio source position catalogues).

where si are the measurement uncertainties. Another modification of AVAR proposed by Malkin (2008a) is intended for processing of multidimensional data. The latter can be considered as a k-dimensional vector þ à yi ¼ y1 ; y2 ; ... ; yk ; i ¼ 1, ..., n, with standard errors i i i si ¼ s1 ; s2; ... ; sk : The corresponding weighted multidii i i mensional AVAR, WMAVAR, is given by s2 ¼ 2
n þ1 1X p d2 ; 2p i ¼ 1 i i i×1

3.1

Station Position Series

d i ¼ yi þ y

;

p¼

n þ1 X i ¼1

pi :

(39.3)

ADEV is an effective tool for investigation of noise characteristics of station position time series (Malkin and Voinov 2001; Le Bail and Feissel-Vernier 2003; Le Bail 2006; Feissel-Vernier et al. 2007). However, the classical ADEV applied to a series with heterogeneous


39

Using Modified Allan Variance for Time Series Analysis
20 15

273 Table 39.1 Scatter characteristics for EPN station position time series shown in Fig. 39.2. Unit: mm Station JOZE ISTA HFLK ADEV X 2.9 10.2 5.2 Y 1.1 17.1 1.9 Z 2.7 24.0 5.6 WADEV X 1.6 2.3 20.2 Y 1.0 1.5 2.9 Z 2.0 1.9 21.2 3D WADEV 2.8 3.3 28.9

ADEV = 4.18 WADEV = 2.66

Station height, mm

10 5 0 -5 -10

1

3

5

7

9 11 13 Week number

15

17

19

Fig. 39.1 Difference between ADEV and WADEV in case of station height time series with implied outliers

X
JOZE, X 0.4 1 0.9 0.8 1000 1200 1400 1600 1000 0.3 0.2 0.1

Y
JOZE, Y 0.5

Z
JOZE, Z

JOZE

0.4

scatter statistics is given in Table 39.1. Station JOZE shows stable behavior without jumps and outliers; for this station all the estimates give close results. Station ISTA has two outliers of several decimeters with large uncertainties; for this station using unweighted estimate gives unsatisfactory result, and using WAVAR allows us to practically eliminate outliers. The last case of station HFLK shows in contrast to the previous example unsatisfactory result obtained with WAVAR. The reason is that the HFLK position uncertainties were relatively large in the period when the station showed position stability and became much smaller to the end of the series where the position jump occurred. As a result, the position estimates around the jump epoch were entered to the statistics with large weight.

Station

1200

1400

1600

1000

1200

1400

1600

GPS week ISTA, X

GPS week ISTA, Y

GPS week ISTA, Z

0.6 0.5

ISTA

0.4 0.3 0.2 1000 1200 1400 1600 GPS week HFLK, X
0.3 0.2 0.1

0.8 0.7 0.6 0.5 0.4 0.3 0.2 1000 1200 1400 1600 GPS week HFLK, Y
0.8

0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 1000 1200 1400 1600

3.2
GPS week HFLK, Z

Source Position Series

0.4
0.7

HFLK

0.3 0.2

0 -0.1 1000 1200 1400 1600 1000 1200 1400 1600 1000 1200 1400 1600

GPS week

GPS week

GPS week

Fig. 39.2 EPN station position time series. X, Y, and Z are the station Cartesian geocentric coordinates. Unit: m

uncertainty can give a biased estimate of the actual signal characteristics. As the first example, let us consider a series of weekly station height estimates depicted in Fig. 39.1. One can see that one to three points may be considered as excludable because of large bias. Their inclusion in detailed analysis (e.g. scatter statistics) depends on the method used to detect outliers. Using WADEV mitigates the impact of implied outliers without application of special statistics. Now, let us consider three station position time series provided by the European Permanent GPS Network (EPN) Central Bureau.1 These series are shown in Fig. 39.2, and

1

http://epncb.oma.be/

Investigation of the noise characteristics of source position time series are usually aimed at ranking of sources by time series statistics, and compiling list of sources that are not stable enough to be solved in VLBI global solution, and require special handling. In particular, computation of some quantitative source stability indices are desirable to make a selection of the International Celestial Reference Frame (ICRF) defining sources as objective as possible (Ma et al. 1998; Feissel-Vernier 2003; Gordon et al. 2008; Malkin 2008b). In Malkin (2009) we investigated source position time series submitted by VLBI analysis centers in the framework of preparation of the Second ICRF realization ICRF2 (Ma et al. 2009). These series were computed using different software and/or analysis strategies, which makes their comparison especially interesting and instructive. Each analyzed series consists of source positions obtained from analysis of 24-h observing sessions (one source position estimate for each session). Out of several scatter indices generally used for analysis of source position time series we used the following two: weighted root-mean-square (WRMS) residuals of the session source position with respect to weighted mean position (the weights are inversely proportional to the square of the reported position uncertainty), and WADEV. Our study


274

Z. Malkin

(Malkin 2009) confirmed that both statistics have their own advantages and disadvantages. WRMS estimate is affected by trend-like and low-frequency signal components. In contrast, WADEV does not depend on the slow position variations. However, it may give inadequate estimate of the scatter index in the case when the time series contains jumps but is stable between jumps. In order to get a more general measure of source position stability, a composite index of WRMS and WADEV can be used. It can be noted that some examples given above may seem to be artificial. It is common practice to identify and reject outliers before further statistical analysis. However, as already mentioned above, it can be not a trivial in all the time series considered here, and others we meet in our work. So, our intention is to show that weighted ADEV (WADEV) provides more robust estimate in case of outliers, provided the outliers have large uncertainties, which is quite common, see e.g. Figs. 39.1 and 39.3 (next section). In such a case, the result obtained with WADEV is more independent of the quality of outlier detection procedure used.

3.3

CPO Series Comparison

Investigation of the noise characteristics of CPO time series is important for different tasks, such as quality assessment of EOP series and weighting EOP series during combination (Gambis 2002). Here we consider ten series of the CPO coordinates dX and dY computed at International VLBI Service for Geodesy and Astrometry (IVS) analysis centers,2 including the IVS combined series. They are depicted in Fig. 39.3 in a deliberately small scale to show most of the outliers. A visual inspection of Fig. 39.3 reveals, for example, that OPA series seems to be the noisiest, and IAA series seems to be practically noiseless. However, a deeper look at the data shows that the first impression is wrong. First of all, the series contain different number of sessions. One also notices that most of the outliers have enormously large uncertainties. Both factors should be taken into account. In particular, it requires to compute weighted scatter index and compare subsets of original series consisted of the sessions common for all compared CPO series. In Tables 39.2 and 39.3, the results of computation of several statistics are presented. The former table shows the impact of data weighting only. The latter table provides scatter indices computed using the same data set for all the series, i.e. free of selection effect. The statistics used are WRMS, MADEV (unweighted 2D ADEV estimate) and WMADEV (weighted 2D ADEV estimate). WRMS values given in the table are computed as average of those for dX and dY. In fact, WRMS1 is just the WRMS value of the reported CPO estimates. Every statistics is computed in two variants: for

Fig. 39.3 VLBI CPO time series. Units: mas


39

Using Modified Allan Variance for Time Series Analysis

275

Table 39.2 Statistics of the VLBI CPO time series depicted in Fig. 39.3. The statistics indexed with 1 are computed from raw data, the statistics indexed with 2 are computed after removing CPO model and linear trend. Units: mas Series AUS BKG CGS GSF IAA IVS OPA PUL SPU USN No. of sessions 862 1212 1083 1276 1101 1158 1437 1167 854 1208 WRMS1 251 191 224 197 211 202 189 213 253 195 WRMS2 113 101 129 86 97 80 87 87 115 86 MADEV1 206 245 301 279 173 236 680 215 227 200 MADEV2 206 245 300 279 173 236 680 215 227 200 WMADEV1 177 163 200 137 146 125 149 137 180 132 WMADEV2 177 163 200 137 146 126 148 137 180 132

Table 39.3 The same as Table 39.2, but computed for common sessions Series AUS BKG CGS GSF IAA IVS OPA PUL SPU USN No of sessions 740 740 740 740 740 740 740 740 740 740 WRMS1 251 197 226 199 214 209 189 217 254 192 WRMS2 113 101 128 82 98 79 78 90 114 85 MADEV1 195 180 222 145 165 139 141 152 223 146 MADEV2 195 180 221 145 165 139 141 152 223 146 WMADEV1 176 159 203 131 151 125 128 139 179 132 WMADEV2 176 158 203 130 151 125 127 139 179 132

Table 39.4 Statistics of the two CPO time series computed with two source catalogues ICRF1 and ICRF2. Unit: mas CRF ICRF1 ICRF2 WADEV(dX) 102 92 WADEV(dY) 107 90 WMADEV 149 129

One more result not included in these tables is that 2D WMADEV estimates are very close to the average of WADEV estimates computed for dX and dY multiplied by 2. Thus the multidimensional ADEV provides more compact expression for noise characteristics.

original CPO estimates as reported by authors (indexed with 1), and for values corrected for CPO model ZM2 (Malkin 2007) and linear trend (indexed with 2). A detailed discussion of the results of these computations is beyond the scope of this paper. Let us just mention several conclusions. ­ The comparison of "1" and "2" variants confirms that WRMS depends heavily on the model of systematic errors, whereas WADEV and WMADEV do not. ­ Using weighted ADEV estimates allows us to severely mitigate the influence of outliers. It is especially visible for the OPA series which includes many CPO estimates having large bias and uncertainty (coming from sessions with poor geometry), evidently for completeness. However it affects less the CGS statistics because the latter includes CPO estimates with large bias but small uncertainty (see Fig. 39.3). ­ Weighted ADEV estimate provides robust statistics for ranging CPO series which does not practically depend on low-frequency components of the time series. In particular, it can be mentioned that combined IVS series shows the least noise level.

3.4

Comparison of CRF Realization

Noise level comparison of CPO time series computed with different radio source catalogues is one of very few absolute methods (maybe the only method proposed so far) for quality assessment of the CRF realizations. In Malkin (2008b), Sokolova and Malkin (2007), this method was used for the first time to compare several catalogues. Here we use it for comparison of the first and the second ICRF realizations. For this purpose, we computed two CPO series using ICRF1 (ICRF-Ext.2, Fey et al. 2004) and ICRF2 (Ma et al. 2009) source positions. Then we computed their scatter level using modified ADEV estimates WADEV and 2D WMADEV. From the results presented in Table 39.4, one can see that the scatter of the CPO series computed with ICRF2 is substantially smaller than for series computed with ICRF1. We consider this result as clear evidence of better accuracy of the ICRF2 source positions as compared with the ICRF1.


276

Z. Malkin Bruyninx C, Roosbeek F (2006) The EUREF permanent network: recent achievements. EUREF Publication No. 16, Mitteilungen des Bunde¨ ¨ samtes fur Kartographie und Geodasie, Band, vol 40, pp 105­112 Feissel M, Gontier A-M, Eubanks TM (2000) Spatial variability of compact extragalactic radiosources. Astron Astr 359:1201­1204 Feissel-Vernier M (2003) Selecting stable extragalactic compact radio sources from the permanent astrogeodetic VLBI program. Astron Astr 403:105­110. doi:10.1051/0004-6361:20030348 Feissel-Vernier M, Le Bail K, Berio P et al (2006) Geocentre motion measured with DORIS and SLR, and predicted by geophysical models. J Geodesy 80:637­648. doi:10.1007/s00190-006-0079-z Feissel-Vernier M, de Viron O, Le Bail K (2007) Stability of VLBI, SLR, DORIS, and GPS positioning. Earth Planets Space 59:475­497 Fey AL, Ma C, Arias EF et al (2004) The second extension of the international celestial reference frame: ICRF-Ext.2. Astron J 127:3587­3608. doi:10.1086/420998 Gambis D (2002) Allan variance in earth rotation time series analysis. Adv Space Res 30:207­212. doi:10.1016/S0273-1177(02)00286-7 Gontier A-M, Le Bail K, Feissel M, Eubanks T-M (2001) Stability of the extragalactic VLBI reference frame. Astron Astr 375:661­669. doi:10.1051/0004-6361:20010707 Gordon D, Ma C, Gipson J, Petrov L, MacMillan D (2008) On selection of "defining" sources for ICRF2. In: Finkelstein A, Behrend D (eds) Proceedings of fifth IVS general meeting, pp 261­264 Le Bail K (2006) Estimating the noise in space-geodetic positioning: the case of DORIS. J Geodesy 80:541­565. doi:10.1007/s00190006-0088-y Le Bail K, Feissel-Vernier M (2003) Time series statistics of the DORIS and GPS colocated observations. Geophysical Reseaech Abstracts, EGS-AGU-EUG Joint Assembly, vol 5, Nice, 6­11 April 2003, p 04078 Ma C, Arias EF, Eubanks TM et al (1998) The international celestial reference frame as realized by very long baseline interferometry. Astron J 116:516­546 Ma C, Arias EF, Bianko G et al. (2009) The second realization of the international celestial reference frame by very long baseline interferometry. In: Fey A, Gordon D, Jacobs CS (eds) Presented on behalf of the IERS/IVS working group, IERS technical note no. ¨ 35, Frankfurt am Main, Verlag des Bundesamts fur Kartographie ¨ und Geodasie Malkin ZM (2007) Empiric models of the Earth's free core nutation. Solar Syst Res 41:492­497. doi:10.1134/S0038094607060044 Malkin Z (2008a) On the accuracy assessment of celestial reference frame realizations. J Geodesy 82:325­329. doi:10.1007/s00190007-0181-x Malkin Z (2008b) On construction of ICRF-2. In: Finkelstein A, Behrend D (eds) Proceedings of fifth IVS general meeting, pp 256­260 Malkin Z (2009) Some results of analysis of source position time series. IVS Memorandum 2009-001v01. http://ivscc.gsfc.nasa. gov/publications/memos/ Malkin ZM (2011) Study of astronomical and geodetic series using the allan variance. Kinemat Phys Celest Bodies 27:42­49. doi:10.3103/ S0884591311010053 Malkin ZM, Voinov AV (2001) Preliminary results of processing EUREF network observations using a non-fiducial strategy. Phys Chem Earth (A) 26:579­583. doi:10.1016/S1464-1895(01)00104-1 Schlueter W, Behrend D (2007) The international VLBI service for geodesy and astrometry (IVS): current capabilities and future prospects. J Geodesy 81:379­387. doi:10.1007/s00190-006-0131-z Sokolova J, Malkin Z (2007) On comparison and combination of catalogues of radio source positions. Astron Astr 474:665­670. doi:10.1051/0004-6361:20077450

4

Summary

The modified Allan variation (AVAR) and associated Allan deviation (ADEV) estimators proposed in Malkin (2008a) for processing unevenly weighted and/or multidimensional data is an effective and convenient tool for geodetic and astronomical time series scatter analysis. An important advantage of ADEV is its low sensitivity to low-frequency signal variations as compared to WRMS, which heavily depend on the model used to separate the stochastic and systematic signal components. However, both the original and the modified ADEV may inadequate estimate the noise level when jumps are present the time series. The AVAR (ADEV) modifications, WAVAR (WADEV) and WMAVAR (WMADEV) developed to process unevenly weighted one- and multidimensional measurements, allow us to get noise characteristics that are less sensitive to outliers than the original ADEV estimate, provided the outliers have exaggerated uncertainties, which is usually the case. It is clear that the WMADEV (WMAVAR) definition is the most general as it includes the original AVAR definition and its modifications, WADEV and MADEV, as special cases. For thoroughness' sake, it should be mentioned that the original AVAR is computed on evenly spaced time series. Unfortunately, it is not the case for many geodetic and astronomical applications. Many time series have regular time span however, e. g. daily or weekly station coordinates, daily troposphere parameters, etc. Other series can be made (near)-regular by averaging measured data over equal intervals. For example, such a method was used by FeisselVernier (2003) for radio source position time series. On the other hand, if AVAR and its modifications are used as a series scatter index, uneven spacing should not influence result of analsys, to our mind. And otherwise, averaging the original series may lead to loss of information.
Acknowledgements This work has made use of data and products provided by the International VLBI Service for Geodesy and Astrometry (IVS, Schlueter and Behrend 2007) and European Permanent GPS Network (EPN, Bruyninx and Roosbeek 2006). The author is grateful to two anonymous reviewers for careful reading of the manuscript and helpful comments and suggestions.

References
Allan DW (1966) Statistics of atomic frequency standards. Proc IEEE 54:221­230 Bregni S, Primerano L (2005) Using the modified allan variance for accurate estimation of the hurst parameter of long-range dependent traffic. In: arXiv: cs/0510006, Milano