. Problem of defects in geodetic network

The numerical elaboration of a geodetic network usually consists of two stages. At first, the initial diagnostics of data is performed, which allows detection and removal of various types of observational errors and geometric defects. Secondly, the final elaboration of the network is performed, including adjustment of observations and calculation of the coordinates of points. Defects in a geodetic network, that is, situations requiring identification and elimination, can be divided into several types, as shown in Figure 1.

Figure 1

Defects in geodetic network: outliers and geometric (datum or configuration) defects.

https://rgg.edu.pl/f/fulltexts/173249/j_rgg-2021-0004_fig_001_min.jpg

A metric defect in observations (outliers) means that the absolute value of the observational correction (error estimator) exceeds a certain threshold level (depending on the adopted probabilistic model and the significance level). Statistical methods are used to identify such errors, but the correct result depends on the network's local reliability (local redundance) (Baarda, 1968; Prószyński and Kwaśniak, 2019). The reliability of a geodetic network is related to the conditioning of the normal matrix in the problem of adjustment of the geodetic network using the least-squares method. In general, the point here is that in ill-conditioned systems, detection of a possible metric defect (outlier) can be difficult or impossible. The difficulty in identifying outliers is also due to the principle of least squares because this method causes some kind of ‘spilling’ of errors over the network elements located in a specific vicinity of the outlier observation. In addition, the identification task becomes more difficult as the size of the observational system increases. Therefore, other rules for observation adjustment are formulated, and the so-called robust estimation, in contrast to the averaging feature of the least-squares method, causes some differences between observations. These replace the least squares function by an objective function with spline components (Hampel (1971); Huber (1964) and many others’ proposals) or based on regular, but non-square components of the objective function (Holland and Welsch, 1977; Kadaj, 1978, 1988). A comprehensive review of the robust estimation methods can be found in the book by Wiśniewski (2013).

The outlier phenomena may be related to the coordinates of the tie points. A straightforward method of detecting such errors was presented by Prof. Hausbrandt (Hausbrandt, 1954; Janusz, 1958). It assumes that the coordinate reference points are treated as pseudo-observations with specific mean errors, estimated on the basis of information on the accuracy of the determination of points in the network of the appropriate class. The magnitude of the correction to such pseudo-observation indicates a possible existence of an outlier. The principles of robust estimations can also be used for the same purpose, although the Hausbrandt method seems to be sufficiently effective.

This article is not a complete dissertation on the above subject. We will focus mainly on the problem of detecting configuration defects, which – as shown by practical experience – are a very common phenomenon and difficult to identify, especially in large data sets.

The geometric (topology) defect (u = uref + uobs) can be classified into two types, namely datum defect (missing of reference elements) and configuration defects (network observational structure defects). The first type implies that a certain number uref of elementary data is missing to determine the position of the object (network) in the coordinate system. It may be a situation intended to perform the so-called free network adjustment (Mittermayer, 1972) and may also be the result of an unintentional mistake. A configuration defect means the lack of a certain amount (uobs) of observations for network determination. The reason may not only be the omission of certain observations, but also numbering errors of points in measurement plans. While the reference defect is easily identifiable, as it is a matter of several parameters, the configuration defect is usually difficult to spot, especially in the case of geodetic networks of considerable size in terms of the number of points and observations. The main topic of this publication is the ability to identify network configuration defects automatically.

. The functional and stochastic model of a geodetic network

The functional model of a geodetic network can be presented in the form of a vector equation (col = single-column matrix):

(1)
L+V=F(X),
where:

  • L = col[Li : i = 1, 2, . . . , m] is an observational vector, with some covariance matrix C,

  • P = C−1 is the weight matrix,

  • V = col[Vi : i = 1, 2, . . . , m] is a vector of unknown corrections,

  • X = col[Xj : j = 1, 2, . . . , n] is a vector of unknown parameters (point coordinates),

  • F(X) = col[Fi(X) : i = 1, 2, . . . , m] is a vector function of the vector X.

Let X0 denote the vector of approximate coordinates – the approximation of the vector X. The linearized Equation (1) has the form:

(2)
v=A⋅x−l,
where:

  • A = (dF/dX)(0) = [aij](m×n); (m ≥ n); aij = (dFi/dXj)(0)

  • x = X − X0; l = L − F(X0)

  • v = vector of observational corrections in the linearized system (2).

The least-squares solution for Equation (2) is the minimization of the function:

(3)
Φ(x)=vT⋅P⋅v=(A⋅x−l)T⋅P⋅(A⋅x−l).

From the necessary condition for the extreme, we get the normal equation:

(4)
AT⋅P⋅A⋅x=AT⋅P⋅l   or   B⋅x=w,
where:
(5)
B=AT⋅P⋅A: is the (n×n) symmetric matrix, rank(B)−n−u,u≥0,w=AT⋅P⋅l:   is the (n×1)vector.

If u = uref + uobs > 0, then the normal equation is singular, generating a linear subspace of solutions. In particular, if uref > 0 (datum defect) and the value of the defect is known, while uobs = 0, then it is possible to apply ‘free adjustment’ (e.g. Mittermayer (1972)) using the so-called pseudo-inversion (generalized Moore–Penrose inversion; see more details in Bjerhammar (1973); Moore (1920); Penrose (1955)). An important feature of such a situation is that the size of the defect is usually known in advance. A completely different case occurs when a defect is unknown, both in terms of its value and its situational location in the network. We only know that it exists because the matrix of the system of equations is singular. Detection and localization of such a defect are the main issues of this work.

. Least-squares adjustment in case if a defect is known

In the case of a non-zero and known defect, an unequivocal solution of the normal equation (4) may take place after applying an additional condition (conditions) to the vector of unknowns x. The universal condition in applications has the following form:

(6)
xT⋅x=min.

Then the unequivocal solution is written in matrix notation using the so-called generalized inversion (pseudo-inversion or Moore–Penrose inversion; see Moore (1920); Penrose (1955)):

(7)
x=B†⋅w.

The B† pseudo-inversion, where the solution (7) satisfies the condition (6), is defined as follows (Deutsch, 1965). Let M, N be matrices of the size (r × n) and rank(M) = rank(N) = r and satisfy the condition B = MT · N. Then

(8)
B†=NT⋅(N⋅NT)−1⋅(M⋅MT)−1⋅M.

In particular, if we assume M = N = R then M = RT · R and

(9)
B†=RT⋅(R⋅RT)−1⋅(R⋅RT)−1⋅R.

The above condition is fulfilled by the R matrix resulting from the Cholesky–Banachiewicz factorization of the B matrix. It has a ‘trapezoidal’ form (Figure 2).

Figure 2

Structure of the R matrix (0 – elements assumed to be zero, X – other elements of the matrix), r – rank of the normal matrix B, u = n–r = defect of the normal matrix.

https://rgg.edu.pl/f/fulltexts/173249/j_rgg-2021-0004_fig_002_min.jpg

Note that in the case of zero defect (u = 0) B† = B−1, then formula (9) takes the form:

(10)
B†=RT⋅(RT)−1⋅R−1⋅(RT)−1⋅R−1⋅R==I⋅R−1⋅(RT)−1⋅I=(RT⋅R)−1=B−1.

Another way of determining pseudo-inversion is the spectral decomposition of the normal matrix. This matrix is symmetric and non-negative definite, so the spectral decomposition takes the form:

(11)
B=S⋅Λ⋅ST;   Λ=diag(λ1,λ2,…,λr,o1,o2,…,ou),
where

  • 01 = 02 = . . . = 0u = 0,

  • S = matrix of eigenvectors as an orthonormal matrix: ST · S = I,

  • Λ = diagonal matrix of eigenvalues: λ1 ≥ λ2 ≥. . . ≥ λr > 0,

then
(12)
B†=ST⋅Λ†⋅S;   Λ†=diag(λ1−1,λ2−1,…,λr−1,o1,o2,…,ou).

Thus, using formula (8) or (9), we can adjust the network if the defect value is known. Typically, this situation is related to the datum defect and is intentional in the aspect of the so-called free network adjustment, consisting in the optimal fit of the network into a figure defined by approximate coordinates with simultaneous adjustment of the observations (Mittermayer, 1972). However, a troublesome problem may be if a non-zero defect of a network is unknown and is a result of mistakes in the network adjustment process. This will be mainly a network configuration defect (uobs > 0); however, missing observational data or errors in point numbering may also result in missing reference elements, that is, datum defect. Usually, such a situation is of an error nature, for example, as a result of accidental omission of certain observations in the data sets, necessary to obtain the effect of network determination, or as a result of errors of point identifiers in the observation plans of the geodetic network.

Similar to gross observation error, a network configuration defect will not be generally easy to identify and locate, especially in large networks with a variety of observation data. For this purpose, procedures using the exact pseudo-inversion formula (8) or (9) will not be (in principle) effective. First of all, they require knowledge of the defect value itself, not to mention its location in the network. Is it possible, then, to resolve the above issues in an ‘automatic’ manner?

The problem of detecting and locating configuration defects can be solved using the Tikhonov regularization method. This method was formulated to solve singular, linear algebraic equations (Tikhonov, 1963, 1965). It should be mentioned here that the problem of regularization also refers to the original works: Levenberg (1944); Marquardt (1963, 1970); Phillips (1962). Phillips’ work formulated the problem of regularization when solving integral equations. The work of Levenberg and the first of Marquardt's works concerned the iterative solution of the systems of non-linear equations as a method optimising the convergence of the iterative process. Another work by Marquardt (1970) already referred to the problems of regularization in the iterative process.

. Tikhonov procedure as alternative to generalized inversion and as a method for detecting configuration defects

General formula for Tikhonov regularization

Tikhonov regularization (Tikhonov, 1963, 1965) is based on the condition:

(13)
vT⋅P⋅v +α⋅ xT⋅x=min
from where:
(14)
(B+αI)⋅xα=w⇒xα=(B+αI)−1⋅w=Bα−1⋅w,Bα=B+αI;   I=unit matrix,
where α > 0 is a ‘small’ number on which the difference between solutions (14) and (7) depends.

If α decreases to zero, then solution (14) approaches (7) (Bjerhammar, 1973):

(15)
limα→0Bα−1⋅w=B†⋅w
but the Bα−1 matrix is not, however, an appropriate approximation of the B† pseudo-inversion. This approximation is the matrix
(16)
Qα=Bα−1⋅B⋅Bα−1=Bα−1⋅(I−α⋅Bα−1)=Bα−1−α⋅(Bα−1)2≈B†
based on a theorem, compatible with (15):
(17)
limα→0Qα=B†=Q

(matrix B and vector w are expressed by the formulas (4) and (5)).

The Q and Qα matrices are also theoretical covariance matrices of the estimated parameters (coordinates of the network points), respectively, for the exact and approximate (regularized) solution of the system, assuming that the standard deviation of the observation with a unit weight s0 = 1.

Influence of the regularization parameter on results of Equation (13)

By subtracting the normal equation (4) (least squares) from the normal equation (14) (least squares with Tikhonov regularization), we will estimate the size of the difference between the solutions using the relative error (η):

(18)
B(xα−x)+α⋅xα=0   or   B⋅Δx=−α⋅xα
(19)
Δx=−α⋅B†⋅xα
(20)
η=‖Δx‖E/‖xα‖E≤α‖B†‖S=α⋅‖Q‖S=α⋅λr−1,
where: ‖ . ‖E is the vector Euclidean norm (associated with the spectral norm ‖ . ‖S of the matrix) and λr is the smallest non-zero eigenvalue of the matrix B. Then, λr−1 is the largest eigenvalue of the matrix Q = B†. We can see that the relative error depends significantly on the value of the regularization parameter in relation to the smallest, non-zero eigenvalue. In the literature, e.g. Bell and Roberts (1973), we also find the estimate of η ≤ α(λr + α)−1, which, including the condition α ≪ λr, should not differ significantly from (20). In any case, however, the estimate using the eigenvalue λr is practically not verifiable unless we know the spectral decomposition of the matrix B. Such an operation is not necessary as it is too expensive. The value of λr−1 can be estimated with some excess, assuming an average value of the expected maximum error of the mean coordinate of the network point. Let us denote it as μx. From the trace Tr properties of the matrix, we have:
(21)
λr−1≤∑i=1,…,r(λr−1)=Tr(Q)≤nλr−1μx2.

We will now return to (20) using the last estimate:

(22)
η≤α⋅n⋅μx2.

If the relative error was to be smaller than ηmax, the parameter α should be taken as

(23)
α<ηmax⋅(n⋅μx2)−1=αmax.

For example, let ηmax = 0.01, μx = 0.01 m, n = 1000 (number of unknowns), then α < 0.001.

The parameter α should also be limited from the bottom due to the rounding error in a digital machine (Wilkinson, 1994). For adding α to the diagonal element of the matrix B to be a significant operation, its value should exceed the rounding error level of the diagonal elements of the matrix. Taking into account (23), we will get a selection range (Kadaj, 1979):

(24)
α∈(max{bii}⋅2−t,αmax),
where bii is the diagonal element of the matrix B and t is the number of significant binary digits in the floating-point notation of a real number.

An estimation of the relative error of the Tikhonov regularization method concerns the case of a one-time solution according to formula (14) or, otherwise, formula (13) in relation to formulas (4) and (6). In practice, regularization is used in iterative processes that correct a one-time solution (e.g. Bakushinskii (1992); Bell and Roberts (1973); George (2010); Hanke and Groetsch (1998); Marquardt (1963); Mehsner (2013); Scherzer (1993)). In particular, this applies to non-linear least-squares tasks using the Gauss–Newton iterative procedure, and thus to the adjustment of geodetic networks also.

Regardless of whether the solution (14) is used once or is refined in the iterative process, the covariance matrix of the unknown vector determined by formula (13) is approximate. The assessment of the relative error of this matrix in terms of the relative increase of the spectral norm leads to a restriction on the parameter α analogous to (24), but with a ‘sharper’ upper limit (Kadaj, 1979): αmax=(12)⋅ηmax⋅(n⋅μx2)−1 .

Using Tikhonov's regularization, is it possible detect and locate a configuration defect in a geodetic network? The covariance matrix Qα expressed by formula (16) and the pseudo-inversion B† do not contain information that would allow for the identification and location of the defect. Such information is provided by the Bα−1 matrix. Using the spectral decompositions, traces of these matrices will be compared:

(25)
Tr(B†)=λ1−1+…+λr−1+0+0+…+0=C1+0
(26)
Tr(Qα)=λ1(λ1+α)−2+…+λr⋅(λr+α)−2+0++ 0+…+0=C2+0
(27)
Tr(Bα−1)=(λ1+α)−1+…+(λr+α)−1+α−1++ α−1+ … α−1=C3+u⋅α−1
where u is the network defect. Since α is a small number in relation to λr, the approximate equality holds C1 ≈ C2 ≈ C3 and the dominant big value component u · α−1 appears only in the trace of the matrix Bα−1 . It is located in the diagonal elements of the matrix inverse to Bα, or more precisely in the elements corresponding to the coordinates of undetermined points. The configuration defect is then identified in large values of the mean errors a posteriori of the respective coordinates (the diagonal elements of the Bα−1 matrix are squares of mean errors). The next section focuses particularly on the choice of parameter for this purpose.

Selecting the parameter α for detection and location of a configuration defect

Let the regularization parameter α be interpreted as α=1/μc2 , where μc denotes a priori expected mean coordinate errors for the unknown parameters (point coordinates). In a probabilistic and statistical sense, it would also mean applying the idea of Bayesian estimation. Contrary to the commonly used Gauss–Markov model for observation adjustment, Bayesian estimation additionally uses a priori information (stochastic model) for the unknown parameters. In our case, on the one hand, it will free the task from the defect and singularity of the normal matrix, while on the other hand, it will allow to identify missing observations or those given with erroneous point numbering. In the general model of Bayesian estimation (e.g. Bossler (1972)), instead of the scalar matrix α · I used here, a full CX covariance matrix may appear as additional a priori information about the determined unknowns.

Suppose the matrix of a system of normal equations B(n x n) is singular due to a network defect, rank(B) = r < n. Assume, for example, the coordinates of all points determined in the network μx = μy = μ c = 100 m, that is, with the weight α = 0.0001 do not significantly affect the results of the network adjustment, we make the non-singular system of normal equations (the network becomes determinable, and the normal matrix becomes Bα = B + αI and is significantly positively defined with the order r = n). Let us add that all eigenvalues (including zero values) of matrix B are increased by the added parameter (λi = λi + α; i = 1, 2, . . . , n), and the inverse matrix Bα−1=(B+α⋅I)−1 exists and has positive eigenvalues 1/λi = 1/(λi + α).

In our case, the upper limit results from the assumption that the diagonal elements of the inverse Bα−1 matrix for indeterminate points are significantly different in size in relation to the diagonal elements of fully determinable points. If the network was completely determinable (with zero u-defect), then the diagonal elements of the matrix (assuming the use of standardized observations with the unit mean error value a priori) would be the squares of the coordinate mean errors. For example, for a horizontal geodetic network, the use of the regularization formula means that we add the equations dx = 0, dy = 0 to the existing system (in the previously adopted notations, dx, dy are some components of the vector x) with assumed a priori (artificially) large standard deviations (mean errors), for example, 100 m, and with appropriately small weights (as above) to reveal undetermined points with sufficiently large values of the resulting mean errors. For this purpose, we assume the condition for selecting a regularization parameter (left part as in (24)):

(28)
max{bii}2−1<(k⋅μmax)−2≤α≪μmax−2          ('≪' signifies a strong inequality),
where μmax is the expected maximum (approximately) error of the coordinate of the determinable point in the network, k is the multiple of the assumed mean error for pseudo-observation dx = 0, dy = 0, for example, k = 10, 100, 1000, . . . (see attached numeric example) and t is the number of significant binary digits of a real number.

. Illustrative example of detection and location of configuration defects in a large traverse geodetic network

The following example illustrates the possibility of detecting configuration defects using the Tikhonov regularization. The object of the test is a relatively large geodetic network with a traverse structure (Figures 3 and 4). The network includes 6227 points with 448 reference points, 6298 angular observations, 1487 directional observations at 391 stations and 6829 length measurements.

Figure 3

A symbolic sketch of an extensive geodetic network in a district in Poland. In the Polish classification, it is a third-class network related to the network of classes I and II. The sketch shows a place in the network where configuration defects have been artificially introduced.

https://rgg.edu.pl/f/fulltexts/173249/j_rgg-2021-0004_fig_003_min.jpg
Figure 4

Localized configuration defects (numerical results in Table 1).

https://rgg.edu.pl/f/fulltexts/173249/j_rgg-2021-0004_fig_004_min.jpg

In the professional adjustment of this network, the following standard deviations of observations were assumed: for directional observation: 28 [cc], for length measurements: constant standard deviation 0.016 m and proportional to the length of 0.0002 m/km. The average length of the side is approx. 350 m. In the resulting report, the estimated average standard deviation of the point was 0.038 m.

Artificial configuration defects were introduced in the network by excluding five angles at the vertices of the traverse marked in Figure 4. The classical least-squares adjustment algorithm signals the singularity of the normal matrix and the program is stopped. We know that the defect exists, but we are unable to locate it and quantify it. In practical situations, detecting configuration defects in such a large observation system would not be an easy task.

The Tikhonov regularization procedure becomes helpful here. In accordance with the example provided in section 4.3, we assume the parameter α = 1/10, 000 (according to the assumption: μx = μy = 100 m). Table 1 shows a fragment of the network adjustment report, where it is clearly visible that the resultant point error takes significantly large values just for the points for which there are deficiencies in the centre angles. These values, however, are smaller than 100 m because they occur in the vicinity of 5 points and the observations (angles and lengths) measured at the neighbouring points ‘work’ in a sense for the points with configuration defects.

Table 1

Fragment of the result set (coordinates and accuracy parameters) of the network adjustment containing configuration defects using Tikhonov regularization with the parameter = 1/10,000 (assuming a priori standard deviations of the coordinates μx = μy = 100 m)

Point Nox [m]y [m]s [m]a [m]b [m]θ[m]
113510515734128.73137593963.79170.02290.01650.0159114.89
113510525733836.11227593815.79400.03370.02650.0207125.33
113510315734659.96227593795.21900.03380.02890.017547.35
113510305734794.18187593605.84200.05550.05000.024143.80
113510295734961.87977593605.26040.06910.06330.027757.73
113510285735150.90547593600.77690.08830.08200.032969.19
113510275735344.55307593600.107199.999399.99930.050699.78
113510265735593.18987593599.205790.012390.01230.055899.77
113510255735845.27087593593.379143.595243.59520.0460101.09
113510245736163.62107593598.83950.08480.07920.030274.12
113510235736508.71767593598.48210.05690.05120.024955.53
113510225736792.56647593337.16970.02760.02140.017461.38
113511125737672.61617592856.10530.02730.02070.017856.84
113511025737876.89707592817.78280.03040.02320.019649.43
113511015738021.95207592509.61900.03570.02760.022649.04

[i] Notations: a, b, θ – ellipse parameters standard deviations; a, b (b ≤ a) – the semi-axes, estimated standard deviation of the point: s = s=(μx2+μy2)1/2 ; θ – directional angle of the semi-axis a; x, y – Cartesian Gauss-Krüger coordinates in PL-1992 system.

Table 2 shows the analogous results of the actual network adjustment (without defects). Comparing the errors of points that do not have a direct observational relationship with the ‘defect’ points in Tables 1 and 2, we find that the Tikhonov regularization practically does not change the accuracy parameters of the network, except for the defective places.

Table 2

Fragment of the result set (corresponding to the data in Table 1) without defects and without Tikhonov regularization

Point Nox [m]y [m]s [m]a [m]b [m]θ[m]
113510515734128.73077593963.79360.02090.01530.014320.80
113510525733836.11137593815.79650.03100.02370.0199126.34
113510315734659.96057593795.21310.02630.02080.016253.54
113510305734794.17607593605.82990.03890.03220.021848.54
113510295734961.87667593605.24000.04450.03790.023458.27
113510285735150.90487593600.74380.05000.04310.025466.20
113510275735344.55537593600.10440.05380.04670.026671.41
113510265735593.19557593599.38100.05580.04880.027074.54
113510255735845.26707593593.03630.05450.04770.026374.50
113510245736163.61397593598.84980.04870.04210.024569.68
113510235736508.71377593598.48770.03920.03240.022157.19
113510225736792.56387593337.17180.02300.01660.016087.27
113511125737672.61617592856.10530.02620.01990.017156.84
113511025737876.89707592817.78280.02920.02230.018949.43
113511015738021.95207592509.61900.03430.02650.021749.04

[i] Designations as in Table 1.

The described method of automatic detection and localization of configuration defects, based on the Tikhonov regularization idea, was implemented in 1993 in the GEONET geodetic calculation system (www.geonet.net.pl). The test calculations presented in this paper were performed by means of the programs of this system.

. Conclusions

The article presents the methodology of automatic detection of configuration defects in geodetic networks. This issue is important, especially in large networks, at the stage of initial measurement processing, where detection of possible data defects is usually a very labour-intensive task. Computer programs for adjusting geodetic networks ‘recognise’ the geometric structure of the network based on the identifiers (numbers) of points in the observation plans, in the list of approximate coordinates and in the list of tie point coordinates. In case of errors in point identifiers, the network structure is incorrectly defined. In particular, errors in identifiers may result in situations similar to the case of missing data.

Both the theoretical basis and numerous practical tests performed with the GEONET system programs show that the described method based on the Tikhonov regularization idea allows to automatically detect possible configuration defects in the geodetic network, interpreted by the computer program as gaps in data sets. In ‘raw’ observational data, as practice teaches, this will often be the result of errors in point identifiers.

When choosing a regularization parameter, we follow its pseudo-probabilistic (pseudo-Bayesian) interpretation. In this interpretation, the regularization parameter means very small value of the weight of the coordinate of the determined point (e.g. 0.0001), which corresponds to ‘excessively’ large value of the standard deviation of the coordinate of the determined point (e.g. 100 m). Such an assumption does not significantly affect the results of aligning the part of the network that is free from configuration defects, but allows to identify defect sites by showing large values of diagonal elements of the inverse of the regularized normal matrix.

The problem of configuration defects may occur in any real observation system, including hybrid networks (satellite-terrestrial, integrated). The given example of the classical network is only an illustration of the problem and its solution using the Tikhonov regularization.

Regardless of the possibility of applying the Tikhonov regularization, there are other methods supporting the identification of configuration defects in the geodetic network. One of them, implemented in the GEONET system, consists in the analysis of the topological structure of the network in terms of counting the number of independent observations that mark each point of the network in relation to other points of this network. However, it is easy to give examples of when such an approach will not work. It may be, for example, a network fragment not related to the rest of the network and not having tie points, but meeting the condition of correct (relative) determinability in relation to other points of this fragment.