On the Helmert-blocking technique: its acceleration by block Choleski decomposition and formulae to insert observations into an adjusted network

The Helmert-blocking technique is a common approach to adjust large geodetic networks like Europeans and Brazilians. The technique is based upon a division of the network into partial networks called blocks. This way, the global network adjustment can be done by manipulating these blocks. Here we show alternatives to solve the block system that arises from the application of the technique. We show an alternative that optimizes its implementation as the elapsed processing time is decreased by about 33%. We also show that to insert observations into an adjusted network it is not necessary to readjust the whole network. We show the formulae to insert new observations into an adjusted network that are more efficient than simply readjusting the whole new network.

2. Introduction

2.1 Review of the Helmert-blocking technique

This work deals with the Helmert-blocking (HB) technique, an approach designed for the adjustment of large geodetic networks, like Europeans and Brazilians. Developed in the nineteenth century by the geodesist Helmert [1], some recent examples of the use of the HB technique are as follows: to adjust the North American Datum of 1983 (NAD83) and to integrate the Canadian geodetic network to NAD83 [2], where Canada developed for this purpose, in the 1980s, the computer system GHOST (Geodetic Adjustment using Helmert Blocking of Space and Terrestrial Data) [3]. This software tool had been used in Brazil since the 1990s by the Institute for Geography and Statistics. The transition to GHOST in Brazil occurred when computer systems based on classical methods of adjustment started to become unfeasible [4].

GHOST was also used in Brazil to adjust the horizontal network of the Geodetic Brazilian network (RGB) to SIRGAS2000, and it is also being employed in the adjustment of the levelling network of RGB to SIRGAS2000 [4–6]. Furthermore, the HB technique has also been recently used in other parts of the world. The International Association of Geodesy Reference Frame Sub-Commission for Europe used the HB technique in the adjustment of a network of continuous operating GPS stations (EPN) [7]. Nocquet et al. [8] also used the HB technique on the generation of a plate kinematics model for Nubia–Somalia region.

The HB technique consists in subdividing the global geodetic network into partial networks referred to as blocks. This way, the adjustment of the global network can be done through the manipulation of these smaller blocks. Regarding the division of the global geodetic network into blocks, several criteria may be used [1,2], but this issue is out of the scope of this article.

2.2 Outline of the paper

HB was and is still currently used on the adjustment of geodetic networks. Recent uses of HB include other kinds of data as well (e.g. GPS and plate kinematics). Thus, considering the small number of references to the mathematical formulae of this technique and considering that the deductions are taken for granted or simplified in these instances, here we formally deduce the formulae of the technique in detail but concisely. In contrast to Wolf [1], here we derive HB general formulae from least squares (LS) formula based on properties of block matrices (§3). Then we derive formulae to insert observations in a geodetic network adjusted in advance based on this a priori adjusted solution (§4). Of note, the derivation in §3 serves as a support for the proposition of the formulae for the insertion of observations.

In previous work, we, together with other colleagues, developed a framework to deduce and implement the technique under Matlab [9]. The proposed framework led us to derive alternatives to implement the HB technique and to conceive softwares for the adjustment of geodetic networks based both on the classical HB approach [1] and on our alternatives which are: HB by block Choleski decomposition (HBC) and HB by block Gaussian elimination (HBG). However, the tests developed and presented by Lema et al. [9] were very restricted due to the nature of the generation of random matrices and especially the small size of the junction unknowns vector (hereafter referred to as junction vector). Here we take the simulations to another level and we present broader tests that validate the deductions and implementation and show once more that the proposed modification in the technique optimizes it and that this optimization is enhanced when the junction vector is very large (§5), which is the case for large geodetic networks.

Previously, in [9], we verified that the proposed changes optimized the implementation but without quantifying this improvement (due to the restricted data generated in this previous simulation). So, here, we extensively tested the technique under several geodetic networks of varying configurations.

3. Deduction of the Helmert-blocking's formulae

While handling observations of large geodetic networks, it is inevitable to face thousands or even millions of equations involving thousands of variables. For instance, Schwarz & Wade [2] handled 1 785 772 observations and 928 735 unknowns in the adjustment of NAD83, in which Wolf's [1] formulation of the HB technique was adopted. As a consequence, it is inviable to employ the classic LS method to the adjustment of the geodetic data of the whole network due to the great dimension of the matrices involved in it.

On the other hand, in many instances, these equations are previously organized in observation sets called blocks (or can be organized as such), giving particular characteristics to this system of equations [1,2]. Such block configuration is more adequate to exploit the sparseness of the network's matrix. And the HB technique manipulates these blocks to give an adjusted solution to the global network while taking advantage of this sparseness of the network's matrix [2].

We will include here the latest derivation of HB's formulae which goes after [9]—we include these derivations here because they will be important for §4 and because [9] is in Portuguese only.

The variables of each block may be classified according to their participation in neighbouring blocks.

• — Junction unknowns. These are the variables corresponding to stations belonging to more than one block of the geodetic network.

• — Interior unknowns. These are the variables belonging to only a single block of the geodetic network.

Assume the dataset of a given geodetic network is divided in n blocks, each having at least one benchmark, all benchmarks referring to the same reference system. Assuming that the network as a whole has at least one junction unknown (the adjustment in the absence of junction unknowns is as outlined in §4.5 which means that each block can be adjusted individually as they do not have any influence on the adjusted solution of the other blocks) and that all the observations have the same quality (the solution considering individual weights for the observations is analogous and will be outlined later, in §3.2), we have, for a given block i of this network [1]:

3.1

where x_i=vector corresponding to the interior unknowns of block i; y=vector corresponding to the junction unknowns; A_i=matrix of the interior unknowns of block i; B_i=matrix of the junction unknowns of block i; l_i=vector of the observations corresponding to block i; v_i=vector of the residuals of block i's observation.

Figure 1 shows the defined terms, some of which are involved in equation (3.1). In this figure can be observed a geodetic network partitioned in three blocks, its junction stations, its interior stations and the measured height differences. For further examples of the block's division and especially the different levels of junction unknowns (this concerns its division into subvectors, see §3.1) we suggest the reader to refer to [2,10].

Rewriting equation (3.1), an alternative representation through block matrices of the global geodetic network Ax−l=v follows:

3.2

Applying LS method to the global geodetic network provides an adjusted solution x given by A^tAx=A^tl [11]. Equation (3.2) will give rise to the corresponding block representations of A^tA and A^tl.

3.3

The HB technique corresponds to the LS method applied not to Ax−l=v but to its block representation given by equation (3.2). Thus, from equations (3.2) and (3.3) the block representations for A^tA and A^tl follow:

3.4

and

3.5

The solution by the HB technique is achieved by solving the matrix equation A^tAx=A^tl, where x=(x₁x₂⋯x_ny)^T. Hence, it is identical to the solution of LS when applied directly to the global geodetic network. Yet, further algebraic manipulations will yield the final HB expression.

Equations (3.4) and (3.5) yield

3.6

Making the following substitutions:

3.7

equation (3.6) is expressed in a shorter expression:

3.8

Isolating x_i in the first line of equation (3.8) and substituting it on the second line of equation (3.8), the HB solution of the geodetic network arises

3.9

The network's adjusted solution is obtained by calculating primarily y from the first line of equation (3.9). Then, each x_i is determined by substituting this computed y-value onto the second line of equation (3.9). Usually, the Choleski method is to be applied to solve the system in the junction vector y [1]. After all, is a symmetric matrix and weighted in its diagonal, a feature of symmetric positive definite matrices [11]. This precludes the need for pivoting and allows a fast Choleski algorithm to solve the system.

3.1 Junction unknowns represented by subvectors

In most cases (as for large geodetic networks), the great dimension of the junction vector, together with particular properties of the geodetic network in consideration, makes it expedient to divide it in subvectors [1]. For instance, NAD83 adjustment had a sizable junction vector for each of its blocks, generally comparable to the blocks' interior vector [2]. In such a case, the junction vector, y, grows too large and this becomes a hindrance to its direct manipulation. Without alternative approaches, the junction vector would extend to hundreds of thousands of unknowns, and thus its corresponding system would have equations to this same amount and the initial division of the network into blocks would become meaningless. Meanwhile, despite the network's block division, the blocks B_i are still sparse, as actually, very few, or no junction unknown, belongs to all blocks. So, here we will show how to achieve an adjusted solution starting from a division of the junction vector into subvectors:

3.10

This for instance, yields the following substitutions:

3.11

3.12

3.13

3.14

The computation of the junction vector y is still done by equation (3.9). But, in this case the system to be solved is a block system given by equation (3.15), according to the representations adopted in equations (3.10)–(3.14):

3.15

where

3.16

An additional change of variables yields a concise expression

3.17

The block system corresponding to the subvectors of the junction vector y is thus as follows:

3.18

To solve the block system of equation (3.18), one can turn to the literature, and use standardized approaches. The one presented by Wolf [1] resembles a Gauss–Jordan elimination adapted for block matrices; however, it is not formally presented. In fact, the term block system for this equation is not mentioned nor is the technique used (Gauss–Jordan elimination). Therefore, we consider the quality of Wolf's [1] approach inferior to a Gaussian elimination, i.e. a block Gaussian elimination [11], as a Gaussian elimination can be from 1.5 to 3 times more efficient than Gauss–Jordan elimination [12]. Furthermore, here we have formally divided the junction vector into subvectors, making a thorough deduction, as Wolf's [1] is concise.

As presented in the previous subsection, the matrix that makes up the system for the junction unknowns is symmetric and generally positive definite as its elements which are greater in absolute value are spread along its diagonal. So here, to take advantage of this characteristic, we also propose to use the Choleski block-triangular decomposition to solve this block system, since block Choleski's algorithm demands n³/3 floating point operations (flops), instead of 2n³/3 demanded by a block Gaussian algorithm [11]. In §5, we present a numerical experiment to attest the validity of such HBC approach. We evaluate its efficiency and accuracy, showing that as expected, it is more efficient than HBG while being just as accurate.

For details on how to implement HBC and HBG refer to appendix A. Another option yet unexplored in the HB context is the use of iterative methods, e.g. a block implementation of the conjugate gradient method might also be of help in future investigations.

3.2 Introduction of weights in the observations

In most cases, each observation has a weight corresponding to the quality of its measurement—e.g. the inverse of its variance (note that usually the covariance between observations is null). Let P_i be the weight matrix corresponding to the set of block i's observations, which has as observations vector l_i, equation (3.1). Therefore, the geodetic network as a whole has the following weight matrix:

3.19

Taking P as the weight matrix for the geodetic network, the adjusted solution xby LS method becomes the solution to A^tPAx=A^tPl. A way to find the adjusted solution x by the HB technique is to consider the equivalent system A^t(PA)x=A^t(Pl), where

3.20

and

3.21

Then, equations (3.20) and (3.21) may be achieved from equation (3.2) substituting each A_i for P_iA_i, each B_i for P_iB_i and each l_i for P_il_i. Similarly, to obtain the adjusted solution to the global geodetic network by the HB technique considering weighted observations, it suffices to make the corresponding changes of variables in the deduction already done. This way, the adjusted solution by the HB technique is also given by equations (3.7), (3.12)–(3.14), where the following changes of variables are to be made:

3.22

3.23

3.24

3.25

4. Insertion of observations in an adjusted geodetic network

Now, given a network whose adjusted solution was already computed, the problem is to insert new blocks into it without readjusting the entire network. Based on the pre-determined solution, to compute the adjusted solution of the unknowns of the inserted blocks and the corrections to be added to the original pre-determined solution.

Concerning the blocks to be inserted, given one of them, there are four possibilities for the unknowns therein: first case, they comprise unknowns of the original network only; second case, they have their own interior unknowns but junction unknowns in common with the original network; third case, they have their own junction unknowns but interior unknowns in common with the original network; and fourth case, they have no unknown in common with the original network.

4.1 A general routine to insert new blocks in an adjusted network

A general routine to adjust the network based on the adjusted solution of the original one without reprocessing the entire network is as follows:

• 1. Organize the blocks to be inserted into four sets, a set for first case blocks, a set for second case blocks, a set for third case blocks and a set for fourth case blocks.

• 2. Compute corrections inserting only the first case blocks based on equations (4.2) and (4.3).

• 3. Insert the second case blocks and compute the adjusted solution for its unknowns and the second corrections based on equations (4.4) and (4.5) (second case blocks are inserted in the network comprising the original one plus first case blocks).

• 4. Insert the third case blocks and compute the adjusted solution for its unknowns and the third corrections based on equations (4.6) and (4.7) (third case blocks are inserted in the network comprising the original one plus first and second case blocks).

• 5. Compute the adjusted solution for the unknowns of the blocks enclosed in the fourth case separately using standard HB or HBC/HBG formulae of the previous section (no corrections to be added in this case, see §4.5).

In the following subsections, the formulae to compute the adjusted solution for the new unknowns and the corrections of the original ones for each case will be deduced. If the network does not have junction unknowns the possibilities are reduced to two: first case and fourth case. The corresponding solution is obtained by setting B and D to zero in these cases. Moreover steps 3 and 4 of the general routine introduced have no place in this situation but only steps 1, 2 and 5.

The original network has n blocks and its adjusted interior and junction vectors are x₀ and y₀, respectively. To this network, N new blocks are to be inserted in each case and the corrections to the interior and to the junction vector are dxand dy, respectively. For first and third case blocks N≤n as they are linked with the interior unknowns of the original network; if N<n, the blocks of unmatched unknowns are to be filled with corresponding null matrices instead of C_i and D_i. For second and fourth case blocks N is boundless.

The original network is represented by a matrix A for the interior unknowns and a matrix B for the junction unknowns and in like manner, the blocks to be inserted are represented by a matrix of interior unknowns C and a matrix of junction unknowns D. The matrices A, B, C and D relate to the matrices for each block as follows: