Extending geodetic networks for geo-monitoring by supervised point cloud matching

2.1 Network-based monitoring

Using a geodetic network consisting of measurements to selected and signalised single points is the conventional approach in geodetic monitoring [10]. Based on the observation of identical points in each epoch, a rigorous deformation analysis including statistical significance tests can be performed. Hereby, the calculated deformation values are tested to determine whether they are statistically significant or not relative to the accuracy of the point determination and based on the selected confidence level [2]. The network consists of a set of discrete (i) control points, which are assumed to be stable, and (ii) object points, located in the moving area characterising the deformation process. Creating a network structure of geodetic measurements between the points leads to redundancy and ensures a realistic uncertainty estimation for the determined point coordinates. In each epoch, the identical point field is observed and evaluated in a free network adjustment. If the variances of the two measurement epochs are comparable, point deformations are calculated and localised according to the congruency model. The congruent control points define the datum and thus the coordinate system for the deformation analysis, whereas some control points may also be eliminated and added to the set of object points as they are detected unstable. As a result, the 3D displacement vector for each object point (including eliminated control points) and its proof of significance is determined.

However, when performing network-based monitoring, the deformation process is generalised by a number of object points. Prior assumptions about the deformation process have to be made in order to select representative points. The limited spatial resolution also limits the detailed modelling and understanding of the deformation process.

2.2 Point-cloud-based monitoring

Instead of discretising the deformation area by a number of carefully selected and signalised points, the surface geometry of the whole area can be represented by a dense point cloud with a high spatial resolution captured by TLS. The main challenge in comparing two point clouds is overcoming the problem of missing identical points across the two epochs. This is because the scanning raster is not reproducible in the different epochs, especially not if the object surface has changed between the epochs due to deformations and thus, the measurement points will not fall strictly on the same locations. Nevertheless, various methods have been developed for deriving deformations from two point clouds. According to ref. [7], the different approaches can be grouped into three categories: point-cloud-based, point-based and parameter-based strategies. Not all methods should be addressed here, but some are mentioned in the following.

Point-cloud-based approaches usually use the complete point clouds of both epochs to determine changes between them. The result is often a colour-coded heat map which depicts the magnitude of geometrical differences between the two point clouds. The comparison is either based on the point cloud itself, like in the M3C2 algorithm [12], 13] or on best-fitting surfaces from triangulation [14], 15] or free-form parametrisation [16].

Parameter-based methods approximate the point cloud or parts of it using analytical surfaces like geometrical primitives [17]–20]. Deformations are then identified analysing changes in the estimated parameters.

Point-based strategies reduce both point clouds to a number of identical points, which then can be compared directly across the epochs. Comparable points are created by reducing and interpolating the irregular point clouds to a specific grid [21], by estimating and intersecting planes [22] or by extracting feature points based on machine learning [23] or on additional RGB images [24]. Furthermore, features can be extracted from the 2D representation of a digital terrain model (DEM) created from the laser scans for each epoch [25], 26].

In geo-monitoring, mostly point-cloud-based and point-based strategies are used as natural objects are difficult to describe by geometric primitives. However, only point-based strategies result in actual 3D deformation vectors (magnitude and displacement direction). Other approaches have to accept a loss in dimension because, in most cases, only one-dimensional offsets between surfaces, hence in the direction of the local surface normal, are calculated. Moreover, it should be mentioned that a precondition for most of the methods is that the two point clouds are already located in the same coordinate system. Therefore, an accurate point cloud registration over two epochs is another challenge in point-cloud-based monitoring, because registration errors later directly add to the actual displacement values [27].

2.3 Patch-based monitoring

A special case of point-cloud-based monitoring involves subdividing the point cloud into tiles or even using only selected patches [28]–30]. In this case, a rigid body movement is calculated individually for each point cloud patch, for example, using the iterative closest point (ICP) algorithm [31]. This process is known as ICP matching. Although the number of object points and thus the spatial resolution is reduced in comparison to other point-cloud-based strategies, ICP matching gives an actual 3D vector for each patch [32] and the patches can even be seen as virtual target points [22], 33].

Using patches makes the deformation analysis more resistant to measurement noise depending on the number of scan points per patch. However, even for point-based strategies where identical points are created between epochs, the proof of statistical significance remains challenging due to the unknown accuracy information of the derived points. Given the numerous influences on measurements and ICP matching, a precise estimation of the stochastic model is nearly impossible. Therefore, we propose a novel monitoring strategy combining network-based monitoring and patch-based monitoring.

3.1 Virtual targets from ICP matching

The ICP algorithm determines the 3D transformation parameters between two point clouds of the same patch P. In case the two point clouds were captured in two different measurement epochs k, the parameter set can be used to transform any point from the first epoch k = 0 to a follow-up epoch k > 0. Based on this concept, for each patch P _i a virtual target P _i is created in all epochs k (Figure 2a):

1. Define virtual point P _i,0 in zero epoch. We use the patch’s centroid as the virtual point, but any other point can also be chosen.

2. Determine transformation T _i,k between zero epoch and epoch k using ICP

   (1) T ( t | R ) i , k = I C P ( P i , 0 → P i , k )

3. Transform virtual point P _i,0 into epoch k → updated virtual point P _i,k

   (2) P i , k = T ( t | R ) i , k P i , 0

Figure 2:

Procedure for creating corresponding points using ICP matching (a) across two epochs and updating the virtual point, and (b) between two stations and integrating them into the geodetic network consisting of control points (triangles) and object points (circles).

Given that the ICP algorithm determines the rigid body transformation between two point clouds, our workflow assumes there are no deformations within one patch. Typically, we use manually selected point cloud patches with a size of approximately 30–50 cm, depending on the surface structure.

The ICP matching cannot only be performed across two epochs but also between two stations (Figure 2b). If the same patch is scanned from two different stations within the same epoch, ICP matching can be used to create redundant observations. The procedure remains the same: the virtual point of the reference patch is transformed into the point cloud of the second station using the transformation parameters determined by ICP. The resulting redundant single-point observations of the patches improve the overall network geometry and, thus, the accuracy of the points (see Section 3.3).

However, despite different viewing angles on the patch, the ICP still has to give reliable and correct results. Regarding the network geometry, it is desirable to have very different viewing angles on the object, but for ICP it is better to have a very similar viewing angle. This may cause that not for all patches redundant observations can be created.

3.2 Supervised selection of structures

The quality of the ICP result depends on different factors [35], especially the measurement configuration, patch geometry, measurement noise and scan resolution. Consequently, for each patch has to be checked if – in case of an arbitrary deformation – the matching is (i) unambiguous and (ii) reproducible. Both conditions are primarily influenced by the structure’s geometry:

1. A lack of geometric variation within the structure leads to the matching being inaccurate in certain directions, for example in the case of a (quasi) planar structure (Figure 3b).

2. Repeating geometries within the structure lead to several local minima to which the ICP converges, depending on the occurring displacement (Figure 3c).

Figure 3:

Example of a (a) suitable, (b) ambiguous and (c) repeating geometry for ICP.

To achieve unambiguous and thus reproducible matching the selected patch has to have a certain geometric variation (Figure 3a). Only structures that allow reproducible matching in all three dimensions should be used. In other words, the ICP should always converge to the exact same local minimum, independent of a possible rigid-body movement between two epochs. Verification can be done either analytically, by analysing the structure’s geometry, or statistically in a simulation.

Using the second approach, besides the geometry also other factors contribute to the result, for example the scan resolution [36]. Hence, we perform a MC simulation for each patch pair P _i and P _j that is matched.

Within the MC simulation, the ICP matching is repeated n times for the same patch pair. After an initial alignment of the point cloud pair, in each iteration m the patch P _j is moved by a certain displacement d _m relative to the reference patch P _i .

The maximum displacement is limited by a bounding box (Figure 4). Within this bounding box, n uniformly distributed displacements are randomly created.

Figure 4:

Principle of the MC simulation in 2D. The second point cloud is successively moved to a number of uniformly distributed starting points within a predefined box and the ICP is performed.

Each moved patch P j , m ′ is then matched on the reference patch P _i using the ICP algorithm.

As already emphasised, the ICP should always convert to the same minimum and thus the initial registration result. The deviation between nominal displacement d _m and translation t _m given by the ICP matching is the registration error e _m

The variation of e over all iterations m can be taken as a measure for the matching precision u _i,j .

For suitable patches u _i,j is expected to be small. However, depending on the geometry, a certain variation in the registration is expected. Thus, a suitable threshold u _max for the registration precision within the MC simulation has to be defined (see Section 4.1), whereby all elements of u _i,j must be below the threshold value u _max. Selected structures that do not fulfil this condition can be eliminated.

The point clouds may be captured from the same station in two different epochs, from two different stations within the same epoch, or even from two different stations in two different epochs. The MC simulation evaluates the actual point clouds later used in the ICP matching. Consequently, the analysis evaluates not only the patch geometry but also the actual measurement configuration, measurement noise, and scan resolution. However, since these latter parameters remain constant throughout the MC simulation, all point cloud pairs P _i , P _j of the same structure must be analysed individually. The box size should be chosen according to the expected displacement to ensure the patch is unambiguous within the border of a maximal expected movement.

An example of supervised structure selection is shown in Figure 5. In this case, a boulder was scanned and the point cloud was divided into rectangular patches with a size of approx. 30 cm by 30 cm. After deformation, the boulder was scanned again and the exact same grid was placed over the point cloud, albeit with a slight shift in grid position relative to the object due to the movement. Nevertheless, the patches can be considered homologous, given that the movement is significantly smaller than the patch size. ICP matching gives a deformation vector for each patch. As the monitoring object is one solid rock, it can be expected to move homogeneously, resulting in similar deformation vectors across all patches. However, Figure 5b reveals several outliers in the results. These can be attributed to a lack of geometric variation in certain patches. The MC simulation can eliminate all unsuitable patches (Figure 5c). In this case, a threshold value of u _max = 1 mm was employed. Ultimately, the vectors of the remaining patches (highlighted in green) show a homologous movement pattern.

Figure 5:

Example for supervised structure selection on a boulder: (a) laser scan of the rock captured in two epochs, (b) random tiling of the block and resulting deformation vectors using ICP matching and (c) filtered patches after the MC simulation.

In contrast to comparable methods, such as feature-based approaches, which often use histograms to filter out false results [26], our proposed method allows for the examination of expected result quality prior to the actual calculation.

3.3 Network adjustment and rigorous deformation analysis

After unsuitable patches have been rejected, still no realistic uncertainty information for the remaining virtual points is available. In the context of point-based monitoring, geodetic network measurements and adjustments allow the most consistent error propagation and, thus, a comprehensive significance calculation. Consequently, the virtual points are converted into tacheometric pseudo observations l _i,k , each consisting of horizontal angle Hz, vertical angle V and slope distance S

where Δx, Δy and Δz are the coordinate differences between the virtual point and measurement reference point of the instrument. Finally, all pseudo observations of the same epoch k can be combined with other single-point measurements from the total station or GNSS and evaluated in a free network adjustment. At the end not only adjusted coordinates x ̂ for all observation points result, but also the corresponding covariance matrix Σ x ̂ x ̂ containing precision information for each point.

Besides the observation accuracy, the uncertainty of the adjusted points is also influenced by the network design. Generally, a higher redundancy improves the reliability and controllability of the network, thereby leading to higher quality. A measure for the controllability of each observation and thus the probability of identifying possible gross errors are the partial redundancies r _i . If within the network adjustment a global test is performed to ensure that the variance factor a priori equals the variance factor a posteriori σ 0 2 = s 0 2 , the partial redundancies r _i can be derived from the diagonal elements of the redundancy matrix according to

with the covariance matrix of the corrections Σ _vv and the covariance matrix of the observations Σ _ll . A partial redundancy of r _i  = 0 indicates an uncontrolled observation, while a value of r _i  = 1 reflects full controllability. By comparing the partial redundancies, the benefit of additionally included pseudo observations on the network design can be examined.

Point deformations between two epochs k − 1 and k are then typically calculated based on the congruency model [2]. This requires the presence of corresponding points in both epochs in order to calculate coordinate differences d and the corresponding covariance matrix of the differences Σ _dd .

The global congruency test is then checking if deformation are statistically significant

with the testing value T _d , the threshold F of the Fisher distribution, the redundancy r of the network and the significance level α.

4.1 Laboratory setup

A small monitoring network was established in the laboratory (Figure 6). It consists of four stable control points and ten moved object points. All object points are installed on a board that can be moved horizontally and thus create displacements. A combination of signalised and non-signalised target points was used. Standard round prisms, mini prisms and reflective tapes are common target signalisations and give reference values for the displacement. Furthermore different structures can be found on the board, suitable for the proposed ICP matching approach. Some elements (white corners) have a more flattering geometry for ICP than others (stones). The scanning of all structures gives a total of six scan patches for ICP matching. Patch 316 is intentionally chosen to be ambiguous for ICP. The points 101–104 are stable and used as datum points in the deformation analysis.

Figure 6:

The setup in the laboratory consists of a number of reflectors and scannable objects on a board. Deformations are induced by moving the whole plate.

In all epochs, all targets and patches are measured/scanned from two stations (201 and 202). The point clouds have a resolution of 3 mm on the object. After measuring the zero epoch, the board is moved four times and each time a follow-up epoch is measured. The first two movements are carried out in x-direction, perpendicular to the line of sight of the two instrument stations. The board is moved by approx. 1 cm each time. The second two movements then are in y-direction, in line of sight. Again, the magnitude of the displacement is approx. 1 cm each.

However, the results of all four epochs showed a similar accuracy, independent of which direction the displacement occurred. Thus, in the following, we only show the result of the deformation analysis between the zero epoch (EP00) and the last follow-up epoch (EP22), in which we have movements in both directions.

First, the MC simulation is performed for each patch. Especially for the rock faces it is not clear if all of them are suited for ICP matching. As the maximal movement in x and y is about 2 cm each, the used box size for the simulation is 5 cm to cover the full range of expected movement within the simulation. Figure 7 shows the results of the simulation for a corner (311), a rock face (313) and an edge (316). As expected, the edge geometry (316) is ambiguous for the ICP along the edge and the patch has to be eliminated. The resulting precision of the corner geometry, on the other hand, is very good. The rock face, however, shows some variation in the matching result, but it is smaller than the nominal observation accuracy of the instrument. All other patches deliver a very high precision as well, similar to patch 311.

Figure 7:

Result of the MC simulation for the patches 311 (left), 313 (center) and 316 (right).

In the following we assume that with a filtering threshold similar to the measurement accuracy of the instrument, a high network accuracy is secured. Meaning, that no observations with an accuracy way worse than the existing total station observations are included. Thus, a threshold

is selected. Consequently, only patch 316 has to be eliminated in the present case.

Next step is the free 3D network adjustment [37]. To investigate how the network quality can be improved by the additional targets, the partial redundancies are evaluated. However, the network geometry in the laboratory setup is already very good. Therefore, additional observations do not bring any decisive improvement. In many applications, however, the network design is limited due to the measurement environment and there is potential for improvement. For this reason, the network design is also degraded in the laboratory to show that, in this case, the virtual observations lead to an improvement. Only two of the four datum points (101, 103) and only one of the four reflectors (303) on the board are used. Figure 8 shows the resulting network configuration in the zero epoch.

Figure 8:

Partial redundancies of the horizontal angle (Hz) observations for a low number of object points (left) and a high number of object points (right).

Afterwards, the network is extended by the observations to the five remaining patches. Figure 8 depicts the partial redundancies of the horizontal angle observations for both configurations. In the extended configuration, not only do the additional included observations have a high partial redundancy, but the partial redundancy of the already existing observations is also improved.

The improvement can also be seen in Table 1 where for each observation type the average partial redundancy is listed. The biggest improvement is made in the distance observations, while the other two groups also show a large improvement. As the measurement conditions in the laboratory are ideal, the point accuracies are with about 0.1 mm standard deviation already very high for the configuration without patches. Therefore, no improvement could be seen in the point accuracies itself.

Finally, the deformation vectors between the two epochs are calculated for the full network configuration and depicted in Figure 9. All vectors show the same magnitude and direction of movement. If not known, it cannot be distinguished between a marked point or a virtual target from a laser scan. The accuracy of the virtual points is similar to the accuracy of the marked points and again it is in the range of 0.1 mm. The point with the worst accuracy (0.2 mm) is point 304, marked by a reflective tape.

Figure 9:

Resulting deformation vectors for the reflectors (blue) and remaining patches (orange). Total station measurements are shown as blue dashed lines.

This is a first proof that the proposed method works. It also confirms that the chosen threshold for filtering unsuited patches is appropriate as the remaining patches delivered correct results. As already addressed in Section 3.2, the quality of the ICP result is not only affected by the patch geometry but also by the scan resolution.

To create a worse scenario, in each epoch all patches are also captured with a resolution of 6 mm. In this case, the resulting ICP precision for patch 313 is above the threshold and also other patches have to be eliminated. If not, the network accuracy gets worse and the resulting deformation vectors for those virtual points do not match the nominal movement given by the reflectors. This confirms as well that the threshold of 1 mm is well suited.

4.2 Study area Mt. Hochvogel

Mt. Hochvogel is the highest peak in the Allgäuer Alps, located on the German-Austrian border (Figure 10). It is exposed to strong erosion and characterised by brittle, tectonically stressed and fissured main dolomite, making it prone to rockfalls and other natural hazards. Over the past decades, an enormous cleft has formed right on the summit, splitting the entire mountain top into two parts. The southern part threatens to break off completely and fall down on the Austrian side, with an estimated volume of up to 260,000 m³. By now, the main crack extends over a length of about 35 m in SW-NE direction and is already gaping more than 5 m at its widest point, while several smaller crevices are showing high opening rates as well, making the modelling of the failure process more complicated. The increasing number of extreme weather events, caused by climate change, further accelerates the failure process and increases the likelihood of rockfalls.

Figure 10:

The summit of Mt. Hochvogel (2,592 m a.s.l.), located on the German-Austrian border, with its enormous cleft. The southern part (left part in the image) is moving several mm per month and threatens to break off completely.

To better understand the present failure process and to implement an early warning system, several measurement systems were set up on the mountain including a geodetic monitoring network based on total station and GNSS measurements [38]. Additionally, terrestrial laser scans of the main crack are captured using the scanning total station.

In this work, our emphasis is on the total station network and on the examination of the network improvement by extending it by virtual targets derived from laser scans. The total station network comprises marked points located on both the stable and moving parts of the summit (see Figure 11). However, not all marked points are visible from both stations, thereby not increasing the redundancy of the network. Figure 11 illustrates only those marked points that are measured from both stations.

Figure 11:

Measurement setup on Mt. Hochvogel and marked points (blue circles) that were measured from both stations (S1, S2).

The first laser scans were captured in September 2018, followed by a second epoch in July 2019. In both epochs, the main crack was scanned from two stations with a large overlap in the scans (Figure 11). Using a scanning total station the scan resolution, however, is limited. The average scan resolution in both epochs is about 1 cm. As discussed in the previous section, a lower scan resolution makes it harder to find suitable patches.

The goal was to find patches that were scanned from both stations in order to generate redundant observations later. For that, a patch that is selected in the scans from both stations is immediately evaluated with the MC simulation. Again, a threshold of 1 mm is used, as the movement of the southern side of the summit is small and thus the accuracy requirement is high. However, due to the variation of the viewing angles on the object by up to 90° between the two stations (see Figure 11) and the limited area of overlap, ICP matching between the stations fails for many patches. Consequently, some of the virtual targets are observed by only one station. In the end, a total of 24 patches were manually extracted from the scans (Figure 12). For example, individual boulders where no internal deformation is to be expected.

Figure 12:

Point cloud of the southern, moving part of the cleft and manually selected patches (blue: station 1, yellow: station 2, green: both stations).

Next, the ICP quality is also checked between homologous patches of both epochs. As the number of points that are only observed from one station is high, this filtering step is crucial. Redundancy-free observations have a partial redundancy of r _i  = 0 and are not controlled by other observations. Thus, a possible gross error within the ICP matching cannot be detected in the network adjustment, inevitably leading to a wrong result in the deformation analysis. However, the MC simulation can minimize this risk. As a result of the MC simulation, three patches were rejected, as they are not suited for reliable matching. The visual analysis of the rejected patches revealed changes in their inner geometry due to erosion. In this case, the ICP matching is no longer reliable.

After calculating the virtual targets for the remaining 21 patches, pseudo observations are generated. Since the scanning total station was used to carry out additional network measurements, the pseudo observations are added to these. Thus, the number of object points in areas that are difficult to access is significantly increased by the additional virtual targets derived from laser scanning, and the spatial resolution is improved.

Also, the overall network configuration is improved. Due to the limited accessibility on Mt. Hochvogel, the design of the monitoring network is not ideal. However, the additional target points help to improve the network quality. Again, the partial redundancies before and after extending the network by virtual observations are analysed (Table 2).

The average improvement, however, seems to be small. This can be explained by the large number of points that are observed from only one station and thus do not raise the redundancy of the network. As their partial redundancy r _i  = 0, the average partial redundancy of the observations is still quite small. The improvement can clearly be seen in Figure 13, where only the observations with partial redundancy r _i  > 0 are shown. All observations with a partial redundancy of r _i  = 0 are not depicted.

Figure 13:

Redundancy component of the horizontal angle (Hz) observations for the original network (left) and after including the virtual points (right).

Finally, the rigorous deformation analysis is performed. The horizontal components of the resulting 3D deformation vectors are shown in Figure 14. Between September 2018 and July 2019, the gap opened by around 21 mm in horizontal direction. In vertical direction, almost no movement happened. The virtual points from the patches (orange points) give the same results as the marketed points (blue points), both in direction and magnitude. This confirms the functionality of the method also in a real use case.

Figure 14:

Horizontal component of the resulting deformation vectors between September 2018 and July 2019.

Furthermore, the enhanced spatial resolution enables the detection of more details regarding the movement. For example, the additional targets uncovered that a specific block, located at approximately 55 m East and 29 m North, is moving differently, ultimately leading to its collapse in 2021.

The average standard deviation of the deformation vectors is less than 2 mm and could also be slightly improved by adding the virtual targets (Table 3). All detected deformations are significant. Generally, the proposed workflow including the supervised point cloud matching ensures a high accuracy of the results.