Robust Geodetic Parameter Estimation

under Least Squares through Weighting on

the Basis of the Mean Square Error
Francis W. O. Aduol

Abstract

A technique for the robust estimation of geodetic parameters under the least squares method when
weights are specified through the use of the mean square error is presented. The mean square error is
considered in the specification of observational weights instead of the conventional approach based on
the observational variance. The practical application of the proposed approach is demonstrated
through computational examples based on a geodetic network. The results indicate that the least
squares estimation with observational weights based on the mean square error is relatively robust
against outliers in the observational set, provided the network (or the system) under consideration has
a good level of reliability, as to make the network (or system) stable under estimation.

Introduction

The classical approach in the estimation of geodetic parameters is through the least squares method
within the framework of the Gauss-Markov model given as:

~
y = Ax + ε ; ε ~ (0, σ 20W −1 ) (1)

where y is an n × 1 vector of observations, A is an n × m design matrix, x is an m × 1 vector of

unknown parameters, ε is an n × 1 vector of observational errors, σ 20 is the variance of unit weight,
and W is an n × n positive-definite weight matrix.

This estimation model is based on the assumption that the observational errors collected in the vector
ε occur randomly and are distributed according to the normal distribution. with this assumption and
under the least squares condition that ε TWε be minimum, the following estimates may be obtained:

x = ( ATWA) −1 ATW~ y (2)

D( x ) = Σ xx  = σ 0 ( ATWA) −1
2

y = Ax = A( ATWA) −1 ATWy~
D( y ) = Σ yy  = AΣ xx  AT = σ 20 A( ATWA) −1 AT
ε = y − y = y − Ax = y − A( ATWA) −1 ATWy
D(ε ) = Σ εε  = Σ ~~
yy − AΣ xx


σ 20 = ε TWε / (n − m)

In the event however that the observational vector y may be contaminated with a bias parameter b
(whereby the bias may be as a result of gross errors, or systematic errors, or a combination of both),

5
then the assumption ε ~ (0, σ 20W −1 ) gets invalidated, in that the errors on y , which now also com-
prise b can no longer be considered to be distributed according to the normal distribution. The conse-
quence of this is that if the estimation of the unknown parameters still be performed according to the
least squares condition, under the Gauss-Markov model as in (1), then the so obtained estimates will
be biased as a result of b . To deal with this problem, two options come into consideration: (i) one
performs the estimation under the least squares under the model (1) but seeks to identify and remove
outliers (biased observations) from the observational dataset in what we may refer to as outlier detec-
tion, or (ii) one adopts estimation techniques that are robust with respect to the biases under robust
estimation.

The propagation of outlier detection in geodesy and surveying was motivated by the works of W
Baarda [2, 3, 4]. Today outlier isolation forms an integral component of any major geodetic data proc-
essing and analysis. However the detection and isolation of outliers within the framework of the
Gauss-Markov model as specified in (1), still suffers from the tendency of the ordinary least squares
method to spread out the effect of outliers among observations, thereby rendering the isolation of the
outliers difficult, and sometimes altogether impossible. To cope with this problem, robust estimation
techniques offer real alternatives.

The objective in robust estimation is to perform an estimation of the parameters from the observations
in such a way that the estimates of the parameters so obtained are virtually unaffected by any biases or
outliers that may be present in the observations. An extensive study of the application of robust esti-
mation in geodesy is reported in [5]. Robust estimation techniques in the estimation of parameters in
general were however brought to the fore through the works of P J Huber [8, 9, 10], while a further
extensive treatment of the subject has been presented by [7]. The core of Huber’s technique is the M-
estimator, which is based on the maximum-likelihood method.

A general characteristic of the robust estimation techniques is that they restrict a range of observa-
tional error within which the observations may be accepted, and observations associated with observa-
tional error outside the specified range are ‘cut off‘ from the estimation process within the process of
‘winsorisation’. The problem in this approach however is that the decision on where the ‘cut-off’ point
itself should be is rather subjective. As an alternative approach in robust estimation, a procedure for
robust estimation based on iterative weighting of observations was suggested in [1]. This was an at-
tempt at a procedure that would avoid excluding any observations from the estimation procedure, but
include all the observations within the estimation procedure except with appropriate weighting.

In this presentation, we extend the concept of iterative weighting by considering it from the point of
view of the observational weights based on the mean square error (MSE), and evaluate the effective-
ness of the method through the computation of a practical network.

The Mean Square Error

Let us consider a parameter vector ξ , whose realisation (obtained through estimation or otherwise) is
ξ , then the mean square error of ξ is given as
M (ξ ) = E[(ξ − ξ )(ξ − ξ ) T ] (3)
In general we have that E(ξ ) = ξ + β where β is a bias vector. Thus we may rewrite (3) as

M (ξ ) = E[(ξ − ( E (ξ ) − β ))(ξ − ( E (ξ ) − β )) T ]

= E[(ξ − E (ξ ))(ξ T − E (ξ ) T )] + ββ T (4)

6
But we have that the dispersion D(ξ ) of ξ is given as

D(ξ ) = E [(ξ − E (ξ ))(ξ T − E (ξ ) T )] (5)

Thus
M (ξ ) = D(ξ ) + ββ T (6)
(see e.g. [11] and [6]).
In the special case that β = 0 , we have then that

M (ξ ) = D(ξ ) (7)

From the fact that the mean square error incorporates the biases in the realisation of a parameter, the
mean square error is a much more effective and efficient estimate of the quality of the parameter in the
sense of accuracy. The dispersion on the other hand, respectively the variance, as is ordinarily known,
gives the precision of the estimate or realisation, which however only becomes also a measure of accu-
racy in the special case when β = 0 , in which case (7) obtains.

We have from (6) that in the special case that it is a single independent parameter being considered,
the mean square error is given as: mean-square-error = variance+bias² .

The Estimation Model

In the event that the observation ~

y in (1) is contaminated with a bias b , then we have that

E(~
y) = y +b (8)

where y is the ’true’ value of the parameter.

But we have that
~
y = E(~
y ) + ε ; ε ~ (0, Σ ~~
yy ) (9)
Then with (8) and (9) we have
~
y = y +b+ε (10)
which with v:= b + ε , becomes
~
y = y+v (11)
For
y = Ax (12)
we then have that ~
y = Ax + b + ε or

~
y = Ax + v , E (v ) = E (b + ε ) = b, M ( ~
y ) = Σ ~~
yy + bb
T
(13)

We adopt this as the model for the estimation of the parameters within the framework of least squares.

We note therefore from (13) that if we can estimate v such that

E (v ) = E (b + ε ) = b + E (ε ) = b , (14)
then we would have been able to obtain an unbiased estimate of x that is relatively free from the in-
fluence of the bias b .

In the conventional least squares approach, whereby the model is defined according to (1), if the
model had a bias parameter as to be described according to (13), but with the stochastic part described

7
through ε ~ ( 0, Σ ~~
yy = σ 0W ) , then the model would have not been appropriately specified, so that
2 −1

the parameters estimated with the model will be biased. We seek to overcome the bias effect in that we
define the estimation model through (14) and weight the observations according to the mean square
error (MSE), which already incorporates the bias effect. We propose then to define the weight W of
the observations as
W = σ 20 M ~~
−1
yy (15)
~
yy = M ( y ) .
in which we have taken M ~~

If we assume independence of observations ~

yi (i = 1,... n) , then we have that for an observation ~y,
the mean square error may be given as
mi = σ i2 + bi2 (16)

for σ i2 and bi being respectively the variance and bias of ~

yi . Then the weight of ~
yi can now be de-
fined as
σ 20
wi = (17)
m
−1
yy will exist due to the fact that M ~~
With the weights so defined, we notice that M ~~ yy has been taken to

be a diagonal matrix, and hence W according to (15) can be evaluated.

The question however is how does one evaluate the mean square error in the first place, when the bias
b itself is in the first instance unknown, and must in any case be evaluated. We seek to deal with this
problem in that we evaluate b iteratively and hence also W .

The Estimation Process

We begin the estimation process by assuming nominally that b = 0 . With this, we notice that we will
simply be having the Gauss-Markov model as described in (1). From this, the first estimates of b as
‘residuals’ will have been obtained. With the residuals vi , a new value for mi is obtained according to
(16), however with σ i being as originally set, since these are the original variances of the observa-
tions, which are assumed known a priori. With the new mean square error values, the estimation proc-
ess is repeated. The process is repeated until convergence for the estimated parameters is achieved at
the specified level of tolerance. In particular, since the main parameters being estimated are the un-
known parameter vector x , the convergence of the x parameters would be more appropriately
adopted as control for the iteration.

Through the iterative process, the mean square error of an observation is estimated for simultaneously
as well and consequently the mean-square-error weight of the observations. The robustness of the
procedure is thus contained in the mean-square-error weight, which is a much more comprehensive
and realistic representation of the observational weights.

8
The Test Example

The test network

A two-dimensional network as shown in Fig. 1 was adopted for the test example. The network com-
prises 9 points, which are linked by distance observations. A single distance observation was consid-
ered to have been measured with a standard error of 3mm+0.5ppm; with this the eventual standard
error for the mean distance adopted was then deduced from the number of individual measurements
from which the particular mean distance is obtained. The network has a total of 30 distance measure-
ments.

Experimental design

Four versions of the network were computed; these were designated as Net-0, Net-1, Net-2, and Net-4.
The networks were specified according to the numbers of gross errors they contained as follows: Net-0
- no gross errors; Net-1 - one gross error; Net-2 - two gross errors; and Net-4 - four gross errors. The
gross errors were simulated into the networks as given in Table 1.

Table 1: The simulated gross errors

Line Error [metres] Network
4-7 +0,780 1,2,4
2-3 -5,067 2,4
7-11 +0,355 4
3-4 -0,055 4

Each version of the network was then computed on the basis of both the ordinary least
squares and the least squares method with mean-square-error weights as proposed here. The
network was computed throughout in free-network mode.

Results

In the results presented below, X and Y are estimated point coordinates in metres; σ X and
σ Y are estimated positional standard errors in metres; a and b are the major and minor axes
of the positional error ellipse in metres, while ϕ is the orientation of the major axis of the
ellipse in degrees taken with respect to the X axis.

Net-0

Table 2: Conventional least squares

Point X Y σX σY a b ϕ
1 5428972.186 3462429.367 0.0044 0.0059 0.0063 0.0039 155.4
2 5439065.854 3468259.524 0.0044 0.0053 0.0056 0.0040 150.5
3 5457025.454 3476522.257 0.0054 0.0059 0.0059 0.0054 170.3
4 5465079.259 3433374.141 0.0045 0.0064 0.0065 0.0042 163.0
5 5448374.040 3427727.520 0.0039 0.0043 0.0047 0.0034 143.1
6 5439527.166 3423319.106 0.0053 0.0066 0.0076 0.0036 145.4
7 5447601.042 3443324.505 0.0060 0.0038 0.0060 0.0037 9.8
8 5411104.687 3454335.360 0.0064 0.0092 0.0097 0.0056 158.0
11 5464986.157 3457965.548 0.0086 0.0053 0.0086 0.0053 179.8

Table 3: Robustified least squares (3 iterations)

9
 Point X Y σX σY a b ϕ
1 5428972.186 3462429.370 0.0030 0.0042 0.0045 0.0026 155.4
2 5439065.857 3468259.522 0.0030 0.0043 0.0045 0.0026 157.1
3 5457025.455 3476522.260 0.0036 0.0040 0.0040 0.0036 3.4
4 5465079.258 3433374.142 0.0030 0.0041 0.0043 0.0028 161.1
5 5448374.040 3427727.518 0.0027 0.0032 0.0034 0.0024 146.8
6 5439527.168 3423319.103 0.0036 0.0049 0.0055 0.0025 147.8
7 5447601.039 3443324.504 0.0041 0.0028 0.0042 0.0027 12.9
8 5411104.688 3454335.360 0.0042 0.0060 0.0064 0.0036 156.0
11 5464986.155 3457965.548 0.0057 0.0035 0.0057 0.0035 2.1

Observations treated as containing gross errors in the adjustment: Nil

Net-1

Table 4: Conventional least squares

Point X Y σX σY a b ϕ
1 5428972.119 3462429.461 0.0616 0.0834 0.0882 0.0545 155.4
2 5439065.822 3468259.569 0.0625 0.0745 0.0796 0.0559 150.5
3 5457025.475 3476522.203 0.0767 0.0828 0.0830 0.0765 170.3
4 5465079.471 3433373.962 0.0631 0.0899 0.0922 0.0597 163.0
5 5448374.111 3427727.603 0.0556 0.0606 0.0666 0.0484 143.1
6 5439527.284 3423319.082 0.0741 0.0932 0.1078 0.0505 145.4
7 5447600.816 3443324.577 0.0841 0.0537 0.0848 0.0525 9.8
8 5411104.631 3454335.397 0.0898 0.1302 0.1367 0.0795 158.0
11 5464986.118 3457965.474 0.1209 0.0752 0.1209 0.0752 179.8

Table 5: Robustified least squares (6 iterations)

Point X Y σX σY a b ϕ
1 5428972.189 3462429.366 0.0033 0.0047 0.0051 0.0027 153.8
2 5439065.857 3468259.524 0.0031 0.0042 0.0044 0.0028 156.4
3 5457025.454 3476522.262 0.0038 0.0043 0.0043 0.0038 178.6
4 5465079.254 3433374.146 0.0042 0.0049 0.0056 0.0032 143.1
5 5448374.039 3427727.516 0.0031 0.0035 0.0036 0.0029 153.1
6 5439527.165 3423319.103 0.0043 0.0053 0.0061 0.0029 144.4
7 5447601.044 3443324.502 0.0055 0.0030 0.0055 0.0030 177.8
8 5411104.690 3454335.360 0.0045 0.0065 0.0069 0.0038 155.0
11 5464986.154 3457965.550 0.0061 0.0038 0.0061 0.0038 4.8

Observations treated as containing gross errors in the adjustment

Line Gross error as isolated Redundancy

(metres)

4–7 -0.7952 100%

10
Net-2

Table 6: Conventional least squares

Point X Y σX σY a b ϕ
1 5428972.886 3462429.822 0.4272 0.5778 0.6114 0.3775 155.4
2 5439067.003 3468259.524 0.4332 0.5164 0.5515 0.3876 150.5
3 5457022.867 3476521.657 0.5314 0.5738 0.5750 0.5301 170.3
4 5465079.528 3433373.759 0.4374 0.6229 0.6391 0.4135 163.0
5 5448374.140 3427727.373 0.3855 0.4203 0.4615 0.3351 143.1
6 5439527.264 3423318.920 0.5135 0.6457 0.7470 0.3502 145.4
7 5447600.915 3443324.313 0.5827 0.3720 0.5880 0.3637 9.8
8 5411105.229 3454336.092 0.6223 0.9021 0.9475 0.5507 158.0
11 5464986.014 3457965.650 0.8379 0.5210 0.8379 0.5210 179.8

Table 7: Robustified least squares (8 iterations)

Point X Y σX σY a b ϕ
1 5428972.187 3462429.365 0.0037 0.0049 0.0052 0.0032 153.5
2 5439065.854 3468259.524 0.0038 0.0043 0.0045 0.0035 147.7
3 5457025.460 3476522.263 0.0061 0.0045 0.0062 0.0043 14.5
4 5465079.254 3433374.146 0.0043 0.0051 0.0058 0.0033 144.4
5 5448374.039 3427727.516 0.0030 0.0035 0.0036 0.0029 157.4
6 5439527.165 3423319.103 0.0043 0.0052 0.0061 0.0028 143.5
7 5447601.044 3443324.503 0.0054 0.0030 0.0054 0.0030 176.4
8 5411104.689 3454335.358 0.0046 0.0065 0.0068 0.0042 158.9
11 5464986.154 3457965.550 0.0061 0.0039 0.0061 0.0039 0.5

Observations treated as containing gross errors in the adjustment

Line Gross error as isolated Redundancy

(metres)

4–7 -0.7951 100%

2–3 +5.0783 100%

Net-4

Table 8: Conventional least squares

Point X Y σX σY a b ϕ
1 5428972.864 3462429.832 0.4259 0.5760 0.6096 0.3764 155.4
2 5439066.998 3468259.726 0.4319 0.5149 0.5498 0.3864 150.5
3 5457022.859 3476521.637 0.5298 0.5721 0.5733 0.5284 170.3
4 5465079.501 3433373.799 0.4361 0.6210 0.6371 0.4122 163.0
5 5448374.131 3427727.363 0.3843 0.4190 0.4601 0.3341 143.1
6 5439527.268 3423318.877 0.5119 0.6437 0.7447 0.3491 145.4
7 5447600.832 3443324.287 0.5810 0.3709 0.5862 0.3625 9.8
8 5411105.206 3454336.094 0.6204 0.8993 0.9446 0.5491 158.0
11 5464986.187 3457965.713 0.8354 0.5194 0.8354 0.5194 179.8

Table 9: Robustified least squares (8 iterations)

11
 Point X Y σX σY a b ϕ
1 5428972.136 3462429.354 0.0061 0.0066 0.0073 0.0052 142.2
2 5439065.800 3468259.517 0.0071 0.0058 0.0077 0.0049 148.5
3 5457025.431 3476522.214 0.0085 0.0089 0.0091 0.0083 29.9
4 5465079.216 3433374.153 0.0060 0.0082 0.0088 0.0051 152.9
5 5448374.005 3427727.512 0.0042 0.0060 0.0060 0.0042 2.0
6 5439527.134 3423319.092 0.0061 0.0074 0.0085 0.0046 145.3
7 5447600.999 3443324.497 0.0077 0.0055 0.0077 0.0055 176.9
8 5411104.642 3454335.337 0.0065 0.0085 0.0086 0.0063 163.0
11 5464986.482 3457965.652 0.0123 0.0083 0.0124 0.0081 167.4

Observations treated as containing gross errors in the adjustment

Line Gross error as isolated Redundancy

(metres)

4–7 -0.7951 100%

2–3 +5.0783 100%
2 –11 +0.3184 99%
5 –11 +0.2669 100%

12
13
Discussion

The technique described here, like most techniques for robust estimation and management of outliers
in observations, depends considerably for is effectiveness on the reliability of the network. The tech-
nique is only able to isolate outliers and damp their effects on the estimation process through the fact
that the bias-free observations are in a position to estimate effectively the unknown parameters and at
the same time resist the influences from the outlying observations. This way, the effects of the outly-
ing observations on the estimated parameters are rendered minimal.

If however the bias-free observations should be overwhelmed by the outlying observations, either
through sheer numbers or through geometric distribution within the observational set, then an ade-
quate solution of the estimates may be rendered difficult, or altogether impossible. For instance, in the
present study, in the case with four gross errors in the network a converging solution was only ob-
tained after eight iterations. However, although the results indicate that the estimated parameters have
been obtained with relatively acceptable precision, the space of convergence of the parameters is bi-
ased, as can be ascertained through comparing the results in Table 9 with those in Table 3. This bias
has been caused by the fact that the network was not sufficiently robust in configuration (i.e. in ge-
ometry, as well as observational type, number and quality) as to able to isolate the observations con-
taining gross errors, which in the first place were rather ‘unsuitably’ distributed. The gross errors were
here distributed such that out of the five network points, 2,3,4,7,11, connected with gross-error-
contaminated observations three of the points, namely 3,4,7, were each connected with gross-error-
contaminated observations. The result of this was that the gross errors in lines 3-4 and 7-11 could not
adequately be isolated, and instead lines 2-11 and 5-11 were interpreted as the ones containing the
gross errors.

In the cases with one and three gross errors, whose results are presented in Tables 5 and 7, the biases
were effectively isolated, even though in this case point 7 was still connected by two gross-error-
contaminated observations. The results for these two cases were found to be even more precise than

14
those from the ordinary gross-error-free least squares case presented in Table 2. In the initial case with
no gross errors we notice from Table 3 that the results for the robustified least squares technique are
considerably more precise than the case with ordinary least squares. Thus we have that even with ob-
servations that are effectively gross-error free one obtains more efficient estimates than with the ordi-
nary least squares approach.

Conclusion

The results of this study demonstrate that the definition of the observational weights through the mean
square error results in robustified least squares estimates. The technique tested was able to cope effec-
tively with outliers in the observational set. The effectiveness of the technique however, as can be
expected, is dependent on the reliability of the network, and especially on the particular observations
contaminated with outliers. When the network reliability is sufficiently high, the technique of weight-
ing observations on the basis of the mean square error instead of the variance can be relied on to yield
fairly reliable estimates even with gross errors in the observational set. The computational process is
rendered rather slower than in the case of weights based on variances, due to the fact that the mean
square error has essentially to be determined iteratively.

Acknowledgement

This work was completed when the author was a visiting DAAD (German Academic Exchange Ser-
vice) scholar at the Geodätisches Institut of the University Karlsruhe, with Prof Günter Schmitt as his
host. The author is grateful to DAAD for this support and to Prof Schmitt for kindly hosting him while
at the Geodätisches Institut. The test network adopted for the study was also kindly made available to
the author by the Geodätisches Institut, Karlsruhe, and this assistance is also acknowledged. The au-
thor also wishes to thank Dr S M Musyoka, Department of Surveying, University of Nairobi, for his
assistance in the preparation of the diagrams.

References

1. Aduol, F.W.O., 1994. Robust geodetic parameter estimation through iterative weighting. Survey
Review, 32, 252: 359-367.
2. Baarda, W., 1967. Statistical concepts in geodesy. Netherlands Geodetic Commission, publica-
tions on Geodesy, New Series, Vol. 2, No. 4, Delft.
3. Baarda, W., 1968a. Statistics - a compass for the land surveyor. Computing Centre of the Delft
Geodetic Institute.
4. Baarda, W., 1968b. A testing procedure for use in geodetic networks. Netherlands Geodetic
Commission, publications on Geodesy, New Series, Vol. 2, No. 5, Delft.
5. Borutta, H. 1988. Robuste Schätzverfahren für geodätische Anwendungen. Schriftenreihe Stu-
diengang Vermessungswesen Universität der Bundeswehr München, Heft 33. München.
6. Grafarend, E.W., Schaffrin, B., 1993. Ausgleichungsrechnung in linearen Modellen. Wissen-
schaftsverlag, Mannheim. Pp. 116 - 117.
7. Hampel, F.R., Ronchetti, E.M., Rousseeuw, P., and Stahel, W.A., 1986. Robust Statistics - the
Approach based on Influence Functions. John Wiley & Sons, New York.
8. Huber, P.J., 1964. Robust estimation of a location parameter. Annals of Mathematical Statistics,
35: 73-101.
9. Huber, P.J., 1972. Robust statistics - A review. Annals of Mathematical Statistics, 43: 1041-1067.
10. Huber, P.J., 1981. Robust Statistics. John Wiley & Sons, New York.
11. Toutenburg, H., 1992. Lineare Modelle. Physica Verlag, Heidelberg. Pp 35 - 36.

15
16