Sensitivity quantification of spatially averaged displacement estimates in TLS-based geomonitoring

2.1 Theoretical background

A monitoring system (including both the deployed instruments and processing methods) is considered sensitive when it can detect a displacement with a specified probability of false alarm (significance level α) and a specified probability of detection (test power 1 − β), where β represents the probability is the probability of a missed detection (type II error) [26].

For geomonitoring, among the available techniques, the M3C2 algorithm [4] and its variants [27], 28] are the most widely deployed approaches to estimate geometric surface changes. For the sensitivity analysis of displacement estimates, M3C2 relies on the LoD, which is considered to be the sensitivity measure and can be used to indicate the minimum detectable changes [4]. The LoD is calculated by considering the variances of the two surface models or point clouds (reference and target) in the direction of the surface normal and the assumed co-registration error. The formula of LoD95 (for a 95 % confidence interval) [4] is given by:

where σ₁(d)² and σ₂(d)² denote the local surface variability (i.e. the combined measurement noise and the effect of the surface roughness) in each point cloud, which is computed from a subsection of the point cloud within a cylindrical neighborhood with diameter d around the query point. The additional parameter σ_reg quantifies an isotropic registration uncertainty. For an α other than 0.05, the factor 1.96 is replaced by the corresponding critical value.

Despite its broad applicability, the LoD shows three key limitations. (i) The standard LoD formulation does not account for the probability of a type II error, (ii) it assumes that all observations are uncorrelated, an assumption which does not hold for TLS points clouds, as discussed in Section 1, and (iii) the LoD computation assumes isotropic registration uncertainty (see Eq. (1)), implying uniform uncertainty in all directions.

To address the first limitation of the LoD, recent studies [29] have started to incorporate hypothesis testing as well-established statistical techniques for sensitivity analysis, with a specific focus on calculating the Minimal Detectable Bias (MDB) as a measure of sensitivity. The MDB was introduced by ref. [30] as part of a hypothesis testing-based approach for outlier detection, known as Baarda’s data snooping procedure. Although introduced for analysis of geodetic networks, it is a more widely applicable method [31], [32], [33]. For point cloud-based deformation analysis, we can define the MDB as a measure of the minimal displacement that has to occur between two epochs in order to be detected with significance level α while allowing for a probability of β that the test fails to detect it.

For the hypothesis testing of an observed displacement estimate d, we have the null hypothesis, H 0 representing the case that no deformation occurred, and the alternative hypothesis H A , which addresses a deformation scenario.

To compute the MDB, we have to consider the alternative hypothesis H A . The MDB is defined as the smallest displacement that leads to the rejection of the null hypothesis H 0 with a probability of at least 1 − β, assuming a given significance level α [34]. By approximating the test statistic of the signed point cloud displacement estimates by a normal distribution due to the Central Limit Theorem, the calculation of the MDB is simplified and leads to the following definition [35]:

where z_1−α/2 is the critical value (two-tailed) for a chosen significance level α (e.g. 0.05), z_1−β for a chosen power 1 − β of the test, and σ _d represents the standard deviation of the computed displacement. In contrast to the LoD formulation, where the uncertainty is derived from the individual surface variabilities of the two point clouds, σ _d in the MDB approach reflects an empirically estimated uncertainty based on the observed differences between epochs. Importantly, this formulation is agnostic to the specific method used to compute the point cloud differences.

In the literature adopting the MDB for point cloud-based deformation analysis (e.g. [29]), the parameter σ _d is still computed as in the case of LoD (term in the parentheses in Eq. (1)), based on the analysis of the local point cloud variance and an additional constant factor requiring prior knowledge and some simplifying assumptions. Hence, although using MDB as a measure of sensitivity eliminates the first of the abovementioned limitations, it still faces the second and third limitation of the standard LoD (does not account for correlations in the data and it assumes an isotropic registration error) if they are not explicitly addressed.

Consequently, the parameter σ _d in Eq. (2) must reflect all relevant uncertainty sources, such as registration error, instrument measurement noise, as well the effect of environmental, material, and geometric influences, as discussed in Section 1. Some of these effects can vary substantially across a point cloud, especially in the case of long-range TLS-based geomonitoring [36]. This also holds for the registration uncertainty, which depends on the estimated transformation parameters and the spatial position of the investigated points. Hence, this makes the third assumption, that the registration error is uniformly distributed, not valid.

Given these limitations and requirements for σ _d , we propose a direct estimation of the uncertainty σ _d based on empirical statistical analysis of locally averaged and sampled displacement estimates.

In the following section, we describe the workflow of the proposed empirical estimation of σ _d , which is later used for investigations in Section 4.

2.2 Workflow for empirical uncertainty estimation

The workflow follows the previously introduced ideas for the estimation of σ _d . In point cloud-based monitoring applications, displacements are most commonly derived from the comparison of reference and target point cloud epochs. These two-epoch displacement estimates are widely regarded as the standard case in geomonitoring and will also serve as the foundation for our subsequent analysis.

The general workflow is realized using the following four steps: (1) reference and target point clouds are subdivided into local segments, (2) displacements are computed between the corresponding segments, (3) displacements are aggregated within a defined local neighborhood, and (4) the empirical σ _d is derived per local segment through sampling of the aggregated displacements and computing the related statistics. This method of estimating σ _d will hereafter be referred to as ELSA (Empirical Local Sampling of Aggregated values). The resulting σ _d can then be used for calculating the Minimal Detectable Bias (MDB) as in Eq. (2).

Within this study, we present only one specific implementation of the workflow for one way of calculating displacement estimates from TLS point clouds acquired from a single viewpoint, which is typical for geomonitoring setups. Although most common algorithms, such as M3C2, compute displacements in the normal direction of the surface, we compute the displacements in the LoS (Line of Sight) direction for the following reasons. First, it aligns with the inherent nature of data acquisition by TLS instruments. Hence, no additional uncertainty is introduced, e.g. through the estimation of the surface normal. Second, it enables processing of the single-viewpoint point cloud in a 2D spherical image representation (a.k.a. range images or depth maps), a common representation in computer vision applications. This representation, in turn, (i) improves data processing efficiency by reducing the dimensionality from 3D to 2D and (ii) facilitates the use of established image-processing algorithms (primarily relevant for Section 4.3). Irrespective of this specific implementation, the core ideas of the proposed workflow can be easily adapted to process displacement estimates in any arbitrary or multiple directions and, hence, it can also be used, for example, for processing the M3C2-based displacement estimates.

The prerequisite for the workflow is that the target point cloud S target is already registered to the reference point cloud S ref , and any user-desired model-based corrections or reductions (e.g. atmospheric corrections) are applied beforehand. Any remaining registration inaccuracies or residual effects in the data will contribute to the final uncertainty and sensitivity estimate.

The upcoming subsections present the aforementioned specific implementation of the proposed workflow, and they follow the four steps presented. Section 2.2.1 describes the subdivision of point clouds into local segments (Step 1) and rasterization, which is necessary for representing 3D point clouds as spherical range images and for computing the LoS displacement estimates (Step 2); and Section 2.2.2 describes definition of local neighborhoods and displacement aggregation strategy (Step 3), as well as the computation of the parameter σ _d (Step 4).

2.2.1 Subdivision and rasterization

Let S be a point cloud represented in spherical coordinates, where each point s k ∈ S is defined as s k = θ k , ϕ k , r k . Here, θ _k denotes the horizontal angle, ϕ _k the vertical angle, and r _k the range for the k-th point.

Assuming that the aforementioned spatial influence factors (Section 1) have a similar or constant influence within a smaller space, we subdivide the point cloud into smaller, local segments. In general, this subdivision can be achieved in different ways, such as dividing the point cloud using a regular grid, segmenting along major break lines, or applying advanced segmentation methods like supervoxels [37].

Herein, with the spherical image representation of S in mind, we define the segments as regular subdivisions of the spatial domain spanned by the vertical angles Φ ∈ 0 , π and the horizontal angles Θ ∈ − π , π of S . Each segment is thus determined by two of its segment vertices (ϕ_min, θ_min) and (ϕ_max, θ_max).

In the following text, we present the remaining part of the implemented workflow for a single segment; this process has to be repeated for each segment. To generate the spherical image representation, we first define a regular grid of size n × m with horizontal and vertical angles θ _i and ϕ _j respectively, as depicted in Figure 1, where i = 1, …, n and j = 1, …, m. The step size of the grid corresponds to the scanning resolution Δφ used during the measurement. This ensures an approximately lossless data transformation, where each pixel approximately corresponds to one measured point. The grid coordinates are defined as:

(3) θ i = θ min + i − 1 ⋅ Δ φ , i = 1 , … , n

(4) ϕ j = ϕ min + j − 1 ⋅ Δ φ , j = 1 , … , m

Figure 1:

The rasterization of the acquired point cloud S involves interpolating the range values at each grid cell M_i,j based on the surrounding points, corresponding to the angular coordinates θ i , ϕ j of the predefined regular grid.

Subsequently, we define the two-dimensional matrix M ∈ R n × m that represents a range raster or image and is realized as an interpolated subset of S based on the previously defined regular grid. Each element M_i,j corresponds to an interpolated range value at the grid coordinates (θ _i , ϕ _j ), such that:

(5) M i , j = f int ( S , θ i , ϕ j )

where f int ⋅ denotes the interpolation function applied to S (we implemented the barycentric interpolation). This procedure is done for segments of both the reference and target point clouds ( S ref and S target ), using the same grid definition, resulting in corresponding range rasters M_ref and M_target respectively.

The LoS displacement estimates are computed as the range measurement differences between the reference and target range raster. This can now easily be computed by the cell(pixel)-wise difference of the interpolated range raster, resulting in the displacement (or range difference) raster D:

(6) D = M target − M ref

2.2.2 Spatial aggregation and variability estimation

As mentioned in Section 1, we assume that spatial influence factors exhibit similar or constant effects within each such segment, allowing us to also assume that the computed range differences can be represented as a stationary stochastic process. Based on this stationarity assumption, we estimate σ _d for a given segment using an approach inspired by the bootstrapping approach.

We follow a similar path as the LoD computation in M3C2; however, in M3C2 first computes the average within a local (cylindrical) neighborhood and then derives the difference, whereas we first compute the difference and then average within a local neighborhood.

For that, we first partition the raster D into non-overlapping tiles T of size t × t as shown in Figure 2. This tile size t serves the same purpose of neighborhood selection for the averaging as, e.g. the cylindrical neighborhood diameter d for M3C2.

Figure 2:

The aggregation is based on the range difference raster D, which is partitioned into non-overlapping tiles T_k,l based on the tile size t.

The tiles T_k,l ⊆ D are extracted from the segment, where k and l index the tiles along the horizontal and vertical dimensions, respectively. The indices k and l are defined as:

(7) k = 0,1 , … , n t − 1   l = 0,1 , … , m t − 1

A tile T_k,l is respectively:

(8) T k , l = D ( k t + 1 ) : ( k + 1 ) t , ( l t + 1 ) : ( l + 1 ) t

As the sensitivity of the displacement estimates is commonly increased through local averaging (i.e. computing local mean of displacement estimates), we also use the average as the aggregation function f _A (⋅) within our workflow and apply it to each tile:

(9) D t , A ( k , l ) = f A T k , l

However, any other function (e.g. the median or weighted average) can also be used depending on the specific application requirements. Therefore, we refer to this process as aggregation to maintain generalizability. Each aggregated value (one per tile) represents a final displacement estimate of improved sensitivity and can then be treated as an independent sample of a stationary process.

The estimate of the σ _d , is then calculated as σ d t , A (where each resulting d is dependent on tile t and averaging strategy A) across the set of aggregated values D t , A = D t , A k , l = f A T k , l , providing the estimate of the uncertainty within the segment:

(10) σ d t , A = 1 | D t , A | − 1 ∑ k , l D t , A k , l − μ D t , A

where |D_t,A| denotes the total numbers of tiles, and μ D t , A is the mean of the aggregated values:

(11) μ D t , A = 1 | D t , A | ∑ k , l D t , A k , l

While the underlying computations differ, the resulting uncertainty σ d t , A can still be compared to the combined uncertainty term σ 1 ( d ) 2 n 1 + σ 2 ( d ) 2 n 2 + σ reg in Eq. (1). By computing the uncertainty of aggregated (i.e., averaged) displacements by empirical analysis and sampling rather than simplified variance propagation, this approach provides a more realistic uncertainty measure by implicitly capturing any residual correlation in σ _d , with the number of aggregated tiles directly defining the effective degrees of freedom. In contrast, the standard error used to compute the uncertainty of averaged displacements assumes that all observations are independent.

4.1 Uncertainty estimates of averaged displacements

With the first analysis, we want to investigate how the uncertainty σ _d of spatially averaged displacement estimate with our method (ELSA) behaves with different aggregation tile sizes t, as the general assumption is that larger tile sizes will reduce the uncertainty respectively. Additionally, we compare this behavior to the standard error of the mean (SE), which represents the conventional approach for estimating uncertainty in aggregated data, e.g. in LoD. The SE for a tile size t is given by SE ( t ) = s / t 2 , with s = Σ ( d i − d ̄ ) 2 / ( t 2 − 1 ) 1 / 2 describing the sample standard deviation of the d _i displacement estimates, d ̄ their mean and t² is the number of individual displacement estimates contributing to a single tile. For tile size t = 1, the SE is equivalent to Eq. (10). The results are shown for the three segments of each dataset as described in Section 3. Figure 4 presents the results for the El Capitan dataset across the segments A, B, and C, and Figure 5 shows the corresponding results for segments D, E, and F in the Säntis dataset. In each plot, the solid line represents the uncertainty σ _d estimated by the SE, and the dashed line depicts our method ELSA, both as functions of aggregation tile size t.

Figure 4:

Standard deviation estimated using the standard error of the mean (SE) and ELSA for the three segments A, B, and C of the El Capitan dataset, as a function of the tiling size t.

Figure 5:

Standard deviation estimated using the standard error of the mean (SE) and ELSA for the three segments D, E, and F of the Säntis dataset, as a function of the tiling size t.

The results show that the standard deviation σ _d computed by ELSA is consistently higher than that computed using the SE for any given aggregation tile size. This trend persists as the tile size increases, whereas for uncorrelated measurements, one would expect σ _d to decrease proportionally to 1 / n . This highlights the notable influence of spatial correlation within the data. Unlike the SE, which assumes independence among observations and thus underestimates uncertainty of averaged displacements, our approach implicitly incorporates the effect of spatial correlations on the uncertainty of the averaged displacement estimates (as discussed in Section 1).

Furthermore, as demonstrated by ref. [15], spatial aggregation generally reduces the uncertainty of displacement estimates. However, this reduction is less pronounced than expected based on the typical SE estimate and does not tend toward 0 as expected based on the simplified variance propagation. In contrary, in many cases it reaches a certain plateau, after which further averaging is meaningless. It can be observed that our uncertainty estimates can be used to better analyze up to which point it is meaningful to average the displacements, and after which point the averaging will just reduce resolution without increasing sensitivity. For our particular datasets, we identify a suitable breaking point at tile size t = 20, which corresponds of averaging approximately 400 points. For t > 20 the improvements in σ _d are less than 10 % whereas the resolution continues with its quadratic loss.

The differences between the curves related to different segments underscore the importance of estimating σ _d locally. For example, the σ _d of segment F of the Säntis dataset is substantially higher compared to the other segments. One contributing factor is the larger distance between the instrument and the rockface, which scales the influence of the various effects, such as refraction or registration errors, for example. Additionally, the scene contains significantly more vegetation than segments D or E, further increasing uncertainty due to increased surface roughness and, therefore, reducing sensitivity. As we can observe, the uncertainty represented by SE, resp. LoD95, also accounts for these differences, as indicated by the different initial uncertainty (without averaging). However, it would fail to capture deviations with lower spatial frequency (e.g. due to spatial correlations over larger distances), as it relies only on the limited neighborhood size defined by the cylinder radius used for M3C2 computation. In contrast, ELSA leverages the entire segment size to compute uncertainty and can capture and reflect such deviations.

The overall offset between our uncertainty estimate σ _d and the SE could theoretically be corrected by incorporating the standard deviation of registration σ_reg as an additional uncertainty term as included in LoD95 (see Eq. (1)). However, this formulation is based on the assumption that σ_reg is constant across the entire scene, which does not align with the reality. This is well demonstrated in our results, as the offsets between the curves for different segments are different. Hence, a constant σ_reg cannot adequately account for this. Furthermore, even with σ_reg applied and potentially reducing the offset between σ _d and the uncertainty based on SE, the approach would not address the different slopes of the curves observed across varying aggregation tile sizes, esp. for smaller tile sizes (i.e. less averaging). Therefore, using the proposed ELSA method can provide better uncertainty estimates, and consequently improve the sensitivity estimates.

4.2 Validation of MDB derived from ELSA

In the previous section, we demonstrated a notable difference between the common and the proposed way (ELSA) of computing the empirical standard deviation σ _d of averaged displacement estimates and argued on the relevance of these differences. In this section, we verify the validity of σ _d estimates via the proposed ELSA method and the subsequent estimation of the MDB as the sensitivity metric by utilizing synthetically simulated ground truth displacements.

To evaluate the correctness of the computed σ _d and subsequently MDB, we simulate displacement estimates for two cases: Case 1 ( H 0 ) no displacement added and Case 2 ( H A ) displacement of magnitude MDB added. These simulated displacements were added to the range rasters D computed based on our real-world datasets presented in Section 3. A total of 10,000 synthetic datasets are generated for each case. The σ _d required for the MDB and critical value calculation is derived according to ELSA introduced in Section 2.2.2 using an aggregation tile size t = 20, which matches the spatial extent of the simulated displacements (for Case 2). The MDB for determining the magnitude of the simulated displacements is then calculated with a nominal significance level α = 0.05 and type II error probability β = 0.2.

In order to analyze the correctness of the MDB calculated in this way, we analyze all computed aggregated (averaged) displacements (i.e. pixels in the range raster D after averaging) and decide whether they correspond to the simulated displacements (alternative hypothesis H A – significant displacement occurred) or not (null hypothesis H 0 – displacement not occurred) based on the critical value computed using σ _d and α = 0.05 as introduced above. The expectation is that, if the calculated value of MDB is correct, we will be able to correctly classify displacements into significant and not significant with probabilities of type I and type II error of α = 0.05 and β = 0.2. During this process, we record the false positives (displacement predicted under Case 1) and false negatives (no displacement detected under Case 2). Based on these outcomes, we calculate the empirical values of α and β. We exemplary present the results in Table 1 for one of the segments used in the previous analysis (segment A, dataset: El Capitan).

The results from simulation-based validation show that the empirically calculated values for α of 0.056 and β of 0.195 fall closely to the nominal value of 0.2. This indicates that the MDB calculated using the proposed method correctly reflects the minimal displacement for the chosen type I and type II error probabilities, making it a suitable metric for assessing the sensitivity of the implemented deformation monitoring approach.

4.3 Sensitivity improvement by additional preprocessing

Although σ _d calculated using ELSA is a better representative of the true uncertainty of the averaged displacements than commonly used values (as demonstrated in the previous two sections), it still falls short in representing the full scope of the displacement uncertainties present in the real-world data. This is primarily because we use the standard deviation as a metric, which does not account for systematic effects such as constant bias and may underestimate uncertainty in cases where the distribution of values deviates significantly from normality (e.g. heavy tails). When we analyzed the specific segments of our real-world datasets, we observed such systematic effects in the displacement estimates. In order to improve the MDB based on ELSA a better representative of the sensitivity, and in order to improve the sensitivity of our specific displacement estimation workflow, we identified three preprocessing steps designed to address these issues, capable of mitigating their detrimental impact on the uncertainty estimate, and subsequently improving sensitivity. They are implemented as follows:

(i) Image-correlation based alignment By visual inspection of extracted segments, we observed smaller misalignment errors in the lateral direction between the target and reference point cloud. We assume this is primarily caused by inaccuracies in the registration of the point clouds. To address this, we implemented a local alignment of the point cloud based on the intensity information. For that, an intensity image of each individual segment is computed following the method outlined in Section 2.2.1; but with intensity replacing the range value. The resulting intensity image is then used in a frequency-based cross-correlation image alignment method, following the technique proposed by ref. [41]. The estimated subpixel shift is then applied in the respective angular direction to the segment of S target , accounting for the potential errors in the estimated translation and/or rotation parameters used to align S target with S ref . Consequently, M_target also becomes better aligned with M_ref, leading to generally smaller values in the resulting range difference raster D reducing eventual biases and larger distribution tails occurring, e.g. due to misaligned surface edges within the scene. Because the estimated shift is purely in the angular (lateral) directions, it does not absorb or compensate any true LoS displacement, which remains fully in D.

(ii) Edge removal Even when the two point clouds are well-aligned, the range differences at the surface edge tend to be significantly larger than the overall range differences in D. This is likely due to the effect of mixed pixels at surface edges, which commonly result in inaccurate range values. Consequently, we propose removing edge points from the range raster M_ref and M_target.

We remove the edge points by first generating a binary edge mask from the range raster M_ref, as we define the edges as regions characterized by rapid changes in local surface topography. To extract these changes, we first subtract the 2D linear trend from raster M_ref. Next, we compute the local gradients and convert them into local slopes following the approach described by ref. [42]. An edge mask is then generated by thresholding the slopes (we chose a threshold of 75°), classifying all values above this threshold as edges. To refine the mask, a morphological closing operation is applied, involving dilation followed by erosion of the binary edge mask raster. Finally, the edge mask is applied to the range rasters M_ref and M_target to exclude edge-influenced pixels. Since this step removes high-gradient (edge) pixels regardless of their LoS value, any genuine displacement in the non-edged (flat) areas is untouched.

(iii) Planar trend removal We also apply a planar trend removal to the range differences D. This process involves estimating a 2D trend by fitting a planar surface to the range differences. The corresponding values of the fitted surface are then interpolated for each D_i,j and subtracted from the range differences, effectively removing the surface trend.

The described preprocessing steps should ultimately reduce the uncertainty σ _d and simultaneously improve the representativeness of the related MDB estimate. We demonstrate the impact of these steps on the segment A from the El Capitan dataset. Their effect on the values of range difference raster D are shown in Figure 6, while Figure 7 highlights the corresponding uncertainty estimates of aggregated displacement σ _d .

Figure 6:

Results of the preprocessing steps applied to improve the bias and variability of the raw range differences (a). In the first step, the target dataset gets aligned by intensity-based image correlation using the intensity map (c) from the TLS, resulting in range difference (b) with significantly reduced bias and variability. In the second step, edge artifacts based on a local slope map computed from the range map (f) reduce big outliers caused by mixed pixels (d). The final step of removing a plane trend (e) reduces any remaining bias and further reduces the variability.

Figure 7:

Uncertainty estimate using ELSA with different preprocessing steps applied for segment A of the El Capitan dataset. The standard error (of the mean) is shown with the results from the proposed data-driven approach is given with the additional preprocessing steps (intensity-based image correlation alignment, removal of edge effects, and planar trend removal), to further improve the variability.

The raw range differences without any preprocessing (Figure 6a) exhibit a strong bias, with a mean of μ = −0.034 m. The standard deviation of the displacement is σ = 0.073 m. The application of the first preprocessing step (i) results in a range difference raster D, as shown in Figure 6b. This step corresponds to the lateral adjustment of the point cloud S target (or derivative, M_target) along the direction of horizontal and vertical angle measurements. Consequently, it does not incorrectly compensate for or absorb any deformations in the LoS direction. As can be observed, this alignment step significantly reduces the bias in the range differences to μ = −0.001 m. Additionally, this step substantially reduces the uncertainty of the range difference by 28 % to σ = 0.052 m, which means that also the sensitivity can be increased. This improvement is further corroborated by the corresponding significant reduction in the computed σ _d for any averaging (i.e. tiling) size as shown in Figure 7.

The next preprocessing step applies the edge removal. Using the extracted edge mask, we eliminate all pixels from the range difference raster D where the corresponding slopes exceed the threshold. As shown in Figure 6d, this step removes areas with larger range differences. The effect on the σ _d in Figure 7 is, therefore, most evident for cases with no aggregation or very small aggregation size ( < 20 ) , which is also shown in the uncertainty of the range differences with no aggregation of σ = 0.022 m. For larger aggregation sizes, the impact diminishes.

The final preprocessing step is the planar trend removal, with the results displayed in Figure 6e. This step eliminates any remaining planar trend and bias, further reducing uncertainty. The cause of this trend could, for example, be residual effects from registration in the line-of-sight (LoS) direction, which were not specifically accounted for in the lateral alignment of the first preprocessing step. Although the observed changes in this case were minor, they are still reflected in the slight reduction of σ _d in Figure 7. It should be noted that the application of this preprocessing step could partially absorb and remove the impact of deformations on the computed displacement values in the LoS direction, corrupting the implemented monitoring workflow. Hence, this last preprocessing step should only be applied in the cases for which it is known that the expected area under deformation is small compared to the chosen point cloud segment sizes.

Examining the resulting σ _d in Figure 7, we can assess the uncertainty before and after each preprocessing step and the corresponding effect when the displacement estimates get averaged. The corresponding uncertainties for the tile size t = 20 are listed in Table 2 along with the respective MDB computed by Eq. (2) with α = 5 % and β = 20 %.

The estimated averaged displacement uncertainty without preprocessing σ_d,w/o is 0.026 m, while with all preprocessing applied σ_d,w/ is reduced to 0.005 m. Consequently, the calculated MDB by Eq. (2) are MDB_w/o = 0.073 m and with all preprocessing steps applied, MDB_w/ = 0.014 m, resulting in a 80 % reduction of uncertainty and a corresponding increase in sensitivity. The most significant improvement is achieved by the image-correlation-based alignment, accounting for nearly 70 % of the total reduction with all preprocessing steps applied.

The investigated preprocessing steps demonstrated their potential for significantly improving the monitoring sensitivity. Based on this preliminary evidence, we have a first indication that the achievable sensitivity of detecting LoS displacements with long-range TLS in geomonitoring can approach approximately 1 cm. This finding will be further explored in the following section using simulated displacements for this specific point cloud segment, while future studies will need to investigate the generalizability of these indications.

4.4 Significant displacement detection using critical value based on ELSA

After assessment of the proposed method to estimate σ _d by ELSA, we now evaluate how these uncertainties affect displacement classification. For this purpose, we compare classification results based on σ _d obtained using the standard error (SE), our ELSA approach, and ELSA with additional preprocessing steps. We simulate three displacement scenarios with increasing magnitudes (0.005 m, 0.02 m, and 0.1 m) and add this to the real displacement estimates of segment A of the El Capitan dataset.

Figure 8 summarizes the results, where each row corresponds to one of the three displacement magnitudes, and each column shows the classification outcome for each method used to estimate σ _d . The red dashed line outlines the patch in which the (true) displacement was introduced. Pixels classified as not showing significant displacement are shown in grey, whereas pixels classified as significantly displaced are shown in blue.

Figure 8:

Classification results of the displacements based on the uncertainty estimate based on three methods (columns) across three simulated displacement sizes (rows). The first column represents the classification based on the estimation computed by the standard error (SE), the second is our method ELSA, and the third is ELSA with the additional preprocessing steps. The outline of the simulated displacement patch is indicated by the red dashed line.

First, we examine the case where the uncertainty is estimated using the SE approach. Using the standard deviation derived from the SE method for significance testing leads to nearly the entire patch being classified as significantly deformed across all simulated displacements, resulting in a large number of false positives. Overall, these results demonstrate that using the SE to estimate uncertainty leads to poorer classification performance due to its unrealistic small estimated measurement uncertainty σ _d , which is in accordance with the data presented in Figures 4 and 5.

When the uncertainty σ _d is estimated using ELSA without preprocessing (middle column), we observe a significant improvement in the classification results by a substantial reduction of false positive classification outside the simulated displacement area. The displacements of size 0.01 m and 0.02 m are not classified as significant, while the displaced patch of 0.1 m is classified as significant. This aligns well when compared to the respective σ _d and resulting MDB for the case of no preprocessing in Table 2. The artifacts, appearing as falsely detected displacements outside the region of simulated deformations (consistent blue regions on the left-hand side), are caused by outliers, primarily located at the edges of the rock surface.

When we combine the uncertainty estimation of ELSA with the previously discussed preprocessing steps, we are capable of fully correctly classifying larger displacements ( > 0.01 m ) . The edge removal as part of the preprocessing steps significantly reduced the observed artifacts, although it did not completely eliminate them. While a lower, more aggressive slope threshold could potentially remove these outliers entirely, it may result in excessive data removal in areas with naturally rougher surfaces. Future work is needed to investigate the optimal definition of this threshold and develop a systematic method for its determination.

To evaluate the classification performance, we computed standard performance metrics for binary classification – recall, precision, accuracy, and F1 score – based on the pixel-wise comparison with simulated displacement ground truth [43]. The results are presented in Table 3. Recall measures the ability to classify actual displacements correctly. A higher recall corresponds to fewer missed detections. When we base the classification on the uncertainty measure using the standard error (SE) and our methods (ELSA), both exhibit good recall. However, relying solely on the recall of SE can be misleading, as the classification based on SE also labels a large portion of non-displaced areas as significantly displaced, as previously observed. The precision metric reflects the proportion of correctly detected displacements among all locations that were classified as significant displaced, with higher precision indicating fewer false positives. The results based on the SE method demonstrate very low precision (0.02–0.06), whereas only the ELSA method with preprocessing and for displacements of at least 1 cm achieves substantially higher precision (0.39–0.74), reaching levels that may be practically useful. Accuracy provides an overall performance measure of correct classifications, where ELSA outperforms the classification based on the SE. To conclude, these results demonstrate that using σ _d computed using ELSA provides a more suitable base for calculating the critical value and subsequent hypothesis testing with the goal of detecting significant deformations within the monitored scene.

4.5 Spatial assessment of MDB using ELSA

Expanding on the previous evaluation, we aim to investigate the spatial extent of the classified displacement and to get further indication of the achievable sensitivity of LoS displacement estimates using the proposed workflow.

To achieve this, we use the same dataset and setup as in Section 4.4. We iteratively simulate displacements of the predefined patch (outlined in Figure 8) with increasing magnitudes from 0.001 m to 0.1 m, by 0.001 m step size. After each iteration, we classify each pixel as significantly deformed or not based on their displacements. Then, we check if any pixel is newly correctly classified as significantly displaced within the patch compared to the previous iteration. If yes the respective magnitude of the simulated displacement is stored for that pixel location. Hence, in the end, we get a heat map-type visualization of the smallest displacements that we detected as significant at a particular pixel location, which indicates the smallest deformation that we could have detected given a critical value (α of 5 %) as a threshold. Additionally, we evaluate the case without any simulated displacement to identify pixels that are incorrectly classified as significantly displaced, labeling them as outliers.

The results of this evaluation are presented in Figure 9, where each pixel is colored according to the smallest displacement magnitude that leads to its classification as significantly deformed, given it is not an outlier (grey). Again, the results are presented for three cases, for the classification using critical value based on SE and based on ELSA, with and without preprocessing. The results for the SE reveal a high rate of outliers, with nearly the entire patch being classified as significantly displaced, even in the absence of any simulated displacement. This again confirms that critical value based on SE is overly optimistic (too small) and can lead to a large type I error rate. In contrast, our method starts classifying displacements in the range of approximately 0.06 m–0.08 m as significant, which is in accordance to the results presented in previous section for the case of analyzing the displacement estimates without preprocessing. When additional preprocessing steps are applied, the results show a further improvement in sensitivity, with significant displacements being detected at much smaller magnitudes, ranging between 0.01 m and 0.02 m. This highlights the substantial enhancement for the classification of small-scale displacements achieved through preprocessing a local point cloud segment using the proposed workflow. Finally, this indicates that sub-centimeter to a few centimeter displacement detection is achievable for long-range TLS geomonitoring in the case of LoS displacement monitoring using the proposed monitoring setup and workflow.

Figure 9:

Result of the classification result using uncertainty estimates by (a) standard error, (b) our data-driven ELSA (tile size t = 20), and (c) ELSA with the additional preprocessing steps. Per pixel, the minimal simulated displacement is shown, which led to classification as significantly displaced. In grey, the areas that were already classified as significantly displaced when no deformation was simulated. The red dashed line outlines the simulated displacement patch.