Introduction
The EUREF Permanent Network (EPN) is a well-dis-tributed and dense network of more than 100 stations. The EPN GPS data are being routinely analyzed by several European Analysis Centers (AC). In addition to the weekly coordinates of the EPN stations, the tropospheric parameters are also estimated. Each AC computes individual solutions for the coordinates and the troposphere parameters, which are subsequently combined in order to yield the final EUREF product, the weekly solutions of combined coordinate and tro-posphere. The EUREF analysis strategy described in Gendt (1997) uses the batch type of processing and is applied in post-processing. Recently new activities aim at hourly GPS data processing from distributed regional networks and focus on meteorological
applications (Ware et al. 2000; Stoew et al. 2001; Gendt et al. 2001; Dousˇa 2001). Thus it seems neces-sary to refine existing algorithms of troposphere solu-tions combination for near real-time applicasolu-tions.
We present a new method for the sequential combi-nation of regional troposphere solutions. This study is motivated by the need of having an internally consistent Zenith Total Delay (ZTD) product on the basis of sev-eral individual troposphere solutions, obtained from near real-time processing. An example is the COST-716 project (COST 2005), in which each AC processed a GPS network in near real-time and delivered ZTD esti-mates every 1 h and 45 min. Using the method presented in this paper, the combination of ZTD data can be performed sequentially as new batches of data become available. To maintain consistency of the combined solutions, appropriate weighting factors are assigned to
M. Keshin
Sequential combination of troposphere
time series
Received: 7 November 2005 Accepted: 9 March 2006 Published online: 7 April 2006 Springer-Verlag 2006
Abstract A new method for sequential combination of tropo-sphere time series is presented. Un-like the EUREF1time series analysis strategy which works in the post-processed mode, the method con-sidered here is capable of combining troposphere solutions sequentially as a new batch of data becomes avail-able. It provides combined tropo-sphere estimates and their standard deviations. In addition to time series biases, the method determines weights to maintain the consistency of combined solutions using vari-ance component estimation. The
time series biases and the weights are determined sequentially using esti-mates obtained at the previous step of time series combination as the a priori information. The mathemati-cal description of the method is presented. The results of experi-mental combination of troposphere solutions obtained by different EUREF and COST2-716 Analysis Centers confirm the ability of the method to combine troposphere time series grouped in daily batches (post-processing mode) and in hourly batches (near real-time mode).
2
The French acronym for European co-operation in the field of scientific and technical research.
1
The IAG (International Association of Geodesy) Reference Frame Sub-Commission for Europe.
M. Keshin
Department of Earth Observation and Space Systems (DEOS),
Delft University of Technology (TU Delft), Kluyverweg,1, P.O.Box 5058,
each contributing ZTD time series. These factors, which can be taken from a regional analysis, are updated sequentially and used to obtain the combined ZTD estimates. Time series weights determination allows one to take into account possible regional and AC depen-dencies, thus providing homogeneous estimates for ZTD and their standard deviations in near real-time.
The performance of the method was tested on a 1-month EUREF dataset for March 2004, by comparing our ZTD estimates with those obtained by the EUREF post-processing combination scheme. We also processed a COST-716 dataset for the same month in order to demonstrate the capability of the method to combine time series grouped in hourly batches in near real-time.
Sequential combination of troposphere time series
Mathematical description
We consider the time series of GPS-derived ZTD esti-mates obtained in-parallel by different AC based on observations from the same GPS station. These ZTD time series become available in batches over a given time interval (hour, day, etc.). One batch contains one or several ZTD estimates from each AC. Ideally the time series within each batch have the same length and sampling interval, and there are no gaps in the data re-cord. Such an assumption is too restrictive because in practice the time series do not necessarily have the same structure. One can specify a set of equidistant time moments (knots) at which to perform time series com-bination. The time interval between knots can be determined, for instance, as the median of the sampling intervals of the time series to be combined. For some time series it may first be necessary to interpolate ZTD estimates to the new sampling interval.
Consider one batch of data comprising K different troposphere time series y(1),y(2),...,y(K). We use bold to indicate a vector quantity. Each time series includes troposphere estimates for N common times t1,t2,..., tN. Denote a troposphere estimate belonging to the solution y(k)for the time moment tnby yn(k)and the corresponding variance by rðkÞ2
n :Let y (k)
be the N·1 vector containing all troposphere estimates from the kth time series and let y be the KÆ N· 1 vector containing all the troposphere estimates taken for combination. The vector y is straightforwardly interpreted as the vector of observa-tional data with the KÆ N · KÆN covariance matrix Qy consisting of the variances of individual ZTD estimates, namely,
Qy ¼ diag ðrð1Þ1 2; . . . ;rð1ÞN 2; . . . ;rðKÞ1 2; . . . ;rðKÞN 2Þ: ð1Þ The unknowns are the combined ZTD estimates (com-bined troposphere time series) and the differences
between the combined and the troposphere time series taken for combination (biases). The observations (ZTD estimates) and the unknowns (combined solutions, biases) are linearly related. Indeed, for the kth time series we have
ynðkÞ¼ Ynþ bðkÞþ eðkÞn ; n¼ 1; . . . ; N ; ð2Þ where Yn is the combined ZTD estimate for the time moment tn; b(k) is the bias estimate for the kth time series, the same for all N times in Eq. (2); en(k)stands for the misclosure error. If we collect all such equations for all the time series, we get the linear functional model. This model can be represented as a Gauss–Markov model assuming observational data with random and Gaussian noise,
E yf g ¼ AYYþ Abb; D yf g ¼ Qy;
ð3Þ
where E and D stand for the expectation and dispersion operators, respectively; Y is the N ·1 vector of the combined troposphere solution; b is the K· 1 vector of the time series biases; AY and Abare the KÆ N · N and KÆN · K design matrices for the parameters Y and b, respectively. They are expressed as
AY¼ INN INN INN 0 B B @ 1 C C A ð4Þ and Ab ¼ eN 0 eN .. . 0 eN 0 B B B @ 1 C C C A; ð5Þ
where INN is the N· N unit matrix and eN is an N-dimensional column vector of ones. The matrices AY and Ab have full rank. Their structure and size remain the same for any combination step, provided the number of time series and the number of common time moments do not change.
We base the sequential time series combination on the Kalman filtering procedure for estimating time series biases. The random walk model (Gelb1974) is adopted to model the time series biases. We have
bi¼ bi1þ li1pffiffiffiffiffiDt; ð6Þ
biases. Therefore, the time series biases are free to move away on average from their initial estimates.
The time series biases estimates from the previous combination step are used as apriori information at the next combination step. The time update equation for the biases then reads
bi;i1¼ bi1; ð7Þ
where bi, i-1 denotes the vector of predicted time series biases at step i. Propagation of the covariance matrix of
bi- 1, Qbi1; is performed as follows:
Qbi;i1¼ Qbi1þ Q; ð8Þ
where, from Eq. (6), Q=Ql Dt=rb2 DtÆ IKK. Using the propagation equations (7) and (8), we can write the partitioned model to compute the combined ZTD solu-tion Yi and the filtered estimate for the vector of time series biases biat step i as
E yi bi;i1 ¼ AY 0 Yiþ Ab IKK bi; D yi bi;i1 ¼ Qyi 0 0 Qbi;i1 ! ; ð9Þ where IKK is the K · K unit matrix. The measurement update equation for the time series biases derived from this model reads
bi¼ bi;i1þ Kiðyi AYYi Abbi;i1Þ ð10Þ with the gain matrix
Ki¼ Qbi;i1A T bðQyiþ AYQYiA T Yþ AbQbi;i1A T bÞ: ð11Þ
The overall validity of the model (9) can be checked by means of the Global overall model (GOM) test sta-tistics, which is expressed in terms of innovations
mi¼ yi AYYi Abbi;i1 ð12Þ
with covariance matrix Qmi ¼ Qyiþ AYQYiA
T
Yþ AbQbi;i1A
T
b: ð13Þ
The GOM test statistics reads (Tiberius1998)
Tj;l¼X l
i¼j
mTiQ1mi mi: ð14Þ
The indexes j and l define a set of consecutive time series combination steps for which the model is validated. The GOM test statistics has a central v2-distribution with P
i=j l
midegrees of freedom; miis the degree of freedom of a functional model for combination step i. If T exceeds the critical value ka, a model error occurred. If
j=l, the local overall model (LOM) test statistics for a single combination step j is obtained.
Time series weights
In the previous section we implicitly assumed that all the Ktime series are taken with the same (unitary) weight. This implies an equivalency that is not the case in practice. Each AC is free to use any GPS data processing software and choose its favorite estimation strategies. This makes ZTD estimates obtained by different ACs relatively biased. In addition, different ACs potentially process different GPS networks which further contrib-utes to the relative biasedness of ZTD estimates. Some factors like non-optimality of strategies or improper modeling of observables may lead to less accurate or even erroneous troposphere estimates. Therefore, the relative offsets between troposphere estimates are mostly AC-dependent. With this in mind, we assign weights to each K time series used in the previous section. We denote below the weight factor for the kth time series by wk. These weights are unknown parameters to be determined in addition to the combined ZTD time series and the time series biases. The estimated weights are then used in the time series combination method as additional factors for the covariance matrix of obser-vations.
We will determine the time series weights by applying the principles of variance component estimation (Koch 1999). In order to show how to perform this, rewrite the functional model (3) in a slightly modified form. After lumping the vectors Y and b into a K+N· 1 vector x and the matrices Aband AYinto a KÆN· K+N matrix A we get
E yf g ¼ Ax; Enðy AxÞðy AxÞTo¼ Qy: ð15Þ Despite the fact that the weight factors are dimension-less, we now formally treat them as ‘variance compo-nents’ of the covariance matrix Qy. It allows us to write the covariance matrix of observations as a sum of K components
Qy¼X K
k¼1
wkQk; ð16Þ
where Qk is the KÆN · KÆN diagonal matrix with N diagonal entries containing the variances of the indi-vidual ZTD estimates from the kth ZTD time series and having the rest entries equal to zero, so that
The matrix Qk is non-negative definite, and the condi-tion forPk=1K Qkbeing positive definite is met. There-fore, Qk can play the role of a cofactor matrix. The stochastic model in Eq. (15) is now given by
Enðy AxÞðy AxÞTo¼X K
k¼1
wkQk ð18Þ
and can be interpret as a linear system of (KÆ N)2 observation equations. For this system we derive esti-mates for wk based on the least-squares principle (Teunissen1988). Least-squares estimators are known to be unbiased and of minimum variance, provided that observables are normal distributed. In addition, due to the use of the least-squares principle we can obtain the covariance matrix of the weight factors, which is the inverse of the normal matrix in this case.
The least-squares solution of the model (18) then reads (Teunissen1988) XK b¼1 nabwb^ ¼ la; a¼ 1; . . . ; K ð19Þ with nab¼1 2trace Q 1 y Q^eQ1y QaQ1y Q^eQ1y Qb ð20Þ and la¼1 2^e T Q1y QaQ1y ^e ð21Þ
where ^e is the residual vector with the covariance matrix
Q^e; ^w¼ ^ wb b¼1;...;K is the least-squares estimator for
the time series weights. The inverse of the matrix whose entries are given by (20) is the covariance matrix of the weight factors Qw^:
Now we apply Kalman filtering to estimate the weights sequentially. Let us consider the K · 1 vector of weight factors wi-1 and the corresponding covari-ance matrix Qwi1 obtained after the (i)1)th combina-tion step. To propagate these estimates to the next ith step, we apply the random walk model (6) in which the vectors bi-1 and bi should now be replaced by wi-1 and wi, respectively. Adoption of the random walk model means that the weights are free to move away on average from their initial values. Indeed, possible future modifications in data analysis strategies or other improvements in the troposphere solutions at contributing ACs will cause changes in the corre-sponding weight factors. Moreover, because such events cannot be predicted, it is reasonable to assume that the initial uncertainty of weight estimates grows unlimitedly with time. This assumption is in agreement with the fact that the random walk is not stationary and has an infinite variance.
The predicted vector of the weights is expressed as
wi;i1¼ wi1: ð22Þ
Propagation of the covariance matrix Qwi1 is performed similar to that of the covariance matrix Qbi1[see Eq. (8)]
Qwi;i1 ¼ Qwi1þ Q ¼ Qwi1þ rw2Dt IKK; ð23Þ
where rwis the random walk process noise drift for the time series weights.
The measurement update equation for the weights can be expressed in the same form as (19). For the ith time series batch, we get
XK b¼1 ni ½ ab½ wi^ b¼ li½ a; a¼ 1; . . . ; K ð24Þ with ni ½ ab¼1 2trace Q 1 y Q^eQ 1 y QaQ 1 y Q^eQ 1 y Qb þ Q1w i;i1 h i ab; a; b¼ 1; . . . ; K ð25Þ and li ½ a¼1 2^e T Q1y QaQ1y ^eþ XK b¼1 Q1w i;i1 h i ab wi;i1 b; a¼ 1; . . . ; K: ð26Þ
Here [•]b, [•]a bstand for the bth component of a vector and the (a b)th component of a matrix, respectively; Qa and Qb are the cofactor matrices corresponding to the ath and bth time series, respectively. That is, compared to (20) and (21), Eqs. (25) and (26) include additional terms to reflect the impact of the a priori information from the previous combination step.
Summary of the sequential combination method Below we give a short summary of the method to com-bine troposphere time series. The combination is per-formed in the following stages:
1. set-up all relevant apriori information; check input ZTD estimates and exclude those with a standard deviation larger than 10 mm;
is terminated and all relevant apriori information kept unchanged for the next combination step; 3. perform time series weights estimation using a priori
information as well as the residuals along with the corresponding covariance matrix estimated at stage 2. The weights estimation process is iterated until con-vergence;
4. the time series biases, the weights and the corre-sponding covariance matrices are stored and subse-quently used at the next combination step as the a priori information.
Results
EUREF troposphere time series combination
For the first test of the method, we combined weekly ZTD time series for GPS weeks 1,260–1,263 (March, 2004) obtained by different ACs. The following two stations were considered: (1) BORK (Germany) pro-cessed by BKG3, GOP and LPT analysis centers, and (2) ISTA (Turkey) processed by BEK, BKG, OLG and SGO analysis centers. The troposphere solutions for these two EUREF stations were combined separately in daily batches over a time of 4 weeks. We assumed that the time series biases vary slightly around the (unknown) mean and do not change sharply between two consecu-tive combination steps. From this point of view, we found it sufficient to set the random walk drift of the time series bias to 0.32 m per sqrt (day). A further in-crease of the bias drift did not affect the time series bias estimates. In contrary, too low drift values did not allow the biases to change over a whole combination time span. As for the time series weights, additional experi-ments showed a minor impact of the random walk drift on the sequential weights estimates. We set the random walk drift for the weight factors to 10)5per sqrt (day), because low weight drift values delivered lesser values for the GOM test statistics. The time series with a 1-h sampling interval were considered. Therefore, the weekly time series taken for combination contained 168 ZTD estimates. No measurements (ZTDs) were rejected as outliers during the combination. As an outcome, we obtained combined ZTD estimates and their standard deviations.
Figure1 demonstrates typical histograms of ZTD differences between the weekly contributing time series and the corresponding combined solutions. A pre-liminary analysis of the ZTD differences for different weeks and ACs showed that they are mostly non-nor-mally distributed both in terms of skewness and
kurtosis. Further studies are required to establish pos-sible factors affecting the distributions of the differences. Figure2 shows the accuracy of the combined ZTD estimates. It is worth noting here that the standard deviations of the combined estimates are of similar
–15 –10 –5 0 5 10
0 10 20
#
Station name = BORK (EUREF).
–15 –10 –5 0 5 10 0 10 20 # –15 –10 –5 0 5 10 0 10 20 # ZTD differences [mm] AC: BKG mean = 1.4 mm AC: GOP mean = 0.6 mm AC: LPT mean =2.2 mm
Fig. 1 Histograms of the ZTD differences between the troposphere time series obtained by BKG (upper), GOP (middle), and LPT (lower) analysis ceters (ACs) and the combined troposphere solution for station BORK and GPS week 1,260
1260 1261 0.6 0.7 0.8 0.9 1.0 1.1 1.2 GPS week stdev in mm
Station name = BORK (EUREF)
Fig. 2 The standard deviations of the combined ZTD estimates shown in Fig.1for station BORK and GPS week 1,260. The units are in millimeter
magnitude as the standard deviations of the contributing estimates (about 1 mm).
The differences between the contributing troposphere time series and the combined solutions estimated with our method were averaged over a 1-week interval and the resulting means (offsets) were then compared to the differences between the same contributing time series and the official EUREF combined troposphere solu-tions, (e.g., EUREF2005). The results are presented in Table1.
This table reveals a reasonably good agreement between the combined solutions provided by our near real-time combination method and those provided by the post-processed EUREF method. In all cases the offsets agree to within 0.1 mm. Such agreement is to be expected because, as a matter of fact, at each sequential time series combination step the EUREF analysis strategy was applied. As one can see from Table 1, the weekly means exhibit quite significant week-to-week variations what allows us to conclude that low values of the random walk drift are not suitable to represent the real behavior of the time series biases because any time variations are suppressed in this case.
We validated the mathematical model (9) for sequential combination of the weekly time series by means of statistical testing. The null hypothesis H0 that the model is completely and correctly specified was tested. No particular alternative hypothesis was
considered. The validity of the null hypothesis was then checked by means of the Global overall model (GOM) test statistics (14). We did not take into consideration the first batches of data when computed the GOM test statistics because of the absence of predicted time series biases. The critical values kawere computed from a v2 distribution with a 5% level of significance. The results are shown in Table2.
It is easy to see that in several cases for station ISTA the GOM test statistics values exceeded the critical val-ues, indicating the presence of model errors. In order to study this effect in more detail, we computed the LOM test statistics for each time series combination, except for the first. The results presented in Fig.3 show that the LOM test statistics values for station BORK are all below the critical value. In contrary, the LOM test sta-tistics values for station ISTA tend to decrease, and generally drop below the critical value only after a few sequential combination steps. This allows us to conclude that the LOM test statistics may need extra time to converge and more than six components may be required to obtain a reliable GOM test statistics.
COST-716 troposphere time series combination
One of the main goals of the European COST-716 project (COST2005) was to demonstrate the capability
Table 1 The means for the sequential combination method (column nrt) and the weekly EUREF combined solution (column eur) for stations BORK and ISTA
Station Analysis center (AC)
GPS week
1,260 1,261 1,262 1,263
nrt eur nrt eur nrt eur nrt eur
BORK BKG 1.4±1.3 1.5 2.3±2.4 2.4 0.8±2.4 0.9 1.5±1.6 1.5 GOP 0.6±1.3 0.6 1.2±1.6 1.3 0.5±1.8 0.6 1.0±1.8 1.1 LPT )2.2±2.2 )2.1 )3.7±3.3 )3.7 )1.5±2.9 )1.5 )2.6±2.2 )2.6 ISTA BEK 0.8±3.8 0.8 )0.3±3.4 )0.3 )1.0±3.4 )0.9 0.0±3.7 )0.1 BKG 1.4±2.0 1.4 0.9±1.9 1.0 0.5±2.0 0.6 0.8±1.8 0.7 OLG 0.0±4.4 0.0 )0.7±4.1 )0.7 )1.5±5.3 )1.4 0.5±4.4 0.4 SGO )2.3±2.2 )2.3 0.0±2.8 0.0 1.6±1.9 1.7 )1.0±2.5 )1.1
The units are in millimeter
Table 2 Global overall model (GOM) test statistics values computed after combination of the weekly time series for stations BORK and ISTA along with the corresponding v2-critical values k
a
Station GPS week
1,260 1,261 1,262 1,263
GOM ka GOM ka GOM ka GOM ka
BORK 253.6 503.3 477.3 503.3 389.9 503.3 342.0 503.3
of GPS regional networks to provide ZTD estimates in near real-time. Within the context of this project, ten European ACs used their best processing strategies in order to obtain properly validated ZTD estimates.
The main purpose for using the COST-716 time series was to demonstrate that the method can be used for sequential combination of troposphere time series that are available in hourly batches. We used a 1-month COST-716 dataset (March, 2004), and selected three groups of ACs with at least nine jointly processed GPS stations (Table 3).
In general, the various COST-716 ACs use different ZTD sampling intervals. We interpolate troposphere estimates of the contributing time series to equidistant
time moments, the time interval being the median of the ZTD sampling intervals of the time series to be com-bined. If one of the contributing time series happened to have only one ZTD estimate per hour, the combined troposphere solution is computed at only one time mo-ment. As a result, the system of observational equations becomes underdetermined, and the combined tropo-sphere solution is obtained from a priori information.
The troposphere times series for the three groups were combined in hourly batches as if in near real-time. As in the previous case, we assumed that the time series biases vary slightly around the (unknown) mean and cannot change sharply between two consecutive combi-nation steps. We checked different values for the bias random walk drift and found that high drift values (more that 0.1 m per sqrt (hour)) lead to noise-like behavior of the time series bias dynamics and cause significant batch-to-batch variations (1.5–2 cm). Small values for the bias random walk drift, e.g. 10)4, cause batch-to-batch bias variations typically to within 1–2 mm. A further decrease of the bias drift eliminates any batch-to-batch variations. Bearing this in mind, we set the bias random walk drift for our tests to 10)4m per sqrt (hour). As for the time series weights, a minor impact of the random walk drift on the sequential weights estimates was observed. Therefore, we used the same magnitude for the time series weights drift as in the case of the EUREF time series combination, i.e. 10)5, but now the drift is expressed in units per sqrt (hour).
Tables4, 5, and 6show statistics of sequential time series combination in the three groups. It should be noted that for stations CACE, MALL and VALE, marked with asterisk in Table5, the first 300 batches were excluded from the subsequent analysis because of erroneous biases estimates for several contributing time series. Tables4, 5, and 6 further show that the per-centage of successfully processed batches is about 40–50 %. Time series combination failures occurred due to the absence of ZTD estimates at a given hour for one of the GPS stations or ACs. In such cases a priori time series biases and weighting factors were kept fixed for sub-sequent combination steps. No measurements (ZTDs) were rejected as outliers during the combination.
The GOM test statistics values computed for each station included in the combination groups 1–3 exceeded in almost all cases the critical values (v2distribution with a 5% level of significance) indicating the presence of model errors. In order to consider this situation in more detail, we computed the LOM test statistics for all (except for the first) time series combination steps. The results in Tables4, 5, and 6 reveal that model errors occurred for only 10–20% of the batches (see column LOM > ka). An additional analysis showed that batches with too high LOM test statistics values spread over the entire combination time intervals. This indicates that the sequential combination method has a poorer
0 3 6
40 60 80 100
Station name = BORK (EUREF)
day of week LOM teststatistic w1260 w1261 w1262 w1263 0 3 6 50 100 150 200
Station name = ISTA (EUREF)
day of week LOM teststatistic w1260 w1261 w1262 w1263
performance for the COST-716 hourly time series com-bination than for the EUREF time series comcom-bination, apparently due to a lack of observations (ZTD esti-mates) in the former case.
We next demonstrate in Figs.4,5, and6the results of sequential time series combination for the ACs and GPS
stations included in the combination groups 1–3. Shown are the mean time series biases and the weighting factors. The individual bias estimates within each combination group for the same AC agree to within 2–4 mm. The mutual agreement of the time series biases depend on the combination group. The agreement is the best for the first combination group when the AC-dependent offsets are clearly visible. In general, the AC-dependent effects in the time series biases are not very well pronounced, most likely because of different performance of the sequential combination method for different stations. This circumstance requires further study. The weighting factors demonstrate a noticeable agreement not only between different stations within the same combination group, but also for the same AC used in different com-bination groups. In each step, the number of iterations to obtain weight estimates was about 2–3, although sometimes for the initial combination steps it increased to several tens. This reveals a good reliability of the weighting factors obtained by applying variance component estimation.
Additional studies were carried out to gain further insight into the impact of under-determinedness of the system of observational equations. We chose three ACs according to the following criteria: (1) at least three common stations, and (2) at least four ZTD estimates per hour for each AC. The ACs were IEEC, ACRI and ASI, the ZTD samples are 10, 15, and 15 min, Table 3 Analysis centers and
GPS stations involved in the experimental sequential combination of the COST-716 troposphere time series
Group Analysis centers GPS stations involved
Group 1 ASI, BKG, LPT, SGN BRUS, CAGL, GENO, HERS,
KARL, MEDI, PADO, TORI, WTZR, ZIMM Group 2 ASI, IEEC, GFZ, SGN BRUS, CACE, CAGL, ELBA,
MALL, VALE, WTZR, YEBE, ZIMM Group 3 BKG, GFZ, GOP, NKGS BOR1, HERS, HOFN, MAR6,
MATE, MORP, ONSA, POTS, REYK, VIS0, WROC
Table 4 Statistics of sequential time series combination for the combination group 1
Station Batches
Processed Percentage LOM < ka LOM > ka BRUS 421 57 356 65 CAGL 376 51 305 71 GENO 329 44 275 54 HERS 398 53 333 65 KARL 412 55 358 54 MEDI 364 49 298 66 PADO 309 42 256 53 TORI 329 44 288 41 WTZR 394 53 294 100 ZIMM 401 54 343 58
Presented are: the number (processed) and the percentage (per-centage) of successfully processed batches, the number of batches with no significant model errors (LOM < ka), and the number of
batches with significant model errors (LOM > ka)
Table 5 Statistics of sequential time series combination for the combination group 2
Station Batches
Processed Percentage LOM < ka LOM > ka BRUS 384 52 313 71 CACEa 165 37 140 25 CAGL 335 45 291 44 ELBA 362 49 301 61 MALLa 227 51 194 33 VALEa 123 28 97 26 WTZR 247 33 199 48 YEBE 217 29 184 33 ZIMM 357 48 299 58 a
First 300 batches were not taken into consideration
Table 6 Statistics of sequential time series combination for the combination group 3
Station Batches
Processed Percentage LOM < ka LOM > ka BOR1 357 48 333 24 HERS 424 57 398 26 HOFN 418 56 398 20 MAR6 362 49 342 20 MATE 359 48 330 29 MORP 115 15 95 20 ONSA 95 13 85 10 POTS 417 56 389 28 REYK 341 46 302 39 VIS0 357 48 335 25 WROC 307 41 290 17
respectively. The number of epochs at which to perform time series combination is determined as the median of the sampling intervals of the contributing time series. Therefore, in each combination step we had four epochs at which time series combination was performed, and hence 12 ZTD estimates and 7 unknowns, so that the corresponding system of observational equations was overdetermined. Then we added troposphere data from the LPT analysis center and repeated the combination for the same stations, this time using the troposphere solutions from the four ACs. Because LPT provides ZTD estimates once per hour, in each combination step we had only one common epoch, resulting in an underdetermined system of observational equations.
The random walk drift parameters had the same values as in the previous tests. We then compared the GOM test statistics values obtained after the two combination runs. The results are presented in Fig.7.
When the linear system of observations is redundant (upper plot), all the GOM test statistics values turned out to be below the critical values indicating no signifi-cant model errors. In contrary, with the underdeter-mined linear system of observational equations the GOM test statistics fell into the right-handed critical region (lower plot) for three stations, thereby showing the presence of model errors. As for BOR1 station, the corresponding GOM test statistics value also increased, albeit did not exceed the critical value in the second combination run. Therefore, these results provide some –6 –4 –2 0 2 4 6 8 10 bias ZTD [mm]
BRUS CAGL GENO HERS KARL MEDI PADO TORIWTZR ZIMM
ASI BKG LPT SGN 0 2 4 6 8 10 weighting factor
BRUS CAGL GENO HERS KARL MEDI PADO TORIWTZR ZIMM
ASI BKG LPT SGN
Fig. 4 Mean biases (upper) and weighting factors (lower) for the COST-716 ACs and GPS stations included in the combination group 1 –6 –4 –2 0 2 4 6 8 10 bias ZTD [mm]
BRUS CACE CAGL ELBA MALL VALE WTZR YEBE ZIMM
ASI IEEC GFZ SGN 0 2 4 6 8 10 weighting factor
BRUS CACE CAGL ELBA MALL VALE WTZR YEBE ZIMM
ASI IEEC GFZ SGN
evidence of poorer performance of the sequential com-bination method for underdetermined linear systems of observational equations.
Conclusion
We proposed a method for the sequential combination of troposphere time series. Unlike the EUREF analysis strategy which works in post-processing mode, our method combines time series in near real-time as new batches of information become available. The method uses sequential estimation and Kalman filtering to esti-mate time series biases sequentially. In addition to the
time series biases, the weighting factors for the contributing time series are determined using variance component estimation. Like the biases, the time series weights are determined sequentially using a Kalman filtering. In each combination step the time series biases are computed following the post-processed EUREF time series analysis strategy.
Our method was tested with troposphere time series obtained by different EUREF and COST-716 AC. In the former case the time series were combined in daily –6 –4 –2 0 2 4 6 8 10 bias ZTD [mm]
BOR1 HERS HOFN MAR6 MATEMORP ONSA POTS REYK VIS0WROC
BKG GFZ GOP NKGS 0 2 4 6 8 10 weighting factor
BOR1 HERS HOFN MAR6 MATEMORP ONSA POTS REYK VIS0WROC
BKG GFZ GOP NKGS
Fig. 6 Mean biases (upper) and weighting factors (lower) for the COST-716 ACs and GPS stations included in the combination group 3
0 2000 4000 6000
BOR1 CAGL WTZR ZIMM
GOM testst. critical value 0 1000 2000 3000
BOR1 CAGL WTZR ZIMM
GOM testst. critical value
batches, and in the latter case in hourly batches. The comparison of the combined troposphere solutions obtained using the sequential combination method with the weekly EUREF combined solutions showed an agreement to within 0.1 mm. The results of the COST-716 time series combinations demonstrated the ability of the sequential combination method to deliver combined troposphere solutions in near real-time even in the case of underdetermined systems of observational equations. In both tests, the mathematical models were verified by means of statistical testing.
Additional experiments are required to establish all possible factors affecting the sequential combination method, in particular evaluate the impact of different accuracy of contributing ZTD estimates on the per-formance of the method in case of hourly batches combination. Extensions of the sequential combination method will be investigated in the future, for instance, the possibility to adapt it for modeling absolute biases using tropospheric data available from different weather forecasting models such as ECMF (European Center for Meteorology Forecast) and HIRLAM (HIgh Resolution Limited Area Model) (Gustaffson 2002).
List of acronyms of the analysis centers
ACRI ACRI-ST, France
ASI Agenzia Spaziale Italiana, Italy
BEK Bayerische Kommission fu¨r die Internationale Erdmessung, Germany
BKG Bundesamt fu¨r Kartographie und Geoda¨sie, Germany
GFZ GeoForschungsZentrum, Germany GOP Geodetic Observatory Pechny´, Czech
Republic
IEEC Institut d’Estudis Espacials de Catalunya, Spain
LPT Bundesamt fu¨r Landestopographie, Switzerland
NKGS Nordic Geodetic Commission, Onsala Space Observatory, \Sweden
OLG Institute for Space Research, Austria SGN Institut Ge´ographique National, France SGO FOMI Satellite Geodetic Observatory,
Hungary
Acknowledgements This work was supported by the European Commission, through the 5th Framework Program, contract no. EVG1-CT-2002-00080 in support of the TOUGH project.
References
COST (2005) COST-716 Action homepage. Web site: http://www.oso.chalmers.se/ kge/cost716.html/, accessed Oct 2005 Dousˇa J (2001) Towards an operational
near real-time precipitable water vapor estimation. Phys Chem Earth (A) 26(6– 8):189–194
EUREF (2005) Combined EUREF tropo-sphere solutions. Ftp-site: ftp://igs. ifag.de/EUREF/products, accessed Oct 2005
Gelb A (ed) (1974) Applied optimal esti-mation. The M.I.T. Press, Massachu-setts Institute of Technology, Cambridge, Massachusetts
Gendt G (1997) IGS combination of tropo-spheric estimates (Status Report). In: Pro-ceedings IGS analysis centre workshop, Pasadena
Gendt G, Reigber C, Dick G (2001) Near real-time water vapor estimation in a German GPS network-first results from the ground program of the HCF GASP project. Phys Chem Earth (A) 26(6–8):413–416
Gustaffson N (2002) Assimilation of ground-based GPS data in HIRLAM 3D-Var. HIRLAM workshop on vari-ational data assimilation and remote sensing, Helsinki, pp 89–96
Koch K (1999) Parameter estimation and hypothesis testing in linear models. Springer, Berlin Heidelberg New York Stoew B, Jarlemark P, Johansson J, El-gered G (2001) Real-time processing of GPS data delivered by SWEPOS. Phys Chem Earth (A) 26(6–8):493– 496
Teunissen P (1988) Towards a least-squares framework for adjusting and testing of both functional and stochastic models. Internal re-search memo. Geodetic Computing Center, Delft. A preprint of the original report of 1988 is also available in 2004, vol 26, Mathematical Geodesy and Positioning Series. Delft Uni-versity of Technology, the Netherlands Tiberius C (1998) Recursive data processing
for kinematic GPS surveying. PhD Thesis, Delft University of Technology, the Netherlands
