Robust estimation of systematic errors of satellite laser ... - Springer Link

11 downloads 0 Views 100KB Size Report
Abstract. Methods for analyzing laser-ranging residuals to estimate station-dependent systematic errors and to eliminate outliers in satellite laser ranges are ...
Journal of Geodesy (1999) 73: 345 ± 349

Robust estimation of systematic errors of satellite laser range Y. Yang1 , M. K. Cheng2 , C. K. Shum3 , B. D. Tapley2 1

Xian Research Institute of Surveying and Mapping, No. 1 Mid-Yanta Road, 710054 Xian, PRC e-mail: [email protected]; Tel.: 86 29 553 35 03; Fax: 86 29 552 53 10 2 Center for Space Research, The University of Texas at Austin, Austin, TX 78759-5321, USA 3 Department of Civil and Environmental Engineering and Geodetic Science, Ohio University, OH 43210-1275, USA Received: 18 March 1997 = Accepted: 17 March 1999

Abstract. Methods for analyzing laser-ranging residuals to estimate station-dependent systematic errors and to eliminate outliers in satellite laser ranges are discussed. A robust estimator based on an M-estimation principle is introduced. A practical calculation procedure which provides a robust criterion with high breakdown point and produces robust initial residuals for following iterative robust estimation is presented. Comparison of the results from the least-squares method with those of the robust method shows that the results of the station systematic errors from the robust estimator are more reliable. Key words. Robust ®tting á Systematic errors á SLR

1 Introduction The satellite laser-ranging (SLR) technique has demonstrated a unique capability for studying global solid-Earth dynamics, for example in measuring the time-varying properties of the Earth's gravitational ®eld (Schutz et al. 1994, Tapley et al. 1994). SLR data provide the most accurate measurement for precise orbit determination during satellite altimeter missions. Improvements in the analysis of SLR measurements will enhance its contribution to space geodesy research. Pre-analysis of SLR measurements is required to establish the station-dependent noise level for data weighting in orbit determination using SLR data, and generation of the normal points from the full-rate data collected by a laser station in order to reduce the random or white noise content of laser-ranging measurements and also the computation burden in orbit determination and its application (Tapley et al. 1994). Correspondence to: Yuanxi Yang

The analysis of SLR measurements usually follows the procedure in the Herstmonceux Recommendation outlined at the Fifth International Workshop (1984) on laser-ranging Instrumentation. In brief, this procedure involves ®tting the laser-ranging residuals with two preassumed functional models, and is performed in two steps. First, the laser-ranging residuals, which are the di€erences between the values recorded by stations and those computed from a reference trajectory with a certain accuracy, of every pass for each station are ®tted with a straight line de®ned by two parameters, the ranging bias and time bias, by a least-squares (LS) estimation. The estimated ranging and time bias can be used as a measure of the station-dependent systemic errors contributed by the errors in the station coordinate, or by the measurement acquisition and reduction. An extremely large ranging bias or time bias is an indicator of unacceptable error for the corresponding data pass. The remainder residuals obtained at this step are further ®tted with an mthorder polynomial to give the single-pass internal precision of the laser tracking system at the second step. An extremely large polynomial root mean square (RMS) of a pass indicates an unacceptable pass from the tracking system. It is clear that the reliability of this procedure is fundamental to the successful application of SLR data in precise satellite orbit determination and other applications. In order to obtain a reliable result, we have to pay attention to accessing the station-dependent systematic errors and eliminating the outliers. However, the current LS estimation procedure may encounter a number of problems, as follows. 1. If a certain percentage of data are outliers, the parameters in a pre-assumed functional model, the straight line or polynomial estimated from LS procedure, will be a€ected. Therefore, obtained LS estimates of parameters may not reasonably represent the systematic errors in the observations. 2. The LS procedure has the major function of smoothing the outliers in minimization of the sum

346

of the square residuals. Consequently, some large errors may be reduced; in contrast, some small errors may be enlarged, thus the true outlying observations may not be eciently marked or may be rejected. 3. The conventional method for rejecting poor data uses a criterion of 3r as the rejection point. However, in actual fact, the in¯uences of the observations in the neighborhood of the rejection point, for example those data with 3.1r or 2.99r error, are not signi®cantly di€erent. If those data points with the RMS close to the selected criteria are retained, the robustness and eciency of the estimates cannot be guaranteed. Furthermore, the outliers will also distort the estimate of r, which in turn distorts the criterion itself, when using the LS algorithm. In order to obtain reliable results for the station systematic errors with corresponding accuracy as well as a reliable criterion, a robust M (maximum likelihood type) estimation, proposed by Huber (1964), is applied and modi®ed in this contribution. Robust estimation has been widely studied and applied in geodesy. For example, a Danish method (see Kubik 1982) and an IGG-1 scheme (Zhou 1989) for di€erent weight observations were established. An IGG-3 scheme (Yang 1994) was proposed for the estimation for dependent observations. Considering the a priori information on the model parameters, Yang (1991) developed three types of Bayesian robust estimators. Scha€rin (1985) and Yang (1992) studied robust collocation, based on di€erent presuppositions, and obtained di€erent estimators. Koch and Yang (1998) developed a robust Kalman ®lter for rank-de®cient observation models. Robust estimation does not yet have a place in space geodesy research, but it will play an important role in the future, because the measurement process for satellites is time and space varying, and probably has a much more complex error distribution or outlier contamination than that of other static measurements. In the present paper, a brief review of LS estimation in the analysis of satellite laser ranging is given in Sect. 2. The application of robust estimation for the analysis of satellite laser-ranging residuals and assessment of the station-dependent systematic errors is discussed in Sect. 3. The inherent formulas for a practical calculation are presented in Sect. 4, and the results and comparison are given in Sect. 5.

at epoch ti , s_i ˆ osi =oti ; and ei is a random error. The error equation is u…ti † ˆ a^ ‡ b^s_i ÿ Ds…ti †

…2†

where a^ and b^ are the estimates of a and b respectively; u…ti † is the residual of Ds…ti † due to high frequency and random measurement noise. As a second step in the analysis of the SLR observations, it is desired to use an mth order polynomial to remove the remaining errors, which is expressed as follows: u…ti † ˆ a0 ‡ a1 ti ‡ a2 ti2 ‡    ‡ amÿ1 timÿ1 ‡ e0i

…3†

where ai …i ˆ 0; 1; . . . ; m ÿ 1† are the coecients to be estimated and e0i is the error of u…ti †. The error equation is Du…ti † ˆ a^0 ‡ a^1 ti ‡ a^2 ti2 ‡    ‡ a^mÿ1 timÿ1 ÿ u…ti †

…4†

where Du…ti † is a new residual of the original residual u…ti †. In general, the error equation can be expressed as V ˆ X a^ ÿ S

…5†

where X is a design matrix, a^ is the estimated vector of unknown parameters, S is the so-called observation vector and V is the residual vector of S. The elements vi , xi and si , of the vectors of V , X and S, and a^ for the case of Eq. (2), are as follows: vi ˆ u…ti †;

si ˆ Ds…ti †;

xi ˆ ‰1 s_i ŠT ;

^T a^ ˆ ‰^ a bŠ

and for the case of Eq. (4), they are vi ˆ Du…ti †; xi ˆ ‰1

ti ti2

si ˆ u…ti †; . . .ŠT ;

a^ ˆ ‰^ a0 a^1 . . .ŠT

The LS target function is n X iˆ1

pi v2i ˆ min

…6†

where pi is an observation weight (pi ˆ r2 =r2i ; r2 is a variance of unit weight and r2i is the variance of si ). The LS estimator of a reads. a^ ˆ …X T PX †ÿ1 X T PS

…7†

and the estimator of the variance of unit weight is r^2 ˆ

V T PV nÿm

…8†

2 Functional models and their LS estimators

The posterior covariance matrices are as follows:

The functional models used initially to ®t the ranging bias and time bias in the analysis of SLR observation can be described as follows:

Ra^ ˆ …X T PX †ÿ1 r^2

…9†

RS^ ˆ …X Ra^X T †ÿ1

…10†

Ds…ti † ˆ a ‡ bs_i ‡ ei

…1†

where a and b denote the ranging bias and time bias respectively; Ds…ti † is the residual of SLR observation, si ,

where P is an observation weight matrix. Usually, the solution of Eqs. (1) and (3) is based on LS method with the assumption that the model error ei

347

in Eq. (1) and e0i in Eq. (3) follow normal distributions, and the laser-ranging residuals are required to be reliable or of a certain accuracy in the reference trajectory. Otherwise, the estimated values for systematic errors will be systematically distorted due to the large orbit errors. At present, the accuracy of the best reference trajectory has attained an RMS of 5 cm for the Starlette orbit in a 5-day interval ®t, therefore the second assumption is generally satis®ed. However, the ®rst assumption is rarely satis®ed in practice. To overcome the in¯uence of this, we will present a robust method for ®tting the systematic errors of SLR data, in a general case, without the requirement of normal distribution for SLR residuals. 3 Robust ®tting for the systematic errors An alternative target function to Eq. (6) is introduced. In order to resist the in¯uence on estimated parameters from those data with large error, we replace the squares of the residual v2i with a slightly stable function q…vi †, such that

pi ˆ

8 > > pi > > < k0

pi v0 > > j ij > > :0



2 k1 ÿ j j v0i

k1 ÿk0

0 v v ˆ i  k0 i

ri

k0 < v0i  k1 0 v > k1 i

…16†

where k0 and k1 are constants depending on the actual distribution, and usually chosen as k0 ˆ 1:0 to 1.5 and k1 ˆ 3:0 to 6.0 respectively. The ®rst segment in the weight function given by Eq. (16) is the same as that of the LS case. In practice, most of the SLR residuals will satisfy the condition v0i  k0 , and then pi ˆ pi , which is just the result of LS estimation. The second part of the weight function is a downweighing part, which smoothes the weight function from the LS case to a rejection point. In order to obtain the robust variances, r2i , of the observations, the robust estimate of the variance of unit weight, r2 , or scale, r, is ®rst required. In fact, the initial robust estimate of r can be found as (Hampel et al. 1986, p. 105) …17† r^0 ˆ 1:483 medf v0 g i

or

Assuming that the derivative of q with respect to vi is w…vi †, then by derivation (Yang 1991; Koch and Yang 1998) the robust estimator of the model parameters is obtained as follows:

r^0 ˆ medf v0i g=0:6745 …17a† 0 0 where medf vi g denotes the median of all vi . It is easily understood that the estimate, r^0 , of Eq. (17) is most robust, because it has the highest breakdown point, 0.5 (Huber 1984; Rousseeuw and Leroy 1987, p. 9). The r^i can be expressed as follows, based on Eq. (17)

a^ ˆ …X T PX †ÿ1 X T PS

p r^i ˆ r^0 = pi

n X

pi q…vi † ˆ min

…11†

iˆ1

…12†

where P is an equivalent weight matrix; its ith diagonal element is pi ˆ pi w…vi †=vi

…13†

The error in¯uence function of the estimator of Eq. (12) is (see Yang 1991, 1994) IF …si ; a† ˆ ÿM ÿ1 xTj pj w…xj a ÿ sj †

…14†

M ˆ X T ZX

r^i is used as a criterion in the calculation of equivalent weight. The posterior covariance matrix of the estimated parameters is (Yang 1997) Ra^ ˆ …X T PX †ÿ1 r^2

…15†

where Z is a diagonal matrix, with elements zii ˆ pi Ew0 …xi a ÿ si † and Ew0 …xi a ÿ si † is the expectation of derivative of the w function. The choice of the w function or the equivalent weight function is theoretically based on the q function in Eq. (11) and, in turn, the determination of the q function depends on the actual distribution of observations. In practice, the distribution of the SLR observation is close to a normal one, even if some outliers exist. Therefore, the main components of the equivalent weight function should be equal to those of the LS case. Based on the calculation experiments, the equivalent weight function can be constructed with three segments (Yang 1991, 1994) as follows:

…19†

where r^2 is approximately expressed as (Yang 1997) r^2 ˆ

where

…18†

V T PV nÿm

…20†

These are the basic robust estimators for evaluating the range bias and time bias, and for accessing the polynomial coecients, which model the higher-frequency noise and the irregular systematic errors. From the whole procedure of the robust estimation, we can see that most of the outliers are downweights or zero weights at both steps. If Ds…tj †, for example, is an outlying observation, then its equivalent weight calculated by Eq. (16) will be small, and the residual u…tj † of Ds…tj † will obviously be apparent. In the second step, u…tj † becomes an observation. It is obvious that if Ds…tj † is downweighted, u…tj † is then outlying, and if u…tj † is large enough, then the equivalent weight of u…tj † will be small or close to zero. The outliers, therefore, cannot in¯uence the model parameters in the above data-®tting steps.

348

4 Practical procedure for robust estimation The robust method for estimating the station-dependent systematic errors discussed above is performed as an iterative procedure. The robust estimates of parameters are based on the equivalent weights, in turn dependent on the accuracy of the reference trajectory which is used to compute the laser-ranging residuals. If the residuals are not reliable, then the robust ®tting from the corresponding equivalent weights will not make any sense. Thus it is required that the procedure for obtaining the robust preliminaries with the highest breakdown point for the error equation (2) should be applied as follows. 1. Calculate the initial r0Ds r0Ds

ˆ medfjDs…ti †jg=0:6745

…21†

2. Reject all of those outlying observations, Ds…tj †, by the initial equivalent weight 8 < pDsi Ds…ti †=r0Ds  k0 0 ˆ …22† pDs i : Ds…ti †=r0 > k0 0 Ds where k0 is chosen to be 1.0 to 1.5. This is an alternative step. In a large observational sample, with or without this step the results will be nearly the same. In a small sample which is heavily contaminated by outliers, this step can improve the reliability of the preliminary estimates of the model parameters by ignoring the most suspicious observations, which in turn provides more reliable residuals for the following iterative estimation. 3. Calculate the initial model parameters a^0 and b^0 as well as the preliminary residuals, v0Dsi . 4. Perform the iterative computation procedure by using the equivalent weight function, Eq. (16). The computation procedure for the error equation (4), which is nearly the same as that for the model of Eq. (2), is omitted here. Analyzing the above robust ®tting procedures, we can make the following observations. 1. The initial estimate of the scale, calculated at the ®rst step based on Eq. (21), can guarantee the scales of Ds…ti † and u…ti † to be the most robust since the breakdown point of 0.5 allows the nearly 50% of outliers existing in the observations but having no in¯uence on the estimate. 2. On the basis of the robust scale, all Ds…ti † in Eq. (1) and u…ti † in Eq. (3) which are larger than k0 r0 will be rejected. At this step, we do not consider the eciency of the estimates of the model parameters, only guarantee the robustness of them. 3. Based on the reliable model parameters, reliable observation residuals can be obtained. 4. The reliable residuals guarantee that the equivalent weight will also be reliable. In turn, in relation to the

reliable equivalent weight, the outlying observations will have small equivalent weights or zero weights. Therefore, the outliers will not e€ect the estimates of the model parameters. 5. In the iterative calculation procedure the eciency can be compensated for, since the iterative computation is based on the equivalent weight function Eq. (16). All the observations will be used in the calculation, including those observations that have not been used at the ®rst step.

5 Results and comparison Laser-ranging data sets for the year 1994, for the geodetic satellite Starlette, were analyzed by using the JGM-3 gravity model and other state-of-the art force models. The average RMS for a ®ve-day orbit ®t is 5 cm. The obtained laser-ranging residuals for these data sets were further analyzed using the conventional LS technique and the robust estimation method presented in previous sections. The averaged ranging bias and time bias as well as their standard derivation and RMS for the 31 stations tracking Starlette in 1994 are listed in Table 1. The goals of analyzing the laser-ranging residuals are to pick out the outliers based on a selected rejection criterion, and to compute the ranging bias, time bias and polynomial RMS for each pass, as well as to obtain an estimate of the station-dependent systematic errors. However, the conventional LS technique has a strong masking ability: the errors for some outlying observations were apparently reduced in the corresponding residuals, and passed through the selected rejection criterion. Such mis-retained observations will distort the estimates of the ranging bias and time bias. Table 1 shows that the standard deviations of ranging bias and time bias for 31 tracking stations in 1994 are reduced by the robust estimation by as much as 31 and 70%, respectively. The small deviations of ranging bias and time bias, and the small polynomial RMS, should ensure the reliability of the robust solution in the analysis of satellite laser-ranging residuals. The improvements are due to the fact that the robust technique can eciently isolate the outliers and resist their in¯uence on the solution. Thus, the results of the station-dependent systematic errors are more reliable when using robust estimation. However, in conclusion, it should be pointed out that the robustness and eciency of the robust estimates of the systematic errors in SLR data are signi®cantly improved in large observation samples. For the small sample, the results of robust estimation will not be much better than those of the LS technique, as shown in Table 1. When the errors of observation are well distributed, the results obtained from LS and robust estimation are almost equal.

349 Table 1. Comparison of ranging bias and time bias (with their standard derivations) for the 31 stations tracking Starlette during 1994 using robust and LS methods ID

Robust

LS

Passes

R-bias (cm)

rB (cm)

T-bias (ms)

rTB (ms)

1873 1884 1953 7080 7090 7097 7105 7109 7110 7210 7295 7403 7411 7530 7541 7548 7805 7810 7811 7824 7831 7835 7836 7838 7839 7840 7843 7882 7883 7939 8834

3 5 50 139 431 25 341 373 411 288 10 182 3 69 7 4 15 95 23 25 64 154 83 49 111 307 163 26 5 184 145

)17.26 7.82 12.05 )1.04 )0.14 )1.75 0.68 )0.96 )1.72 )1.76 )3.00 0.40 )3.33 1.33 )2.09 )5.68 8.67 0.31 )2.92 )0.60 )0.23 )0.06 1.89 )1.05 0.66 )0.16 )1.78 )1.89 )2.73 1.23 )2.28

2.11 3.27 2.38 0.32 0.21 0.81 0.20 0.17 0.19 0.32 2.00 0.31 0.15 0.58 0.63 5.35 3.93 0.38 0.73 0.51 0.52 0.20 0.36 0.72 0.29 0.20 0.35 1.10 2.35 0.57 0.37

0.73 4.64 0.71 0.19 0.31 0.16 0.01 )0.06 )0.25 0.14 )0.24 2.51 1.38 )0.04 )1.63 )2.95 )13.97 )0.31 )2.25 1.03 0.51 )0.01 )0.14 )1.16 0.45 0.65 )0.29 )2.92 5.76 )1.21 )0.30

0.15 2.06 4.24 0.34 0.34 0.71 0.24 0.18 0.16 0.26 0.77 0.60 0.53 0.43 3.60 5.03 3.20 0.27 0.86 0.77 0.34 0.11 0.40 0.70 0.21 0.27 0.53 1.54 2.57 0.50 0.24

3 4 40 142 436 29 322 377 415 298 10 190 3 67 8 5 11 89 25 24 65 171 79 50 120 305 164 26 6 177 185

Totals

3790

)0.38

0.20

0.02

0.08

3846

Acknowledgements. This project is sponsored by the National Science Foundation (No. 49825170) and the National Climbing Program of China as well as the Scienti®c Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry.

References Hampel FR, Ronchetti EM, Rousseeuw PJ, Stahel WA (1986) Robust Statistics. The Approach Based on In¯uence Functions. John Wiley, New York Huber PJ (1964) Robust estimation of a location parameter. Ann Math Statist 35: 73±101 Huber PJ (1984) Finite sample breakdown point of M- and Pestimators. Ann Statist 12: 119±126 Koch KR, Yang Y (1998) Robust Kalman ®lter for rank de®cient observation model. J Geod 72: 436±441 Kubik K (1982) An error theory for Danish method. International Symposium on Numerical Methods for Photogrammetric Mapping, Helsinki, 1982, Helsinki Rousseeuw PJ, Leroy AM (1987) Robust regression and outlier detection. John Wiley, New York

Passes

R-bias (cm)

rB (cm) T-bias (ms)

rTB (ms)

)13.54 )9.27 8.03 )0.96 )0.18 )1.63 0.01 )0.95 )1.68 )1.78 )7.74 )0.31 )5.44 )0.06 )0.04 )1.17 6.22 0.78 )2.33 0.12 )1.05 )0.67 1.32 )1.96 )0.19 )0.46 2.00 )0.73 )3.84 1.00 0.33

2.54 2.59 2.45 0.53 0.27 1.16 0.35 0.22 0.28 0.47 4.02 0.56 1.28 1.12 3.28 2.40 3.30 0.52 1.33 0.84 0.82 0.32 0.50 0.83 0.39 0.28 0.53 1.53 1.95 0.70 0.55

19.40 )6.55 )11.12 )3.36 0.76 3.31 0.79 0.37 )1.71 )0.53 1.17 1.31 2.83 2.56 1.12 )37.43 )41.90 1.10 6.26 10.11 )2.84 1.75 2.14 )3.43 0.51 1.50 0.00 )5.97 )4.28 0.98 0.12

7.80 5.60 8.44 1.13 0.81 2.88 0.81 0.56 0.57 0.79 4.38 1.51 2.29 1.75 4.91 11.43 5.89 1.27 3.21 2.49 2.07 0.69 1.32 2.35 1.13 0.74 1.71 3.58 6.60 1.99 1.24

)0.53

0.11

0.13

0.26

Scha€rin B (1985) On robust collocation. In: Sanso F (ed) Proc 1st Marussi Symp Math Geod, Milano, 1986, pp. 343±361 Schutz BE, Cheng MK, Eanes RJ, Shum CK, Tapley BD (1994) Geodynamic results from Starlette orbit analysis. In: Smith DE, Yurcotte AL (eds) Contributions of space geodesy to geodynamics: Earth dynamics, geodynamics series vol 24, pp. 175±190 Tapley BD, Schutz BE, Eanes RJ, Ries JC, Watkins MM (1994) Lageos laser ranging contributions to geodynamics, geodesy, and orbital dynamics. In: Smith DE, Yurcotte AL (eds) Contributions of space geodesy to geodynamics: Earth dynamics, Geodynamics series, vol 24, pp. 147±174 Yang YX (1991) Robust Bayesian estimation. Bull Geod 65: 145±150 Yang YX (1992) Robustifying collocation. Manuscr Geod 66: 21±28 Yang YX (1994) Robust estimation for dependent observations. Manuscr Geod 19: 10±17 Yang YX (1997) Estimators of covariance matrix at robust estimation based on in¯uence functions. Z Vermess 122: 166±174 Zhou J (1989) Classical theory of errors and robust estimation. Acta Geod Cartogr Sin 18(2): 115±120