INFLUENCE OF RANGING UNCERTAINTY OF TERRESTRIAL LASER SCANNING ON CHANGE DETECTION IN TOPOGRAPHIC 3D POINT CLOUDS

Terrestrial laser scanners are commonly used for remotely sensing natural surfaces into 3D point clouds. Time series of such 3D point clouds can be analysed to gain information of surface changes that are induced by Earth surface shaping processes. The atomic unit in time series analysis is a bitemporal change detection and quantification. This should involve an estimation of the minimum quantifiable change, the Level of Detection, to separate signal from noise, e.g. stemming from the measurement. To enable such an estimation through error propagation, a model of the sensing instrument’s measurement uncertainty is required. In this work, we present an investigation on the ranging component of terrestrial laser scanning on this uncertainty and its influence on 3D distances between point clouds of two epochs. Specifically, we analyse the effects of incidence angle, intensity and range for different object materials, and make additional considerations with respect to waveform information returned by the sensor. We estimate a model for the rangefinder uncertainty of a terrestrial laser scanner and apply it on experimental data. The results show that using a sensor-specific model of ranging uncertainty allows an appropriate estimation of the Level of Detection. At a range of 60m and a rotational displacement of 10◦, this Level of Detection ranges between 0.1mm to 1mm for a white and a grey surface and up to 5mm for a black surface. The completeness of the detection of significant change ranges from 60.2% (black) to 89.8% (grey) for the proposed method and from 65.5% to 88.9% for the baseline, when compared to tachymeter measurements. The similarity between the results is expected and suggests the validity of error propagation for the derivation of the Level of Detection.


INTRODUCTION
The availability of high-resolution topographic 3D point cloud time series has enabled the detection and quantification of changes in surface geometry (Eitel et al., 2016;Fey and Wichmann, 2017;Zahs et al., 2019). These 3D point clouds are unordered, irregular sets of points in 3D space acquired, e.g., by terrestrial laser scanning (TLS). TLS sensors are used to repeatedly sample the Earth's surface, but the exact locations of the sampled points differ from scan to scan and therefore from epoch to epoch. To obtain a reliable measure of surface change, points are aggregated to local surface models. Often, a locally planar surface is assumed. The most developed method to quantify distances between two point clouds for geographic applications is the multiscale model-to-model cloud comparison (M3C2, Lague et al., 2013), which directly compares two point clouds by such a locally planar model.
To separate real change from change values stemming from measurement and other noise, statistical tests can be employed. Lague et al. (2013) estimate the required second-order moments by modelling the laser rangefinder uncertainty 1 from the distribution of points in the local point cloud neighbourhood. A coregistration term accounts for the misalignment of the datasets of the two epochs. James et al. (2017) present an extension to this uncertainty model for a photogrammetric point cloud by using "precision maps" which map a spatially variable coregistration uncertainty to the dataset. Still, the per-point uncertainty is only estimated from the data and its fit to a planar model.
On the sensor side, multiple studies have investigated the influence of parameters like incidence angle (Soudarissanane et al., 2007;Kersten et al., 2009;Zámečníková et al., 2014), colour (Clark and Robson, 2004;Mechelke et al., 2007) and intensity (Pfeifer et al., 2007;Wujanz et al., 2017Wujanz et al., , 2018 on the ranging precision of laser scanners. Furthermore, Fey and Wichmann (2017) have quantified how range and incidence angle influence the positional uncertainty of 3D laser scanning points, and have applied their findings to geomorphic change detection in alpine terrain.
These studies suggest that a combination of the models derived from the sensing process and from data analysis could lead to a better understanding of the distribution of measurement uncertainties throughout a dataset. Knowledge of these uncertainties can then be used to estimate the resulting uncertainty in a change analysis, allowing for a clear separation of measurement noise and change signal.
We therefore perform an experiment to investigate the influence of material, range and incidence angle on intensity, pulse shape deviation and, most prominently, ranging uncertainty (i.e., ranging precision), with the goal of detecting change in a bitemporal dataset acquired with TLS. Subsequently, we create a model that allows the estimation of the ranging uncertainty for each point of the point cloud individually. We show how to use this estimation in a change detection analysis and present the applicability of our method by comparing our results to the state-of-the-art method of M3C2.

Experimental design
To evaluate the influence of different object and sensing properties, we conducted an experiment under controlled conditions. This consisted of scanning a locally planar wooden board at different ranges and different incidence angles. Additionally, we covered the board with dull plastic foil (black and white) and toned paper (grey), resulting in different reflectance values. The reflectance at approximately 0 • incidence angle was calibrated by using a target with known reflectance (SphereOptics Zenith LITE with a reflectivity of 92.5 % at λ = 1550 nm) at the same range. The resulting values are listed in Table 1 The wooden board was mounted in an upright (vertical) position, and subsequently rotated from normal to the laser beam (0 • ) to angles of up to 60 • , in increments of 10 • . Via this rotation, we cover a number of different incidence angles of the laser beam. We define the incidence angle as the angle between the inverted beam vector and the local surface normal vector oriented towards the scanner. Furthermore, the rotation acts as controlled displacement for the surface change quantification. The scans of all seven rotations were repeated for ranges of about 30, 60, 90, 160 and 220 m, resulting in a total of 35 scans.
As a reference for the displacement quantification, eight retroreflective targets were installed on the board and measured with a tachymeter (Leica TCRA705power). During the TLS scans, we covered the targets to avoid any interferences. Since displacement only regards relative movement, the total station was not tied in with the TLS coordinate system, but we assume the scale of both systems to match. The experiment was carried out under controlled atmospheric conditions in a basement with artificial lighting, and with two different TLS instruments (RIEGL VZ-400 and RIEGL VZ-2000i). The measurement setup is shown in Figure 1.
To ensure that the recorded noise stems from the sensor and not from roughness in the surface itself, we used a sanded OSB-3 wooden board, which has a tolerance in thickness of 0.3 mm according to EN300 (DIN, 2006). This is an order of magnitude lower than the expected ranging precision given by the TLS manufacturer, which is 3 mm at 100 m range for both employed TLS systems (RIEGL LMS, 2019). Figure 2 shows the board in the measurement environment.

Modelling range uncertainty
We evaluate the ranging uncertainty by calculating the offset of individual point measurements from the locally best fitting plane. First, the 3D points on the board are manually segmented from the full point cloud and separated into the three surface materials. A local neighbourhood of 10 points (corresponding  to a circle with a radius of approx. 3 cm at 60 m range with the used scan parameters) is then used to define the local plane for each point using a least-squares optimization, which allows us to disregard potential low-magnitude but large-scale deformations like a bend of the board.
We then project the quality of the plane fit (σ0, i.e. standard ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume V-2-2020, 2020 XXIV ISPRS Congress (2020 edition) deviation of the points in the local neighbourhood in direction of the plane normal vector) onto the beam vector of the neighbourhood centroid to get an estimate of the ranging precision. Since our ranging distance is large compared to the extent of the board, we average the ranging precision per surface material.
In addition to the range and angular measurements of the TLS, we record the signal intensity Int and the pulse shape deviation Dev. In accordance with RIEGL LMS (2012), we define the intensity to be a unitless [dB] measure of the amplitude of the returned signal without any corrections corresponding to range or incidence angle normalized to 16 bits. The pulse shape deviation is a unitless measure of how similar the outgoing and the incoming pulses are shaped, and can therefore be related to the quality of the range measurement.
With the 35 estimates of ranging uncertainty made at five different ranges and seven different incidence angles respectively, we use a least squares method to create a model of linear combination for the estimation of the ranging uncertainty. The inputs to this model are the cosine of the incidence angle, the deviation and the recorded intensity. In accordance with Wujanz et al. (2017), we rely on the range equation (Eq. 1) to model the change in range implicitly via the change in intensity. To test this assumption, we add the range as a separate input and analyse the correlation between the parameters.
where Pr = received signal power Pt = transmitted signal power Dr = diameter of receiver aperture R = range from sensor to target β = beam divergence ηsys = system efficiency ηatm = atmospheric transmission factor σ = target cross section (Jelalian, 1992)

Change detection and quantification
We first quantify change between point clouds by using the common M3C2 point cloud distance measure (Lague et al., 2013). It aggregates points within a local neighbourhood formed by a cylinder oriented along the local normal vector of the reference epoch point cloud. The points within these local cylinders are then projected onto the cylinder axis for each epoch separately, where the difference between their mean positions is used as a distance measure. The standard M3C2 then uses the standard deviation of the points along the cylinder axis as a precision measure.
In our experiment, we use the board at a 0 • incidence angle, with the board normal to the laser beam, as reference epoch, and quantify change orthogonal to this epoch, as shown in Figure 3.
Subsequently, we detect significant change as a binary label using a two-sided t-Test. This test can be reformulated to find the minimum quantifiable change, referred to as the Level of Detection (LoDetection, Eq. 2). Here, σi refers to the standard deviation of the points along the cylinder axis and ni to the number of points within the cylinder, for epochs i = 1 and i = 2, respectively. The term reg. is an additional constant coregistration uncertainty, which is neglected in this experiment because of the static setup, since the scanner was not moved during the acquisitions at a constant range.
Instead of applying Equation 2, we propose to propagate the ranging uncertainty from the TLS to the individual laser points.
Since it extends the M3C2 algorithm by error propagation 2 , we refer to it as M3C2-EP.
We estimate a full covariance matrix for the local neighbourhood centroids of each epoch. Since we are looking for significant change in the specific direction normal to the planar surface, we subject this matrix to a projection onto the normal vector n and derive a Level of Detection using a test statistic, which is distributed according to Hotelling's t-squared distribution (Rencher, 1998). This distribution can be related to an F-distribution, giving the expression in Equation 3. This corresponds to a two-sided test of multivariate normally distributed means with p = 3 degrees of freedom.
Here,Ĉ is the pooled covariance matrix (Eq. 4), n1 and n2 are the number of points, and C1 and C2 are the covariance matrices for the respective epochs.
To achieve full 3D per-point uncertainties in the form of covariance matrices for each laser point, we assume the scanner to have an angular standard deviation defined mainly by the laser beam divergence. This neglects the influence of the angular measurement resolution, in contrast to the work by Lichti and Jamtsho (2006). For the specific scanner used, the beam divergence accounts for most of the angular uncertainty. According to the datasheet, the beam divergence of the RIEGL VZ-2000i is 0.27 mrad (0.3 mrad for the RIEGL VZ-400), and is defined at the 1/e 2 points of the energy distribution. Assuming a constant reflectance of the board within a single laser footprint, we use a quarter of this beam divergence as standard deviation in scan and yaw angles: σScan = σY aw = 0.27 mrad/4 = 0.0675 mrad.
Since the scanners were not moved when changing the incidence angle, we disregard any terms concerning the coregistration, leading to a simple functional model for the error propagation, as shown in Equation 5. Here, ϕ represents the yaw angle (horizontal), θ the scan angle (vertical), r the range and (x, y, z) the cartesian coordinates for each point i.
(x, y, z) T i = ri · (cos ϕi sin θi, sin ϕi sin θi, cos θi) T By comparing the quantified change with the reference change measured by tachymeter, we can evaluate the detection metric in form of a confusion matrix, as shown in Table 2. We assume change to be significant (i.e. detectable) if the change calculated from the tachymeter measurements is greater than or equal to the estimated LoDetection of the respective method, and to be indistinguishable from noise otherwise. This significance is then compared to the one from the statistical test using the quantified change by means of correctness and completeness (Eq. 6). Here, "TP" refers to a change calculated by tachymeter measurement that is larger than the associated LoDetection and is estimated as significant from the TLS point cloud. "FP" is a change that is estimated as significant from the TLS data, but is not larger than the LoDetection when calculated from tachymetry. Similarly, "TN" and "FN" refer to changes estimated as not significant and are smaller and larger or equal to the LoDetection (when change is calculated from tachymetry), respectively. The LoDetection is always derived from the TLS data.

RESULTS AND DISCUSSION
We first present dependencies between different factors of influence by visual analysis of respective scatterplots in Section 3.1. Subsequently, a ranging uncertainty model is estimated from the data (Section 3.2), which is then used and validated in error propagation for the detection of significant change in Section 3.3.

Dependencies between ranging uncertainty, range, incidence angle, intensity and pulse shape deviation
In order to find a functional model describing the ranging uncertainty, we examine how different measures provided by the rangefinding system itself (intensity, range, and deviation) and by the data (incidence angle) correlate with each other and with the ranging uncertainty. Figure 4 shows a scatter plot between ranging uncertainty and range for the VZ-2000i system at a constant incidence angle of about 0 • . The curve shows a clear trend for all three materials that the ranging uncertainty increases with larger ranges. Similarly, the ranging uncertainty decreases for increasing intensities, as shown in Figure 5.  Similar observations can be made for intensities and incidence angles, where the intensity decreases with an increasing incidence angle (Fig. 7). However, the quality of the plane fit and also of the projected ranging uncertainty exhibit different behaviours for different surface materials. For the black surface, the ranging uncertainties increase with increasing incidence angles as would be expected and was also shown in previous studies (e.g. Soudarissanane et al., 2007), even if the effect is more prominent with phase scanners than with time-of-flight scanners like the RIEGL sensors (Kersten et al., 2009). But for the grey and white surfaces, the uncertainties decrease for incidence angles above about 20 • (Fig. 8). This behaviour is similar for all ranges, and suggests that the type of material influences the ranging precision going beyond just the strength of returned signal, i.e. that the material is not behaving like a Lambertian reflector. This suggestion is further supported by the fact that especially the black foil has specular reflective properties as visible in Figure 2. For the black surface, only three data points are available, since no points were recorded for incidence angles above 20 • at the range of 60 m.  The data recorded with the RIEGL VZ-400 shows a similar pattern (Fig. 9), with the maximum ranging uncertainty shifted towards higher incidence angles. Still, the ranging uncertainty decreases for higher incidence angles. However, this can be explained by the projection of the plane fit onto the beam vector reducing the uncertainty by the cosine of the incidence angle (where the angular uncertainty starts to contribute more to the total positional uncertainty). We also plot the quality of the plane fit against the incidence angle (Fig. 10). This shows an almost linear relationship, which fits the theoretical expectation. For the RIEGL VZ-2000i, this relation is not that prominent, but the decrease of the ranging uncertainty with increasing incidence angles (Fig. 8) is partly explained.   In addition to the polar measurements and the intensity, the scanners deliver information on how the received pulse shape deviates from the outgoing one. Plotting this deviation as a function of the incidence angle shows a similar pattern to the ranging uncertainty (Fig. 11), suggesting that the deviation can be used, at least in part, to explain the behaviour of a decreasing ranging uncertainty for high (flat) incidence angles.
Similarly, a plot of the ranging uncertainty as a function of the pulse shape deviation shows higher ranging uncertainties with a larger pulse shape deviation for the grey and the white surfaces as expected, but an inverted effect for the black surface (Fig. 12), that we again attribute to the non-Lambertian properties of the material.

Ranging uncertainty model
The observations presented in Section 3.1 indicate that a single model is not sufficient to adequately describe the ranging uncertainty for all three surface materials. Instead, each surface individually exhibits its own material properties which influence the ranging uncertainty. Therefore, we estimate a function σr = f (. . .) separately for each surface material. This also covers the difference in colour, i.e. reflectance at λ = 1550 nm, as described by Clark and Robson (2004).The TLS system further has specific characteristics, especially concerning the (uncorrected) intensity signal. Therefore, individual models have to be created for different TLS models (Wujanz et al., 2017(Wujanz et al., , 2018. We choose f = a + b · Int with unknown parameters a and b and the recorded intensity Int as a starting point for our model. In contrast to Wujanz et al. (2018), we do not model an additional parameter in the exponent of Int, because this did not lead to a converging solution with our data. Still, we expect this function to adequately model the influence of different measurement ranges, as these explain much of the variation in the intensity (see Fig. 6). Additionally, we consider the cosine of the incidence angle cos(ϕ) and the pulse shape deviation Dev as additional linear influences with factors c and d. The full model is stated in Equation 7.
The parameters for the different surfaces are estimated for the RIEGL VZ-2000i using a least-squares method on all available data. The results are listed in Table 3  An analysis of the correlation matrix of the estimated parameters for the white surface shows a maximum absolute correlation of 0.46 between parameters a (the constant term) and d (the deviation dependent term). When including the range as an additional linear parameter, the respective coefficient highly correlates (0.92) with parameter b, the amplitude dependent term.
The models for the other surfaces exhibit slightly higher correlations, especially the one for the black surface, where a and c (the constant term and the incidence angle-dependent term) correlate with a coefficient of 0.95, even without including the range dependent term in the model. This is not surprising, as only data for low incidence angles is available, limiting the model's ability to separate the influence of incidence angle from a constant term, validating our assumptions.

Error propagation and assessment of significant change
We apply the model presented in Section 3.2 to the original point cloud, generating a stochastic point cloud in accordance with Wujanz et al. (2017). In order to investigate the ranging uncertainty as the major contribution in the uncertainty of the change quantification, we calculate changes of the rotated board with respect to the 0 • position, and look for changes along the normal vector of this position (see Fig. 3).
Because of the noise in the intensity and deviation signals, the resulting per-point covariance varies in small neighbourhoods. However, since the M3C2 algorithm includes an aggregation, these variations are averaged out when testing for significant changes between the two epochs. Figure 13 shows the per-point covariance by means of the semimajor axis of the error ellipsoid for one scan after this aggregation step. The different materials show specific patters of their respective ranging uncertainty models. Since the incidence angle is assumed constant in this example and local variations in intensity and deviation are averaged out within the search cylinders, the only remaining influence is the range -which is slightly (i.e., by a few mm) larger at the edges of the board. Therefore, Figure 13 shows the dependency of the uncertainty by range, only changing by a few micrometers across the board. We show the differences between the reference change and the quantified change for an increase of incidence angle from 0 • to 10 • in Figure 14. The reference change was calculated by fitting 2.5D quadrics to the tachymeter measurements in the respective epochs, sampling points on these surfaces and then applying M3C2. The same sampled points are then used as neighbourhood query centers for the subsequent, point cloud-based analysis. By comparing the quantified change to the estimated LoDetection, significant change is identified. The quantified change is calculated using the M3C2 algorithm with a search radius of 5 cm. The normal vectors are taken from the artificially sampled point cloud at approximately 0 • incidence angle.
Two different charts of significant change are shown in Figure 15 for a part of the board. The detection of significant change between the standard M3C2 and our method with error propagation is practically equal. Additionally, the difference of the LoDetection is shown in Figure 16, where our method estimates a much more homogeneous LoDetection for all surfaces than the baseline model. By comparing the estimated LoDetection to the change quantified by the tachymeter, we estimate a Figure 14. Difference between change estimated from the tachymetry measurements and the M3C2 change quantification. Differences up to 5 mm appear, especially in the black surface. truth for significance of the change, as presented in Section 2.3. Correctness and completeness for both methods are listed in Table 4. Comparing the state-of-the-art M3C2 model with ours shows very similar values.  While our model shows a generally lower LoDetection, it is slightly higher than the baseline on the black surface material. The baseline method suffers from effects of locally lower point densities especially at the edges and close to data gaps.

CONCLUSIONS
In this contribution, we present an experiment investigating how the precision of the laser rangefinder of a terrestrial laser scanner varies for different surfaces, ranges, and intensities. We propose a novel uncertainty model that makes use of object properties (surface material), scan geometry properties (incidence angle), and properties of the measurement itself (amplitude and pulse shape deviation).
We subsequently apply this uncertainty model to estimate perpoint uncertainties which are used in a change analysis to identify significant change, and compare this to the state-of-theart change detection of M3C2. The herein presented method requires more information about both the data acquisition (system and geometry) and the sensed objects (i.e., which model of surface material to use) than the baseline method. However, it allows a more rigorous analysis of the errors involved and can therefore provide a better prediction of change significance. The model itself may be derived from the data, given that a careful planning of scan positions is undertaken.
In the context of time series analysis of topographic data, the presented method is a two-fold improvement of the state of the art: (1) Since the application of the rangefinder uncertainty model does not have a requirement for planar objects, observations of natural objects with curved or ragged surfaces are not limited by planarity constraints. This also applies for data gaps, where planarity measures suffer from the inexistence of data.
(2) Including the sensing process in the analysis allows the rigorous combination of multiple data sources, e.g. combining airborne laser scanning data with terrestrial laser scanning data. The individual uncertainty components of each sensing method can be modelled and then propagated to the change quantification.
While the findings remain to be shown with real-world, topographic data, we assume that an improved uncertainty model allows the detection of changes that were previously discarded as non-significant, especially on non-planar surfaces. With increased acquisition frequencies and a trend towards 4D point cloud analyses using more than 2 epochs (e.g. Anders et al., 2020), these considerations become increasingly important.
With respect to the model generation, further investigations should show whether (1) models from experiments like the one shown in this paper can be transferred to real data or (2) the scanner's uncertainty can be estimated from patches of overlapping point clouds acquired from different scan positions.