Di ﬀ erent Paths for Developing Terrestrial LiDAR Data for Comparative Analyses of Topographic Surface Changes

: High resolution terrestrial laser scanning data (TLS; terrestrial LiDAR) provide an excellent background for quantitative resource estimation through the comparative analysis of topographic surface changes. However, unlike airborne LiDAR data, which is usually provided as classiﬁed and contains a class of ground points, raw TLS data include all of the points of the scanned space within the speciﬁed scanner range. In e ﬀ ect, utilizing the latter data to estimate the volume of the resource by the di ﬀ erential analysis of digital elevation models (DEMs) requires the data to be specially prepared, i.e., separating from the point cloud only the data that represent the relevant class. In the case of natural resources, e.g., mineral resources, the class is represented by ground points. This paper presents the results that were obtained by di ﬀ erential analysis of high resolution DEMs (DEM of di ﬀ erence (DoD) method) using TLS data that were processed both manually (operator noise removal) and with the use of the automatic Cloth Simulation Filter (CSF) algorithm. Three di ﬀ erent time pairs of DoD data were analyzed for a potential gravel-cobble deposit area of 45,444 m 2 , which was located at the bottom of the mouth section of the Scott River in south-east Svalbard. It was found that the applied method of ground point classiﬁcation had very little inﬂuence on the errors in the range of estimating volumetric parameters of the mineral resources and measurement uncertainty. Moreover, it was shown that the point cloud density had an inﬂuence on the CSF ﬁltering e ﬃ ciency and spatial distribution of errors.


Introduction
High resolution LiDAR surveys have been used as basic data for the inventorying of both natural and man-made resources since the 1990s [1]. The obtained point clouds (PCs) are an accurate, unscaled representation of the scanned space [2]. The way they are used and their usefulness depend on the source, the data [3,4] and their resolution [1]. Point clouds that are obtained from airborne laser scanning (ALS) with densities ranging from a few to several dozen points per square meter (pt/m 2 ) and with an accuracy of about 0.3-0.1 m [5] are normally used to inventory huge spaces as well as high-volume resources such as open pit mines, post-extraction heaps [6], railway or road embankments [7][8][9], landslide investigation [10] and various land cover features [11].
On the contrary, data from terrestrial laser scanning (TLS) with densities that range from several dozen to several hundred or even several thousand pt/m 2 are used to inventory larger-scale resources such as quarries, gravel pits, sandpits [12] and areas that do not exceed 1 km 2 [3,4]. High resolution TLS data can be used to estimate the volume of mineral deposits in order to determine the quantitative potential. The estimation of the resources' cubic capacity (e.g., gravel or sand deposits) not only allows for assessing the economic viability of their exploitation or selecting the appropriate type of transport, but also planning the best course of action aimed at the revitalization of post-mining areas.
Repeat surveys are conducted with the intent of controlling for changes at the assumed time intervals in long-term investments as well as natural processes of geomorphic character such as erosion or deposition, which are extended in time. Such data sets, processed into digital elevation models (DEMs), make it possible to get accurate estimations of both changes in the altitude of the topographic surface and the volume in the assumed surface area where estimating the latter occurs indirectly. One of the methods of DEM differentiation at various time intervals was developed by Wheaton [13], which is the DEM of difference (DoD). Comparative studies using the DoD were used for the volumetric estimation of erosion and deposition as well as changes in the position of the topographic surface, e.g., [14][15][16][17][18].
In this method, the accuracy of estimating changes in the surface and thus the loss or increase in the resource depends on the accuracy of extracting from the raw point cloud only those points that represent the ground point [19]. On the one hand, generally available ALS data are usually provided as classified in accordance with the American Society for Photogrammetry and Remote Sensing (ASPRS) Standard LIDAR Point Classes [20] among which ground points are separated into class 2: Ground [20]. On the other hand, TLS data are usually obtained as raw and unclassified. It is in the case of the aforementioned data that the extraction of ground points from the entire cloud still poses a challenge and remains a process with a high risk of error [21]. The highest accuracy of classification is achieved by manual point cloud cleaning of all off-ground objects. However, this is an extremely time-consuming and tedious process, which can account for up to 80% of the time in the DEM development process [22,23]. Despite the widespread use of LiDAR data for DEM creation and the development of software that supports mass point cloud (PC) data processing, the classification process has proved to be extremely difficult to automate [24]. Large data sets of areas with varied land cover have proven to be a particular problem. Recently published works provide information on new generally available algorithms classifying PCs, which enable the separation of ground points, e.g., developed for the American Forestry Service (Fusion) [25], the Simple Morphological Filter [26], the Multilevel Adaptive Filter (MAF) [27] and the Cloth Simulation Filter (CSF) [28]. For noise removal from high resolution digital surface models (DSMs), deep learning methods are also applied, e.g., to sensitive filtering of raw surface elevation data or to effectively increase the automation of the segmentation process [29,30]. The CSF also handles the classification of TLS data relatively well [31]. Consequently, the accuracy of the quantitative estimation of the resource based on LiDAR data depends on the uncertainty range, which, in turn, is connected to survey errors and classification errors. The first group (survey errors) includes the range of accuracy of the measuring devices (scanner rangefinder) and positioning devices (GPS/GNSS, total station, etc.) as well as the accuracy of point cloud integration and georeferencing. While the first group includes among its input parameters the DoD for estimating the range of the measurement uncertainty, it does usually omit the accuracy of the ground point classification.
The aim and novelty of this study is therefore to answer the question as to what extent the ground point classification method affects the estimation of the quantitative parameters of mineral resources. For this purpose, an analysis was carried out on three collections of PC data sets that were obtained by utilizing one type of laser scanner at different time periods in a single area (area of interest (AoI)) whose surface was subject to significant changes during the survey's three-year period. The classification of the ground points in this data was performed manually whilst using the CSF algorithm to obtain three sets of data to compare the range of measurement uncertainty.

Study Area
The study was carried out in the mouth section of the gravel-cobble Scott River valley floor located in the north-west part of Wedel Jarlsberg Land (south-west Svalbard, 80 km apart from the capital Longyearbyen) in the northern part of Sør-Spitsbergen National Park where 65% of the area is covered by glaciers. The Scottbreen valley glacier that occupies about 40% of the catchment area (10 km 2 ) has been retreating at a rate of about 20 m/year and it is the main source of power (90%) of the proglacial river, which flows into the Recherche Fjord (Bellsundfjorden) ( Figure 1A). The study was carried out on an alluvial fan filled with gravel and cobbles, which extends below the gorge section of the river through a raised marine terrace (25-30 m above sea level (a.s.l.)). In the gorge, the river floor narrows to 50 m and again expands to 500 m within the alluvial fan ( Figure 1B).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 15 capital Longyearbyen) in the northern part of Sør-Spitsbergen National Park where 65% of the area is covered by glaciers. The Scottbreen valley glacier that occupies about 40% of the catchment area (10 km 2 ) has been retreating at a rate of about 20 m/year and it is the main source of power (90%) of the proglacial river, which flows into the Recherche Fjord (Bellsundfjorden) ( Figure 1A). The study was carried out on an alluvial fan filled with gravel and cobbles, which extends below the gorge section of the river through a raised marine terrace (25-30 m above sea level (a.s.l.)). In the gorge, the river floor narrows to 50 m and again expands to 500 m within the alluvial fan ( Figure 1B). At the bottom of the narrow gorge section, the Scott River's channels converge into one, often cutting the steep edges of the slopes of tertiary sandstone outcrops with coal, sandstone, slate and sludge inserts. Below the gorge, the Scott River again branches out widely, occupying the entire surface of the alluvial fan whose bottoms are dominated by Quaternary sea gravels [32] with a large At the bottom of the narrow gorge section, the Scott River's channels converge into one, often cutting the steep edges of the slopes of tertiary sandstone outcrops with coal, sandstone, slate and sludge inserts. Below the gorge, the Scott River again branches out widely, occupying the entire surface of the alluvial fan whose bottoms are dominated by Quaternary sea gravels [32] with a large share of cobbles and boulders, which sometimes reach over 1 m in diameter. In the vicinity of the storm rampart separating the alluvial cone from the fjord, there are numerous driftwood trunks on the surface, which become embedded and partly buried in the bottom sediments during storms.

High Resolution Surveying: LiDAR-Derived Digital Terrain Model (DTM)
The survey used data from three measurement campaigns in July 2010, July 2013 and August 2013. High-precision TLS surveys were performed using a Leica Scan Station C10 terrestrial laser scanner. It is a scanner with a range up to 300 m emitting green laser pulses (532 nm) at 50,000 points per second [pt/s] [33]. During the 2010 and 2013 campaigns, various referencing techniques were used [34]. In 2010, a scan of the test area was performed in a local reference system. The analyzed section of the bottom was scanned from eleven sites that were registered by nine ground control points (GCPs) with an accuracy of ± 0.009 m. After the integration of the entire model by Leica's Cyclone 8.1 software (Leica Geosystems AG, Heerbrugg, Switzerland), it was georeferenced basing it on position five of the proportionally distributed GCPs in the scanned area (measured by a Leica 500 dGPS with an accuracy of ±0.02 m). In 2013, a network of eleven GCPs was first established to determine georeferencing. The instrumentation was based on the use of dual-frequency GNSS TOPCON receivers, namely TopCon Hiper II, which were operated by employing the base-rover survey method. An application of the real-time kinematic Global Navigation Satellite System (rtkGNSS) (GPS+GLONASS; both a Global Positioning System (GPS) and its Russian equivalent named Globalnaya Navigacionaya Sputnikovaya Sistema (GLONASS)) allowed us to make measurements with an accuracy of ±0.02 m. The GCPs network was then used to locate both the scanner and the reference points. In 2010, the majority of the scans were made in high resolution mode (0.05/100 m), while for some part of July and the ensuing month of August in 2013, all measurements were obtained in medium resolution (0.1 m/100 m).

Classification of Ground Points and DEM Development
The measurement data from each campaign were integrated into the DTM by Leica's Cyclone 8.1 software ( Figure 2).
The point clouds in each of the three DTMs were then classified in order to separate the ground points. The classification was performed by adopting two methods; one done manually by the author and the other done automatically by using the CSF algorithm [28].
The 'manual' ground point classification was carried out in Leica Cyclone 8.1 software by removing all off-ground objects from the point cloud ( Figure 2). As the High Arctic does not have any high vegetation (only dwarf vegetation and there is little of it in the valley floor), off-ground objects were limited to the following: living organisms (humans, birds, insects), objects abandoned in the coastal area during storms (fishing buoys, driftwood logs, fragments of a wooden winch), measurement installations (traps, piezometers, meteorological stations) and 'noise', i.e., points representing non-existent objects related to the misinterpretation of the sun's rays directly striking the mirror of the scanner at an angle of about 45 degrees and the refraction of laser beams at the edge of water ( Figure 3). The 'limited fence' function was used to remove off-ground points from the 'view' module, which made it possible to work on a fragment of the cloud that had been selected by the operator (Figure 3). This function significantly reduced the unintentional elimination of points from the cloud by the operator. At the same time, raw DTMs were filtered with the CSF algorithm ( Figure 2) as a plug-in for CloudCompare 2.11. The following advanced parameter was used: the slope processing was set to 'active' (the valley section varies in altitude); then (1) the Cloth resolution was set to 0.1 m (the value refers to the grid size of cloth that was used to cover the terrain); (2) the max iterations were set to 1000 (the value refers to the maximum iteration counts of terrain simulation) and (3) the classification threshold was set to 0.1 m (the value of the threshold to classify the point clouds into ground as well as non-ground parts, which was based on the distances between points and the simulated terrain). In settings (1) and (3), respectively, the lower of the available units was used that was then the same as the unit in the point clouds processing. Lastly, using LP360 7.0 software, six ground point DTMs were converted into triangulated irregular network (TIN) models and further rasterized as a GeoTIFF with 0.1 m resolution to produce three pairs of DEMs ( Figure 2). Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 15    [28], Geomorphic Change Detection (GCD) diagram http://gcd.riverscapes.xyz, [13,15].

Surface Change Detection by DEM of Difference (DoD)
For the purpose of obtaining information about the influence that the classification method has on the range of uncertainty in the estimation of surface and volume differences, three pairs of DEMs (prepared in exactly the same manner as mentioned above) were analyzed in Geomorphic Change Detection 7.4.4 (GCD) plug-in version for ArcGIS 10.8 software ( Figure 2). According to the guidelines, as stated by the author of the DoD method [13], the approach, which is founded on the error assessment based on the spatial variable, was applied in order to calibrate the calculations.
There was an analysis done on data from the same source but differing in the method of separating the point clouds so the accuracy of the DEMs was closely related to the effectiveness of the ground point classification. The manually-classified DEM was assumed to be basal (new) and so least affected by classification errors. The DEM that resulted from Cloth Simulation Filtering was considered to be comparative (old) and more exposed to erroneous points, which had been omitted by the algorithm in Figure 2. Inherent to the calculation of the surface error for both DEMs was the assumption of the minimum detection level (minLoD) at ±0.1 m. A 'mask' for the AoI was applied, thus limiting the calculation results only to the gravel-cobbles surface of the 45,444 m 2 alluvial fan. The final stage included applying 'error propagation' to determine the level of measurement uncertainty [15,33].

Results
The implementation of the various scanning methods in subsequent measurement campaigns as well as the application of two predefined measurement resolutions resulted in both the number of points and the average point density in the cloud being different. Both 'manual' and 'automatic' classification processes led to a reduction in the number of points within the range from 0.01% to 2.7% ( Table 1). Out of the eleven individual point clouds from 2010, nine were made in high resolution. However, a relatively small number of scan stations produced a total of 52,358,443 points in the DTM, while the average density of points for the AoI equaled 345 pt/m 2 . Moreover, both of these parameters were similar to those effectively recorded in the third measurement campaign (August 2013) when the average resolution was predefined for all eleven scan stations. The result was a total number of 45,207,502 points in the DTM and an average point density of 309 pt/m 2 for the AoI. The DTM from July 2013 was a fragment of the model of the whole valley so in addition to the eleven scan stations located at the bottom of the alluvial fan, the DTM supplemented points by taking them from seven other locations on the slopes as well as the lifted marine terraces plateau. Therefore, this last DTM was characterized by a much larger total number of points equal to 70,071,395 as well as a higher average point density for the AoI at 437 pt/m 2 (Table 1). The analysis of surface differences carried out with the use of GCD software [12,14] made it possible to estimate the physical parameters of the gravel-cobble deposit, which was the result of errors in the CSF algorithm classification. The results of the DEMs comparison from 2010 show negligible differences between both methods of classifying the ground points ( Figure 4;    Surfaces that were misinterpreted by the CSF represented only 3% of the AoI. Insofar as lowered areas slightly predominated within the range of surfaces in the raw point cloud, the lowered areas at threshold values (minLoD = 0.1 m) constituted as much as 79% of them ( Table 2). The height change diagram shows that the average total thickness of difference was 0.23 m (± 0.10) with an average net thickness difference of −0.16 m (± 0.08) (Figure 3). All error areas were concentrated in the north-east part of the AoI and occurred within or at the boundaries of water surfaces (Figure 4). The oval shape of these areas indicated that the CSF algorithm omitted small concentrated groups or even single points that were found below the water level.
In relation to the data from July 2013, the area of detected changes was identical to the AoI area for raw point clouds whereas it was only 0.2% of the AoI area for the threshold values (minLoD = 0.1 m) ( Table 3). Among the areas misinterpreted by the CSF, the raw PC was slightly dominated by elevated areas (52%) while for threshold values it was the lowered areas that prevailed (95%) ( Table 3). An analysis of the height change diagram showed that the averages of the total as well as the net thickness of difference were smaller than in 2010 and were from 0.16 m (± 0.10) to −0.15 m (± 0.09), respectively ( Figure 4). As in 2010, all of the error areas were concentrated in the north-east of the AoI where they occupied similar locations and took on comparable shapes ( Figure 5).
In terms of the data from August 2013, the area of detected changes was less than 60% of the AoI area for raw point clouds whereas it was 0.1% of the AoI area for the threshold values (minLoD = 0.1 m), which was similar to the observations made in July (Table 4). Among the areas misinterpreted by the CSF, the raw PC was slightly dominated by elevated areas (52%) while for the threshold values, it was the lowered areas that prevailed (95%) ( Table 4). An analysis of the height change diagram showed that the averages of the total and the net thickness of difference were similar to 2010 and were equal to 0.21 m (± 0.10) and −0.17 m (± 0.09), respectively ( Figure 5). All of the error areas were concentrated on the right bank of the river in a section that was parallel to the storm dike. Lowered areas dominated there. Meanwhile, the largest area and height differences occurred coincidentally at the mouth of the river that flows into the fjord (Figure 6). Table 3. Changes in physical parameters of the Area of Interest in July 2013 calculated volumetrically with reference to the total volume of the gravel-cobble material's net change recorded by the DoD (both lowering and raising) from the manually cleaned-up (new) to the CSF-filtered (old) DTM. The uncertainty analysis was calculated by simply using the minimum level of detection (minLoD) approach while accounting for a uniform error of 0.1 m.

Discussion
In areas where mineral resources are extracted, the key issue is to quickly determine the quantities of extracted or exploitable mineral resources. An accurate estimation of deposit parameters requires taking into account three potential sources of uncertainty: (i) the physical properties of the material in the deposit, (ii) the measurement methodology along with the accuracy of the measuring equipment and also the measurement methods in the case of modern remote sensing (iii) the methods for classifying ground points as well as the interpolation of the ground surface and their respective analysis. In this context, a significant part of the research is focused on estimating the cubature or bulk density of the raw material under in situ conditions, all of the while using conventional and modern remote sensing methods such as digital close-range photogrammetry (structure from motion (SfM) photogrammetry) and terrestrial laser scanning. These works usually document relatively small (about 5%) differences between the parameters determined in a laboratory and those under

Discussion
In areas where mineral resources are extracted, the key issue is to quickly determine the quantities of extracted or exploitable mineral resources. An accurate estimation of deposit parameters requires taking into account three potential sources of uncertainty: (i) the physical properties of the material in the deposit, (ii) the measurement methodology along with the accuracy of the measuring equipment and also the measurement methods in the case of modern remote sensing (iii) the methods for classifying ground points as well as the interpolation of the ground surface and their respective analysis. In this context, a significant part of the research is focused on estimating the cubature or bulk density of the raw material under in situ conditions, all of the while using conventional and modern remote sensing methods such as digital close-range photogrammetry (structure from motion (SfM) photogrammetry) and terrestrial laser scanning. These works usually document relatively small (about 5%) differences between the parameters determined in a laboratory and those under field conditions by TLS [12]. This is all the more so because during the extraction or movement of the sediment resulting from natural processes (landslides, mud and rubble run-off, debris flow, etc.), the cohesion and loosening of the sediment structure is reduced. As a result, it is often found that the volume of the deposited material is larger than that of the erosion niche. Bremer and Sass [35] highlighted an average 7% prevalence of deposited volumes over the detected eroded volumes for two debris flow events in the Halltal (Austrian Alps). Kociuba [36] also pointed to a 6% growth of volume material that was deposited in a small alluvial fan that developed directly below a rill and suggested that this was due to the loosening of the transported material as well as multiple stages of aggradation (multilayer deposition on the cone). In fact, the loosening coefficient for the heterogeneous mineral resources of the perlite deposits was estimated to be 1.38 [12].
A similar uncertainty range (about 2%) is associated with the use of various algorithms for filtering ground points. Studies aimed at the accuracy of the algorithms most often focus on the following: comparing the effectiveness between the algorithms or deviations from differential morphological profiles and reference surfaces [37], discarding returns that exceed a threshold curvature calculated from a surface that is interpolated using a thin plated spline [38], identifying in the meantime the lowest elevation point in every cell that is defined by an operator and, finally, creating a surface by interpolating these lowest points [39].
The key classification problems are presented in many studies on the effectiveness of algorithms that classify ALS data; however, none of these works involve data obtained through terrestrial scanning [23,27,28,37,40,41]. Although some researchers think that the most effective classification of TLS point clouds is performed by manual point removal [42], others consider that automation of the classification process, especially with algorithms that do not use the altitude threshold parameter, significantly increases the time efficiency and eliminates errors resulting from operator bias [21].
Many algorithms that successfully classify ground points from ALS data are much less effective in separating them from TLS clouds. This results, among others, from the following: (i) a different angle of attack of the laser beam, (ii) a significantly higher number of points in the cloud, (iii) a different density of points in the cloud (density decreases with the distance from the scanner) but also (iv) a lack of appropriate recommendations for the selection of classification parameters [24]. Available as a plug-in for Cloud Compare software, the CSF algorithm [28] is one of the few commonly available algorithms that has a graphical interface with tooltips or hints, which are shown in both graphic and text form, in order to help the user understand how the algorithm works and to significantly accelerate the definition of optimal parameters. It also provides relatively good results in the classification of ground points from TLS data due to its use of two factors; the inversion of the ground surface and coating the ground surface with a simulated cloth of specific strength and thickness [28]. This makes it possible to eliminate the previously mentioned differences between ALS and TLS clouds, which influence the less effective results of other algorithms that contain TLS data. Unfortunately, what is a 'strength' is also a minor 'weakness' of this method. The algorithm is sensitive to errors occurring on the boundary between land and water where we often have to deal with noise. As this study shows, it is below the water surface that the noise is not always filtered out by the CSF. Keeping these points leads to an artificial lowering of the surface. In addition, defining the thickness of the canvas results in the possibility of 'shearing' objects that exceed this parameter; for example, large boulders buried at the bottom of a gravel-cobble bed creek. The results of the study show that the applied method of filtering ground points by the CSF to a small extent increases the scope of uncertainty in estimating the quantitative parameters of the analyzed resource. The study adopted a range of MinLoD = 0.1 m as this value referred to the smallest possible key parameters defined in the CSF-the thickness and flexibility of the cloth. It also referred to the raster cell size used in this study. Working with this limitation caused the values that were estimated as the DoD threshold parameterizing the areas that deviated from the reference area to be only about 1% of the AoI area. The assumption of MinLoD = 0 as well as the application of error propagation made the values for raw data and the estimated DoD threshold area overlay in up to 100% of the AoI area. The obtained values for the raw and threshold data were also the same. However, the graphical representation of the DoD of these changes and the height indicators show (Figure 7A,B) that the uncertainty associated with the accuracy of the algorithm is not a parameter that significantly affects the results of the quantitative assessment of the deposit resources.

Conclusions
This paper answers the question as to what extent the method of classifying ground points affects the result of a mineral resource volume estimation. During the analysis of DTM pairs derived from the same area at three different times, it was shown that the errors associated with the applied automatic algorithm of ground points classification did not differ significantly from the model

Conclusions
This paper answers the question as to what extent the method of classifying ground points affects the result of a mineral resource volume estimation. During the analysis of DTM pairs derived from the same area at three different times, it was shown that the errors associated with the applied automatic algorithm of ground points classification did not differ significantly from the model developed by the operator (manual noise removal). The main accomplishments of this research work are the two methods that were used in the classification of ground points (manual removal and CSF). These generated very similar results for cases of mineral resource exploitation (or potential exploitation) areas such as gravel, perlite and sand and in the analyses of changes to the land's natural surface without vegetation cover such as landslides, river valleys, etc.
The CSF algorithm's automatic classification errors minimally affected the scope of the quantitative estimation of the resource. In the assumed range of MinLoD = 0.1 m, the errors only pertained to about 1% of the AoI area while their average height deviation values were determined not to usually exceed 0.2 m.
The limitations of this study are that the analyses that were conducted in the course of the study have broadened the knowledge base in this topic area although it has only been in relation to a specific type of resource (gravel-cobble deposits) and two classification methods (manual and CSF).
In DoD analyses, any uncertainties that are associated with errors in the interpretation of the ground surface are usually not taken into account. We may question whether this is the right approach. It seems that due to the intensive development of current and newly available classification algorithms, this topic should be further explored and it ought to be looked at in terms of testing both new algorithms and other types of resources along with their conditions of occurrence.
Funding: This research received no external funding.