VIRTUAL LASER SCANNING OF DYNAMIC SCENES CREATED FROM REAL 4D TOPOGRAPHIC POINT CLOUD DATA

: Virtual laser scanning (VLS) allows the generation of realistic point cloud data at a fraction of the costs required for real acquisitions. It also allows carrying out experiments that would not be feasible or even impossible in the real world, e.g., due to time constraints or when hardware does not exist. A critical part of a simulation is an adequate substitution of reality. In the case of VLS, this concerns the scanner, the laser-object interaction, and the scene. In this contribution, we present a method to recreate a realistic dynamic scene, where the surface changes over time. We first apply change detection and quantification on a real dataset of an erosion-affected high-mountain slope in Tyrol, Austria, acquired with permanent terrestrial laser scanning (TLS). Then, we model and extract the time series of a single change form, and transfer it to a virtual model scene. The benefit of such a transfer is that no physical modelling of the change processes is required. In our example, we use a Kalman filter with subsequent clustering to extract a set of erosion rills from a time series of high-resolution TLS data. The change magnitudes quantified at the locations of these rills are then transferred to a triangular mesh, representing the virtual scene. Subsequently, we apply VLS to investigate the detectability of such erosion rills from airborne laser scanning at multiple subsequent points in time. This enables us to test if, e.g., a certain flying altitude is appropriate in a disaster response setting for the detection of areas to immediate danger. To ensure a successful transfer, the spatial resolution and the accuracy of the input dataset are much higher than the accuracy and resolution that are being simulated. Furthermore, the investigated change form is detected as significant in the input data. We, therefore, conclude the model of the dynamic scene derived from real TLS data to be an appropriate substitution for reality.


INTRODUCTION
Virtual laser scanning (VLS), or the simulation of laser scanning surveys, is a method to acquire point clouds instead of using real data. Compelling advantages of VLS include the possibility to simulate vast amounts of data at a very low cost, the generation of training data for machine learning with perfectlyknown reference labels, or the simulation of hardware that does not yet exist . VLS may also be used for survey planning (Hämmerle et al., 2017) or validation of signal analysis methods (Richter and Maas, 2022).
A prerequisite for the simulation of realistic point clouds is the existence of an appropriate 3D model of the scene to be scanned. If the target is VLS of the Earth's topography, such a model can be created from scratch in 3D modelling software, by simulation of physical processes (e.g., drip erosion, debris flow, or landslides), or by transfer from existing datasets, usually recorded at a higher quality than what is aimed for in the simulation. For example, the simulation of a UAV-based survey can employ a model derived from terrestrial laser scanning (TLS) data (Weiser et al., 2021).
For the simulation of surface change, physical models are often employed (Chang et al., 2015). These models build on the knowledge of the subsurface dynamics, e.g. deep-seated landslides. These models are difficult to employ, if underlying causes for the surface change are not known sufficiently well * Corresponding author or no such models are available, or if multiple geomorphic processes act on the same position resulting in a signal superimposition at the surface.
We aim to transfer a change form detected in a real 4D point cloud dataset (three spatial and one temporal dimension) and apply its morphometric properties and their temporal changes to a different surface, represented by a mesh model. This mesh model can subsequently be scanned using VLS, and the resulting point cloud can be analysed for the detectability of the modelled change with the simulated acquisition.
The use case presented here transfers an erosion form from a TLS dataset to a virtual scene representing a high-mountain setting and then investigates whether the surface changes can be detected in data acquired with airborne laser scanning (ALS). In a disaster response situation, such a scenario can be imagined when it is necessary to quickly assess slope stabilities of large affected areas, e.g., after heavy rainfall, where immediate danger of debris flow is present. A lower flight altitude with lower flight speeds results in a more detailed sampling of the surface, as the spatial resolution is increased. However, with higher altitudes, more area can be covered per time unit. It is therefore desirable to find an optimal altitude, where change can be reasonably well detected, yet as much area as possible can be covered. Additionally, this altitude varies depending on the magnitude of the observed change, which is assumed to increase over time. Stakeholders need information on (a) the point in time when a change form becomes identifiable for (b) different choices of flight altitudes. We, therefore, develop an approach to apply VLS to a multitemporal dataset.
The objective of this paper is to generate virtual dynamic scenes, adding to the established methods of transferring phenomena from existing datasets to virtual scenes by creating mesh models from processed TLS point cloud time series. We focus on terrain surface changes induced by mass movements, i.e., ignoring effects from vegetation such as trees or tall grass.
Section 2 presents the 4D point cloud dataset used in this study. Subsequently, the methods of processing this 4D dataset by spatial and temporal smoothing, clustering, extraction of a single change form, transfer to a virtual scene, and simulation of laser scanning are explained (Section 3). We present (Section 4) and discuss (Section 5) the results and draw conclusions for real use cases (Section 6).

DATASET
The 4D point cloud dataset used to extract the change form is a near-continuous TLS time series, acquired over one week in Vals, Tyrol, Austria, in August 2020. The time series is part of a larger investigation over multiple months (Schröder and Nowacki, 2021). The scene was monitored after a rockfall event, which caused a large debris cone at the foot of the slope. The time-series features surface activity due to subsequent anthropogenic works as well as erosional mass movement patterns of rill erosion.
The time series was acquired between 2020-08-20 00:00 and 2020-08-25 20:00 (CEST, formatted YYYY-MM-DD HH:MM), every two hours, resulting in 71 epochs. We aligned all subsequent epochs to the first epoch ("null epoch") by an ICP algorithm (Besl and McKay, 1992) operating on manually delineated stable parts outside of the observed slope. We use the OPALS software (Pfeifer et al., 2014, v 2.3.2) for this ICP alignment. After the alignment, vegetation points and outliers were removed from the dataset by applying a statistical outlier filter (k=8, multiplier=10.0; Rusu et al., 2008) and an SMRF filter (cell size=0.5 m, slope=2; Pingel et al., 2013), as well as a filter on the waveform deviation (≤50), all implemented in PDAL (PDAL Contributors, 2018).
We obtain epoch-wise changes and associated uncertainties using the multiscale model-to-model cloud comparison (Lague et al., 2013) combined with error propagation (M3C2-EP) as presented by Winiwarter et al. (2021) for a subset of the points of the null epoch, so-called "core points". The subset was created by resampling from the point cloud to obtain a homogenous point density of about 0.8 points/m 2 . For each core point, a normal vector was estimated from a neighbourhood defined by a 5 m radius, and the M3C2 search cylinder had a radius of 0.5 m and a maximum length of 3 m.
The result of this preprocessing is a point cloud containing the core points, each attributed with change values and uncertainties for every bitemporal comparison (i.e., 2 × 70 = 140 attributes). The core point cloud, coloured by bitemporal changes between the first and the last epoch is displayed in Figure 1.

METHODS
The method and the workflow of the study are presented in the following subsections and are visualized in Figure 2. Accordingly, we

Change analysis using Kalman filtering
To extract the erosion form from the dataset, we follow the method presented by Winiwarter et al. (2022a). We first apply Kalman filtering (Kalman, 1960) with Rauch-Tung-Striebel (RTS) smoothing (Rauch et al., 1965) to the time series of change values at each core point location. The method is based on a mathematical description of the surface change system. In the employed implementation, a constant acceleration, nonconstant velocity and change values are used as the state of the system. The Kalman filter and the RTS smoother use the observations, in our case the epoch-wise bitemporal change quantification, in combination with the system description and assumptions on how much the -otherwise fixed -acceleration can change over time, to arrive at smoothed estimates for change value, velocity and acceleration. In addition to smoothing, it allows interpolation over data gaps and further provides uncertainty measures for the derived values.
As the state vector is transitioned from one point in time to the next, the effects of the estimated acceleration are applied to the velocity value, and from there to the change value. Measurements are then used to update these estimates, taking their individual uncertainties into account. In this way, the Kalman filter finds the best estimates for the state given all preceding measurements. The RTS smoother is subsequently applied in reverse order, starting with the most recent measurement. This enables the reduction of uncertainty and incorporation of all of the available data into the final estimates (cf. Figure 3). The thereby obtained uncertainty increases on the last day of the observation period since no future measurements are available to be included in the estimation. We, therefore, disregard the data of the last day in the following analyses.

Clustering and extraction of change characteristics
We use the smoothed time series of change values in a subsequent clustering method based on Kuschnerus et al. (2021). Therein, kMeans clustering is applied to find core points where the change histories over time are similar and then group these core points. From the clustering result, we can identify the erosion form and extract the time series for the affected core points.
The clustering allows for identifying local change occurrences, which can then be segmented from the dataset as individual change forms. In our case, we manually identify three erosion rills (shown by the dashed polygon in Figure 1), which we transfer to a virtual target scene (Section 3.3). The Kalman smoother estimates of change values for the core points within this area are selected and used for the virtual scene generation.

Virtual scene generation
To approach the transfer of a selected and extracted change form to a different geographic setting, we generate a virtual scene of a V-shaped valley of which we scan one side, representing a typical high-mountain slope. We model this slope as a plane, acting as a canvas for us to transfer the change object onto. The change form is then transferred as follows: An eigenvalue decomposition is applied on the covariance matrix of the point cloud subset with extracted changes. Using the obtained eigenvectors, a principal axis transformation is carried out. The point cloud is then meshed using Delaunay triangulation in the two major dimensions. The minor dimension (normal to the best-fitting plane of the dataset) is disregarded, and replaced by the target terrain. In our application, the target terrain is flat, therefore all values in this dimension are initially set to zero.
For each epoch to be transferred, the results of the Kalman smoother are evaluated. As noise has been largely eliminated from these time series, the estimates for displacement are applied directly to the third dimension of the flattened dataset. In the case of a more complex target terrain, they could be added in the direction of the target normal vectors, as shown in Equation 1, where x1, x2, x3 are the three coordinate components of the target mesh nodes, the index o represents the origin data set coordinates, n is the target normal vector in its three components and d is the displacement value obtained from the Kalman smoother.
In case that the displacement in the first two dimensions is not zero, the triangulation should be carried out after adding the displacement and for each epoch separately, to avoid selfintersections or overlapping triangles.
For the last step in the scene generation, the mesh with its displaced nodes is mapped back to the terrain, in our case it is rotated to form the flank of the V-shaped valley.

Point cloud acquisition using VLS
We perform VLS using the open-source software HELIOS++. It allows the simulation of laser scanning acquisitions from different platforms, including static terrestrial and airborne sensors. A key feature is also the modelling of diverging laser beams in a sub-ray approach, where multiple sub-rays are employed to create a single return waveform (Winiwarter et al., 2022b).
In our use case of a heavy rain event causing erosion on steep slopes in a larger area, we are interested in deriving answers to two main questions: With a given laser scanning system, how high can we fly to still reliably identify the change form as significant, and at what point in time (e.g., days after the rain event) can this identification be carried out. For this, we perform multiple simulations with different flying altitudes over the virtual scenes generated for each epoch, resulting in a time series of 3D point clouds for each acquisition strategy.
As a sample sensor, we use a RIEGL VQ-780i. The sensor properties are taken from the sensor datasheet 1 and used as input to HELIOS++ (    Table 2. Scan rates to produce a uniform point pattern for different settings of altitude above ground level (Alt.), flight speed and pulse repetition rate (PRR). Furthermore, resulting point densities (Pdens.) and laser scanning productivity (Prod.) assuming horizontal terrain and no overlap between adjacent flight strips were calculated for the different settings used in the VLS experiment.
Given the pulse repetition rate (PRR), flight speed, and altitude above ground, the scan rate can be adapted to ensure a uniform point spacing in along-and across-track directions on the ground. The scan rate characterizes the number of scan lines per second and is equal to the rotation speed of the polygonal mirror multiplied by the number of mirror facets. Using RIEGL's Ri-PARAMETER 2 software, we derived scan rates, resulting point densities (at nadir position, assuming horizontal ground) as well as the productivity (area covered per time unit, not considering strip overlap) for different combinations of flight altitude, flight speed and PRR. The PRR was adapted to fit the maximum measurement distance according to the datasheet, the flight speed was set to be a minimum for a typical fixed-wing aircraft at the lowest altitude, increasing up to an altitude of 1500 m, and then remaining constant at 75 m/s (145 knots, 270 km/h) as typical cruising speed for such an aircraft. The resulting acquisition settings used in the simulations are shown in Table 2.

Analysis of the simulated point cloud
The simulated point clouds resulting from each parameter set are analyzed for changes using M3C2 distances to quantify erosion. The first epoch is taken as a reference epoch for the respective analysis of each successive epoch. We apply a sensitivity study to find a set of M3C2 parameters that works best for all investigated flying altitudes, resulting in a 3 m radius for normal vector calculation and 0.5 m radius for the M3C2 cylinder projection. From the distance values obtained for points outside of the change form, which are representing stable surfaces, we can derive the registration error, which we quantify to be 0.01 m.
We select the mesh nodes as core points, as they correspond to the core points employed in the Kalman smoother. Therefore, we can directly compare the attribution of significant change from the bitemporal ALS point cloud differences to the changes obtained from multitemporal Kalman filtering. As the latter acts as a model for the simulation, it serves as a reference in our evaluation.
For each parameter set identified by its altitude, we extract (a) the point in time where the change form becomes visible in the change quantification, (b) the point in time when significant change shows the erosion rills, and (c) measures of precision (number of true positives divided by number of positives in the estimation), recall (number of true positives divided by number of positives in the reference), and F1-score (harmonic mean of precision and recall) to investigate the detection rates of locations where change was found to be significant.

RESULTS
We model the surface displacement for three erosion rills extracted from the Vals dataset. The rills show negative displacement values at the top, where sediment is transported by gravitational mass movement, and positive values towards the bottom, where the material is deposited. A sample time series for one of the locations experiencing erosion is visualised in Figure 3. Individual M3C2 distances are shown (red crosses) along with the Kalman smoother estimate (green line) and its respective uncertainty (green shaded area, converted to a 95% significance level) stemming from the original TLS measurements. The uncertainty of the bitemporal comparison is additionally indicated (red dots).
After extracting the cluster, the core points are mapped onto a planar triangle mesh. As an example, the mesh with applied displacement values after five days (at the end of the simulated timespan) is presented in Figure 4. The displacement values in the visualisation are multiplied by a factor of 20 so that the erosion rills and the deposition areas can be identified. The different flight altitudes result in different point patterns on the ground. The scan rate is adapted for each flying altitude to achieve even point densities in along-strip and across-strip directions according to the parameter settings in Table 2. For each flight altitude, M3C2 distances are calculated for each epoch, using the first epoch as a reference. As we select the original core point locations for M3C2 distance calculation, they are the same for each altitude.
Figure 5(a) shows the quantified M3C2 distances for an altitude of 500 m and the surface topography after five days. In Figure 5(b), the core points with significant change are visualized in red. For this altitude and at this point in time, the pattern of the change form is clearly identifiable. In contrast, Figure 6 depicts the same results for an altitude of 3000 m. The form is not identifiable in the M3C2 differences, as it is superimposed by noise. In the significance plot (Figure 6b), it cannot be identified.
Repeating this analysis for all altitudes and intervals of 6 h allows for identifying an optimal acquisition altitude. The precision, recall, and F1-score of the core points for which change was detected are visualized in  the values of precision, recall, and F1-score are summarized in Table 3. While the recall and F1 scores are highest for the lowest flight altitude, the precision value is highest for an altitude of 1000 m at this point in time.
Finally, we map the core points by the altitude where they are first detected, starting from the highest altitude ( Figure 10). The key idea behind this is that an operator may decide which features need to be detectable, for example, based on the real data that was used in the scene generation. By identifying the colour in which these features are visible, the required survey altitude can be selected. In Figure 10, this information is shown for two different points in time, namely 2020-08-23 00:00 (after the main change event) and 2020-08-25 00:00 (at the end of the time series). The outer two rills (I and III) appear in orange and brown colours at the end of the observation period ( Figure 10b) and are therefore detectable from altitudes below  Table 3. Precision, recall, and F1-score of identified significant change for different flight altitudes, 66 h into the time series (2020-08-23 00:00, just after the main change event). 2500 m, whereas the middle rill (II) only starts to exhibit an identifiable linear shape when including the simulation at 1000 or even 500 m altitude. For the earlier point in time, the outer rills become visible at 2000 m, and the middle one is barely visible even when acquired at 500 m altitude.

DISCUSSION
The time series in Figures 7-9 show how the values of precision, recall, and F1-score of surface change detection using VLS data develop over time for each altitude. As expected, all three metrics improve (i.e., achieve higher values) for lower altitudes, corresponding to better detectability. However, the lowest altitude corresponds to the lowest productivity, that is, the area that can be covered within the time of a single survey flight. Therefore, a small reduction in the scores, e.g., by using 1000 m instead of 500 m altitude, may be an appropriate practical consideration. This still allows for a similar detection performance (the difference in F1-score is less than 0.05), but a three-fold increase in productivity (cf. Table 2).
A drop in performance between 2020-08-21 and 2020-08-22 can be seen especially for lower altitudes (Figure 9). Our investigation shows that this is mainly due to a decrease in the precision of change detection at this time (Fig. 7). The precision is given as the ratio of the number of correctly detected change to the two epochs shown in Figure 10. and the total number of points that were identified as changes.
In our case, the number of points where the change (in the reference) is significant at that point in time is quite small -no more than 60 to 70 points until the major change event occurs on the evening of 2020-08-22. A small number of false positives (i.e., points detected as 'change' where no change was recorded in the reference) therefore leads to a seemingly large decrease in the performance. Similarly, the recall value drops due to a number of false negatives, which result from the generally increased Level of Detection for the bitemporal comparison, in contrast to the Kalman smoother result. After 2020-08-23, the evaluation becomes more stable as the number of points attributed to change increases.
Our results confirm the hypothesis that small-magnitude change can only be detected from relatively low altitudes. This has also been reported in the literature for real data, e.g., by Yurtseven (2019), who found that the altitude of an image-based stereophotogrammetric UAV survey has a large influence on errors in the acquired data. The vertical error scales with the altitude in photogrammetric surveys, however, this is not directly applicable for our use case of laser scanning. Here, the effect of the increased beam footprint, especially for flat incidence angles (cf. Fey and Wichmann, 2017), and the decreased point density play a major role. Julge et al. (2014) investigated a multitemporal ALS and TLS dataset and found that they were unable to identify a certain scarp in a single ALS dataset that was acquired at a higher altitude than previous acquisitions. In a different application, Morsdorf et al. (2008) showed that estimates of tree height, fractional canopy cover and leaf area index can be derived more accurately from lower flight altitudes. Also for these cases, simulation could provide a suitable estimate of the data quality that can be achieved, assuming an appropriate scene model representing the vegetation can be found.
In our application, the creation of dynamic scenes was allowed by employing the full TLS time series. To answer the question of the optimal flight altitude, the point in time -and with that, the magnitude of change -has to be considered for data acquisition. For example, for detection on 2020-08-24, the altitudes of 500 m, 750 m, and 1000 m perform similarly well and a choice of 1000 m would be adequate to maximise coverage. However, when change is to be detected between 2020-08-21 and 2020-08-22, a flight at 1000 m only achieves an F1-score of ≈0.45, while a maximum of ≈0.62 can be achieved for flights at 500 m and 750 m, respectively. At this point in time, a choice of 750 m flight altitude would be appropriate. Additionally, the time series allowed us to reduce noise in the input data, which would otherwise have led to a larger Level of Detection for the reference data.

CONCLUSIONS
Laser scanning simulations allow to quickly assess acquisition plans for their fitness for purpose. In the presented example of a topographic time series, stakeholders may select change forms in a high-resolution dataset, e.g., obtained from terrestrial laser scanning, and transfer them to a different area of interest. Using a real dataset as a template has compelling advantages: a more detailed physical description of the underlying geomorphic process or the surface effects of the originally observed change is not required. However, measurement and alignment errors limit the usability of single high-resolution datasets. By applying spatial and temporal smoothing using M3C2-EP and a Kalman filter, we obtain change data with minimised noise levels. From these data, singular snapshots can easily be extracted and used in simulations of change detection using other acquisition strategies in other geographic settings.
In contrast to a pure estimation of point density for different flying altitudes, our approach enables accurate representations of the change forms and magnitudes that are interpretable for each scenario. Furthermore, the laser scanning simulation allows considering occlusion effects unique to each platform and each change form. These occlusion effects can further change over time if the surface geometry is subject to change (e.g., Williams et al., 2021). The increased ranging error resulting from larger footprints and oblique incidence angles is also represented in the simulated point cloud due to the sub-ray approach in HELIOS++. Figure 10. Points detected as significant for two epochs (a and b), coloured by the maximum altitude where they were detected.
For example, green points were detected as significant change for altitudes of 1500 m and below, but not for higher altitudes. Note that this plot is turned counterclockwise by 90 degrees with respect to the other plots in this paper to allow for a larger display area.
Though an analysis of real data, e.g. UAV-based acquisitions, would allow for an even better estimation of detectability at different altitudes, such data is typically not available as dense time series. Reasons include the high cost and logistical as well as legal constraints with permanent monitoring using airborne sensors. Even if a dense time series is recorded, this is typically done at a single altitude. In contrast, we simulated point clouds acquired from 9 different altitudes.
The presented method may, e.g., be utilised in disaster response settings, where a trade-off between data quality and spatial coverage has to be considered in post-event acquisitions. For such cases, our approach provides the required insights to quickly make informed decisions, and avoid unnecessary costs for data acquisitions at less suitable parameter settings. The use case of the erosion and accumulation patterns showed that the choice of flight altitude should not only be related to targeted spatial resolution and production rate, but it rather also requires consideration of the surface change magnitudes and change forms.
Finally, our approach can be extended to include other methods of 4D change analysis. We could then consider different methods for deriving change in the selection of an optimal survey configuration. As such, the presented method can also be used to compare the detectability of change phenomena across different methods, for different altitudes or other survey parameter settings.