A Workflow to Estimate Topographic and Volumetric Changes and Errors in Channel Sedimentation after Disturbance

Light Detection and Ranging (LiDAR) methods, such as ground-based Terrestrial Laser Scanning (TLS), have enabled collection of high-resolution point clouds of elevation data to calculate changes in fluvial systems after disturbance, but are often accompanied by uncertainty and errors. This paper reviews and compares TLS analysis methods and develops a workflow to estimate topographic and volumetric changes in channel sedimentation after disturbance. Four analytic methods to estimate topographic and volumetric changes were compared by quantifying the uncertainty in TLS-derived products: Digital Elevation Model (DEM) of difference (DOD), Cloud to Cloud (C2C), Cloud to Mesh (C2M), and Multiple Model to Model Cloud Comparison (M3C2). Mean errors across surfaces within each dataset contributed to a propagation error of 0.015–0.016 m and 0.017–0.018 m for the point clouds and derived DEMs, respectively. The estimated error of the total volumetric change implied increased errors in the conversion of point clouds into a surface by C2M and DOD; whereas C2C and M3C2 were generally simpler, efficient, and accurate techniques for evaluating topographic changes. The comparison of methods to analyze TLS data will contribute to applications of remote sensing of hydro-geomorphic processes in stream channels after disturbance. The workflow presented also aids in estimating uncertainties inherent in data collection and analytic methods for topographic and volumetric change analysis.


Introduction
Understanding the impacts of natural and anthropogenic disturbances on ecosystem processes requires reliable and efficient methods to survey and detect topographic and geomorphologic variability in affected regions.Novel methods to evaluate the impacts of disturbances on stream channels include high-resolution three-dimensional remote sensing.Terrestrial Laser Scanning (TLS) is a ground-based Light Detection and Ranging (LiDAR) remote sensing method that collects high-density point clouds, which can be developed into digital elevation models (DEMs).Repetitive point cloud data collection can then be used to quantify disturbance-related topographic and volumetric changes of stream channel sediment over time and aid in understanding geomorphic responses after disturbance [1][2][3][4][5].
Many methods are available to analyze and compare point cloud data and derivative DEMs [6], which are often accompanied by inherent uncertainty and error [7].These spatially distributed errors occur during data collection and processing [7][8][9][10][11][12], which can be enhanced in fluvial systems with large topographic variations [13] or propagated during analysis of multiple scans [14][15][16] discussed in the next section.These errors have significant implications for data analysis and interpretation [7][8][9][10][11][12].Thus, a need exists to understand the challenges and benefits of various TLS analytic methodologies and to develop a framework to guide optimal processing and analysis of TLS information.This paper: (1) reviews techniques and challenges for processing TLS data for analyzing topographic changes in stream channels in terrains after disturbance; (2) reviews and compares methods for estimating changes between multi-temporal TLS datasets; (3) develops an error analysis workflow for analysis of point clouds and DEMs derived from TLS datasets; and (4) tests the proposed workflow to compare four methods to analyze multi-temporal TLS point clouds collected for a burned stream channel following the 2012 Waldo Canyon Fire (Colorado, USA).

Review of Terrestrial Laser Scanning for Disturbance Applications
Numerous studies have applied ground-based and aerial remote sensing to provide data needed to estimate changes in topography after landscape disturbance [17,18].This review focuses on TLS, a LiDAR method useful in the context of quantifying topographic and volumetric changes in channel sedimentation related to topographic variability after disturbances.

Overview of Methods to Detect Topographic Changes
LiDAR is a remote sensing method [19] increasingly used to acquire aerial or ground-based elevation data points for studies estimating topographic or volumetric changes.A scanner emits a laser pulse, which interacts with various media and objects such as ground surface, vegetation, or buildings.The laser pulse is scattered and reflected back to the scanner after interacting with the land surface or vegetation.The travel time of each light pulse and distances are used to estimate a three-dimensional (x, y, h) "point cloud", where x and y indicate horizontal position and h is the elevation.Airborne LiDAR is capable of sampling large areas of 10-1000 km 2 with a resolution of 0.5 to 100 points/m 2 , whereas ground-based LiDAR is typically capable of scanning an area of 0.1-1 km 2 with a higher resolution of 50 to 10,000 points/m 2 [20][21][22].
Terrestrial Laser Scanning is a ground-based LiDAR method that collects high-density point clouds and generally consists of three main processing steps-data acquisition, global co-registration of the point cloud, and estimation of the deformation [23]."Bare Earth" land surfaces are created by filtering vegetation [24,25].The data are then interpolated to triangulate or grid (rasterize) the point cloud information.Multiple datasets repeatedly collected in the same location over time can be compared to develop estimates of topographic changes.In addition to LiDAR, structure-from-motion (SfM) is a low-cost photogrammetric alternative for collecting high-resolution point clouds and creating DEMs, which can be analyzed using processing techniques analogous to the methods discussed in this paper [26].
Inherent uncertainty in LiDAR-based methods exist and are attributed to a combination of instrument error, point data collection techniques, and processing methods.Field sampling methods, surface composition and surface roughness related errors, systematic error (e.g., accuracy of the survey equipment and registration errors), the density and distribution of data points used to represent the surface, and interpolation methods to derive DEMs from the data can contribute to spatially distributed errors in the dataset [7][8][9][10][11][12].Missing data due to shadowing effects from vegetation or other obstructions can result in areas with sparse to no points [22,27,28], and ensuing errors and uncertainty can be enhanced in regions with high geometric complexity and large topographic variation [13,27].Moreover, the aforementioned errors can propagate from analysis of multiple point cloud datasets from the same areas producing DEMs from different collection times [11,[14][15][16]29].

Synthesis of TLS Methods and Development of an Error Analysis Workflow
To aid in quantifying error, we compared methods appropriate for analysis of LiDAR point cloud data in a disturbed system to develop a four-part workflow for data acquisition, processing, interpretation, and quantification of error (Figure 1).A variety of methods and approaches exist for each part of the workflow, such as those mentioned in Heritage et al. [9] and Fuller et al. [30].This review utilized the Riscan Pro and CloudCompare algorithms and methods, which are summarized in Table 1.As a demonstration, we applied the workflow to test four methods to analyze TLS data collected after a wildfire disturbance in Waldo Canyon, Colorado (Section 4). Figure 1 summarizes the workflow for analyzing the LiDAR data and illustrates a proposed hierarchical processing for analysis.
Remote Sens. 2018, 9, x FOR PEER REVIEW 3 of 19 interpretation, and quantification of error (Figure 1).A variety of methods and approaches exist for each part of the workflow, such as those mentioned in Heritage et al. [9] and Fuller et al. [30].This review utilized the Riscan Pro and CloudCompare algorithms and methods, which are summarized in Table 1.As a demonstration, we applied the workflow to test four methods to analyze TLS data collected after a wildfire disturbance in Waldo Canyon, Colorado (Section 4). Figure 1 summarizes the workflow for analyzing the LiDAR data and illustrates a proposed hierarchical processing for analysis.

LiDAR Workflow
After acquiring the LiDAR data (I; Figure 1), data preparation includes registration and georeferencing, data filtering, removing vegetation, and point cloud triangulation detailed subsequently (II; Figure 1).In general, LiDAR data from adjacent scans are registered for alignment and converted into a single dataset of the channel, which is georeferenced into an Earth-centered Earth-fixed coordinate system.Multiple scans from adjacent locations can be registered based on tie points (reference points), which are fixed control points with known location coordinates [31].
The process of data filtering typically manipulates and filters point clouds to create a homogenous point cloud to decrease the size of data, making the processing and analysis more efficient.The management and analysis of large and complex data is extremely cumbersome and inefficient.Several data filtering methods are available, such as a range gate method to filter points by a range of distances from the scanner, an intensity gate to filter by a range of return strength of the reflected pulse for data points, or a point filter to keep each n-th point (n is a user defined factor) [32].The octree structure is a common technique, which is often applied to filter data [33,34].It is based on an initial cube, sub-divided into eight equally sized smaller cubes.The octree filter segments the point cloud into these cubes, with the data within each cube averaged into a single point, which is the center of gravity for each cube [32].The size of the initial cube and the increment value are determined (usually between 0.1 m and 1 m) and represents a stopping point for sub-division [32].
Removing vegetation to define the bare Earth ground surface utilizes a height difference between two adjacent points [24,25].This method separates points above the bare Earth (such as vegetation) from points on the bare Earth (the ground surface) by analyzing the distances of the points from an estimated ground surface.The ground surface is the lowest elevation point within the point cloud.The filter proceeds from coarse to fine grids.For each cell of the grid, a representative point is selected from the point cloud based on a percentile parameter, which defines the percentage of points within a cell that is below the representative cell point (RCP).For example, in the case of a percentile of 1%, the RCP is selected so that 99% of the points are above the RCP and 1% of the points are below the RCP.Then, those RCPs are used to estimate a surface for each non-empty cell by estimating a plane through the RCPs.For each cell, there is a defined tolerance range above or below the estimated plane.For each cell, points that have not yet been classified as "off-terrain" and are outside the tolerance range are marked as "off-terrain" points and excluded.The procedure continues for the remaining points by subdividing the cells into four equal size cells.Then, the workflow uses a fine filtering process based on a tolerance value smaller than the tolerance range in the previous step.The next step defines a maximum slope angle as the maximum slope of the estimated ground surface.This value is applied to the filtering points, preventing the filtering of steep surfaces [32].The grid size, the percentile parameter, the tolerance range, the fine filtering tolerance value, and the maximum slope angle are defined according to the purpose of filtering (e.g., vegetation removal) [32].
After vegetation removal, interpolation and triangulation of the bare-Earth points create high resolution DEMs.A DEM is a set of elevation values, which are recorded as a regular grid in a square, triangular or rectangular form [35] that is developed through data interpolation such as Delaunay triangulation [36,37] or Kriging [9,30].Triangulation is the process of creating a surface from a point cloud by connecting the data points with triangles.During this procedure, maximum triangle edge length, maximum triangle tilt angle (the maximum angle between the triangle normal and the ray of sight based on data orientation) and minimum triangle angle are defined.For example, the two-dimensional Delaunay triangulation algorithm linearly interpolates data points in a plane such that no point is inside the circumcircle (the circle that passes through the triangle vertices) of any triangle [36,37].It is based on the two-dimensional coordinates of the points mapped onto the computer screen.Triangulation is performed from the current point of view, which is defined depending on the feature of interest (e.g., top view for active channel triangulation).Considering the point of view before analysis is important to avoid overlapping and inappropriate triangulation [32].Once a triangulated surface is created, the DEMs are generated and the data are smoothed to reduce noise in the interpolated topography.Smoothing modifies the surface structure of the DEM by smoothing the edges or decimating (reducing) the number of triangles [32].

Change Estimation and Analysis
Several methods exist to compare point cloud data, including Cloud to Cloud (C2C) [38], Multiple Model to Model Cloud Comparison (M3C2) [6,38], Cloud to Mesh (C2M) [6,38], and comparing gridded data sets to produce a DEM of difference (DOD) [14,15] that quantifies volumetric changes [6].Yet, topographic change detection using TLS data remains challenging because of the point cloud complexity, number of data points, and inconsistency of scan positions that can contribute to minimal point-to-point overlap.
The methods of C2C, M3C2, and C2M allow estimation of topographic changes between two successive scans of the same reaches in a three-dimensional environment.The C2C and M3C2 methods compare two point clouds, whereas C2M compares derived DEMs.These methods highlight areas of loss or accumulation as an elevation (topographic) difference.Positive values indicate areas of deposition that raise elevations, whereas negative values indicate areas of erosion that lower elevations.C2C is the simplest and fastest method to compare point clouds, since it does not need gridding or triangulation of the data [34].For each point of the second point cloud, a closest point is defined in the first point cloud, and the surface change is estimated as the vertical distance between the two points.The technique is suitable for rapid change detection in very dense point clouds.It is not, however, an accurate method to measure horizontal changes, or distance between vertical surfaces [34,38].
The M3C2 algorithm measures the local distance between two point clouds in the direction normal to the surface, operating directly on point clouds without triangulation or gridding [38].Lague et al. [38] performed an error analysis based on estimated vertical displacements by an algorithm versus field-based displacement measurements.Their study concluded that, among these methods, C2M and M3C2 are most accurate for estimating distance differences between two surfaces [38], which is important for quantifying erosion and deposition between two successive scans.The most common method, C2M, evaluates the change in elevation of the surface.The vertical distance between data points within a point cloud (or surface mesh) and a reference three-dimensional mesh defines surface change [38].However, interpolating missing data introduces uncertainties that are difficult to quantify [38].The C2M and M3C2 methods are applicable for change estimation of vertical surfaces such as cliffs or riverbanks, while M3C2 provides accurate horizontal displacements [6].
For volumetric change estimation, the DOD method [11] compares two DEMs (reference and base mesh) and quantifies volumetric changes caused by erosion or deposition.Similar to the methods outlined above, the DOD method highlights areas of sediment loss or accumulation with positive values that indicate deposition and negative values that indicate areas of erosion.In the DOD method, two point clouds are gridded to generate DEMs and for comparison.The two DEMs are then differentiated pixel-by-pixel by estimating a vertical distance [38].For DEM differencing, two meshes can be converted to a raster structure in Riscan Pro or a grid in graphing software such as Golden Software Surfer.Volumetric change is calculated by summing the raster cell's volume, which is the cell area multiplied by its difference in height.This method uses a reference plane on top of the DEM surfaces to project the changes and estimate the difference volume as a quantity of erosion or deposition.Areas of the base mesh, which are below the reference mesh, are added to the erosion volume.The deposition volume includes areas of the base mesh that are above the reference mesh [32].
The application of the DEM differencing identifies and quantifies morphological change in river channels, such as volumes of scour and fill at the reach scale [9,11,14,15].It is a 2.5D environment, where DEMs have an assigned mean elevation at the center of each pixel and a topographic data reduction occurs when point cloud data is reduced from 3D to 2.5D [6].This method introduces grid elevation uncertainty in three-dimensional environments and overhanging areas.It is unable to grid TLS point cloud data for rough surfaces with missing data due to occlusion [38].

Error Analysis
It is important to understand the potential for errors introduced by field sampling methods and surface complexity to interpret the results and correctly differentiate between error sources.This includes distinguishing between topographic changes and systematic errors (e.g., accuracy of the survey equipment and registration uncertainty between point clouds), the composition of data points (e.g., density, sparsity, and spatial unevenness representing the surface), or the interpolation methods to derive DEMs from the source data [7][8][9][10][11][12].Error analyses provide estimates of measurement accuracy (difference between the measurement and the actual value) that influences change detection from repeated scans [10].These analyses include: (1) assessing the accuracy of the point cloud data and DEM surface for each TLS dataset; (2) propagating the uncertainty from each dataset into the difference surface; and (3) determining the propagated uncertainty on volumetric changes measurements [11,14,29].
Calculation of the propagated error in both surfaces utilizes the Root Mean Square Error (RMSE) (Equation ( 1)) to establish a minimum level of detection (LoD min ) [11].The minimum level of detection (LoD min ) is the smallest measurable elevation change between two surveys [13,39].The LoD min determines if the elevation difference between two surfaces is real erosion or deposition, rather than uncertainty in measurement.Typically, the elevation changes above this LoD min are considered as real [12].
The RMSE is: where S is the TLS point clouds and O is the reference data (or triangulated surface).RMSE defines the difference between the TLS and reference data, and zero is the optimum value.It is important to note that errors in the control points can exist; however, following Fisher and Tate [35], the reference data from the field are assumed to have minimal errors.

Error in Points
The error of a set of collected point clouds can be determined by comparison with another set of known, usually more accurate measurements.Such reference data are presumably error free [35,40].The uncertainties are either spatially uniform [39] or vary spatially based on wet and dry areas, vegetated and un-vegetated areas, topographic roughness and complexity, and high point density versus low point density areas [9,14,16].Uncertainties in the raw point cloud data can arise in the direction of x, y, and h.Vertical (h) uncertainties significantly influence DEM quality and DOD measurements on flat areas, while horizontal (x and y) uncertainties are more significant in areas with more roughness.
To account for potential sources of uncertainty, the RMSE (Equation ( 1)) provides an estimation of a spatially distributed uniform error (or average error across the reach).A spatially uniform error estimate assumes that error is not a function of location and is constant in space [16].If independent check points are not surveyed, overlapping data points from a previous survey with the same ground control network and coordinate system, using fiducial (or reference surfaces in areas that have not changed such as bedrock outcrops), can provide an alternative [41,42].The δ h (Equation ( 2)) calculates the propagation of errors of point clouds in x, y, and h directions: where 1 and 2 denote the first and the second temporal scan, respectively.

Error in DEMs
Interpolation algorithms approximate the surface elevation based on the elevation of the surrounding measured data.If measured data are sparse in regions of high topographic variation, the topography contributes larger uncertainty to the actual surface [9].Complex and rough surfaces increase the risk of occlusion in data points, which impacts interpolation during the DEM creation, reducing the accuracy for change detection.Working directly on the point clouds avoids this type of uncertainty [38].However, the DOD is required in estimating volumetric changes.Thus, the error in both DEMs can propagate in this method.
After point cloud triangulation, the potential error in DEMs can be evaluated by comparing the surface before and after triangulation [1].To estimate the total propagated error in a DEM of Difference in the h direction, δ h(DoD) , (Equation ( 3)) combines the estimates of errors in the first DEM (RMSE 1 ) and the second DEM (RMSE 2 ) [11,14,16]: The volumetric error, δ v (Equation ( 4)), determines the uncertainties in the volumetric calculations and is defined using the total propagated error in a DEM of Difference [2,14]: where A v is the area of the active channel.

Background for Test Case Study Area
TLS data acquired following a wildfire disturbance in Colorado, USA [5] provided a test case to compare the analysis methods reviewed in this paper and illustrate the workflow presented in Section 3 (Figure 1 and Table 1).The 2012 Waldo Canyon Fire occurred in Pike National Forest, northwest of Colorado Springs, Colorado in El Paso County [43].This study utilized data collected as part of a larger investigation of the response of mountain streams following the Waldo Canyon Fire, described in Chin et al. [44]; Kinoshita et al. [45]; and Chin et al. [5].The test presented here focused on a channel reach in the downstream portion of the Williams watershed (Figure 2), with a drainage area of 4.87 km 2 .This study reach, henceforth referred to as "Williams Downstream," has a length of 165 m and encompasses "Willis Reach" described in Chin et al. [5].Ground-based laser scanning of this study reach provided the data for a comparison of each change estimation method and errors.We demonstrate the steps outlined for the workflow in order: (1) data acquisition, (2) preparation, (3) change estimation and analysis, and (4) error analysis.
Remote Sens. 2018, 9, x FOR PEER REVIEW 8 of 19 Section 3 (Figure 1 and Table 1).The 2012 Waldo Canyon Fire occurred in Pike National Forest, northwest of Colorado Springs, Colorado in El Paso County [43].This study utilized data collected as part of a larger investigation of the response of mountain streams following the Waldo Canyon Fire, described in Chin et al. [44]; Kinoshita et al. [45]; and Chin et al. [5].The test presented here focused on a channel reach in the downstream portion of the Williams watershed (Figure 2), with a drainage area of 4.87 km 2 .This study reach, henceforth referred to as "Williams Downstream," has a length of 165 m and encompasses "Willis Reach" described in Chin et al. [5].Ground-based laser scanning of this study reach provided the data for a comparison of each change estimation method and errors.We demonstrate the steps outlined for the workflow in order: (1) data acquisition, (2) preparation, (3) change estimation and analysis, and (4) error analysis.

TLS Data Acquisition
UNAVCO collected TLS data from Williams Downstream reach using a Riegl VZ-400 scanner with a Nikon D700 20 mm camera mounted on a tripod.The scanner system operated in a near-infrared wavelength (800-1000 nm) and the maximum scanner range of data collection varied from 60 to 600 m, depending on object brightness.This setup allowed a 360-degree panorama scan in the horizontal direction and 100 degree in the vertical direction, which measured at a rate of up to 122,000 measurements/sec.The data coverage depended on the incident angle of the laser, as determined by the height of the scanner above the ground and the clarity of the site with respect to the scanner (obliqueness) [46].The minimum horizontal and vertical step size is 0.0024 degrees.The beam divergence angle and diameter at exit is 0.3 mrad and 7 mm, respectively with a laser beam focal spot at 50 m distance of 18 mm.The maximum measurement rate is 300 kHz.The scanner accuracy (the difference between the true value and measured value) was ± 5 mm, and the precision, or repeatability, was ±3 mm [46].Each scan collected approximately 30-33 million points of the hillslope and channel.Multiple scans collected expanded the range and density of the data points and removed or minimized shadowing effects.The additional scan locations provided overlapping data points, which enabled registration of the point clouds for developing a topographic surface of the study reach with an average point density of 20 points/m 2 .

TLS Data Acquisition
UNAVCO collected TLS data from Williams Downstream reach using a Riegl VZ-400 scanner with a Nikon D700 20 mm camera mounted on a tripod.The scanner system operated in a near-infrared wavelength (800-1000 nm) and the maximum scanner range of data collection varied from 60 to 600 m, depending on object brightness.This setup allowed a 360-degree panorama scan in the horizontal direction and 100 degree in the vertical direction, which measured at a rate of up to 122,000 measurements/sec.The data coverage depended on the incident angle of the laser, as determined by the height of the scanner above the ground and the clarity of the site with respect to the scanner (obliqueness) [46].The minimum horizontal and vertical step size is 0.0024 degrees.The beam divergence angle and diameter at exit is 0.3 mrad and 7 mm, respectively with a laser beam focal spot at 50 m distance of 18 mm.The maximum measurement rate is 300 kHz.The scanner accuracy (the difference between the true value and measured value) was ± 5 mm, and the precision, or repeatability, was ±3 mm [46].Each scan collected approximately 30-33 million points of the hillslope and channel.Multiple scans collected expanded the range and density of the data points and removed or minimized shadowing effects.The additional scan locations provided overlapping data points, which enabled registration of the point clouds for developing a topographic surface of the study reach with an average point density of 20 points/m 2 .
The first set of TLS datasets collected in April 2013 comprised six separate scans of the Williams Downstream reach.The second and third TLS data sets consisted of five separate scans of the Williams Downstream reach.These scans were collected in September 2013 and September 2014, respectively.Scanner instrument locations are illustrated on Figure 3.

Data Preparation
All the scans from Williams Downstream were registered and georeferenced.Five to six scans were registered based on tie points as reference points, which were fixed control points and managed by a tie point list (TPL) or log of the data collection to align the single scans of the stream reach.The scans were located in the same coordinate system (project coordinate system), one point cloud was fixed and the other point clouds were transformed to fit the target points of the fixed scan [31].For georeferencing the data, the scan locations were determined by registering the GPS position of the targets to place the point clouds into a global coordinate system in the WGS84 system, UTM Zone 13N.
The octree structure filtered the data points by sub-dividing the points into eight equally sized cubes.The size of the base cube was −10,000 for min x, y, h, and +10,000 for max x, y, h, and the increment value was 0.1 m.A terrain filter tool filtered and removed vegetation, from the data based on the distances of the points from an estimated ground surface.Based on these distances, the points were classified either as "terrain" or as "off-terrain."The process started from a coarse grid with 0.25 m size.For each cell, the representative point was selected based on the percentile parameter of 1%, and the tolerance factor parameter was ±0.7 m.For fine filtering, the tolerance value was 0.1 m, and the maximum slope angle was 60 degrees.
Triangulating the point clouds utilized a plane triangulation tool to create DEMs.Maximum triangle edge length was 10 m, maximum triangle tilt angle was 90 degrees and minimum triangle angle was 0 degrees.The triangulation was performed for the active channel excluding hillslopes, and with a top view for our case of study.Also, this study avoided smoothing, which would underestimate the ground surface by ~3.5 m.

Data Preparation
All the scans from Williams Downstream were registered and georeferenced.Five to six scans were registered based on tie points as reference points, which were fixed control points and managed by a tie point list (TPL) or log of the data collection to align the single scans of the stream reach.The scans were located in the same coordinate system (project coordinate system), one point cloud was fixed and the other point clouds were transformed to fit the target points of the fixed scan [31].For georeferencing the data, the scan locations were determined by registering the GPS position of the targets to place the point clouds into a global coordinate system in the WGS84 system, UTM Zone 13N.
The octree structure filtered the data points by sub-dividing the points into eight equally sized cubes.The size of the base cube was −10,000 for min x, y, h, and +10,000 for max x, y, h, and the increment value was 0.1 m.A terrain filter tool filtered and removed vegetation, from the data based on the distances of the points from an estimated ground surface.Based on these distances, the points were classified either as "terrain" or as "off-terrain."The process started from a coarse grid with 0.25 m size.For each cell, the representative point was selected based on the percentile parameter of 1%, and the tolerance factor parameter was ±0.7 m.For fine filtering, the tolerance value was 0.1 m, and the maximum slope angle was 60 degrees.
Triangulating the point clouds utilized a plane triangulation tool to create DEMs.Maximum triangle edge length was 10 m, maximum triangle tilt angle was 90 degrees and minimum triangle angle was 0 degrees.The triangulation was performed for the active channel excluding hillslopes, and with a top view for our case of study.Also, this study avoided smoothing, which would underestimate the ground surface by ~3.5 m.

Change Estimation and Analysis
This comparative analysis calculated the volume between two surfaces (base mesh, reference mesh) based on the DOD method.The two meshes were converted to raster structures with 0.1 m cell size.Three DEMs were derived from the TLS and point cloud preparation for Williams Downstream reach (April 2013, September 2013, and September 2014; Figure 3), which showed the topographic variations within the stream channels.Volumetric changes in channel sediment between scans were calculated by subtracting the DEM of the post-fire, pre-storm from the post-storm DEM [2,11].We then compared three methods, C2C, C2M, and M3C2, to estimate topographic changes between two successive scans of the same reaches.

Error Analysis
To quantify the error for two TLS scans between April 2013 and September 2013, two point clouds from the regular TLS survey and fixed targets with high resolution [31] were compared.Surveys of fixed targets were not available to evaluate the error for September 2014.Therefore, empirical error was calculated, which was the uncertainty of two point clouds from 12 immovable objects (large boulders) within the overlapping areas represented by adjacent scans (i.e., point offset in any x, y, or h direction) [41,42].This represented the uniform mean error or average overall error in point clouds including registration error for the reach.The errors in x, y, and h directions were propagated in the temporally separate scans and calculated according to Equation (2).
The DEM error was estimated for each scan (April 2013, September 2013 and September 2014) separately by comparing ten random cross sections of the channels before and after triangulation.The RMSE (Equation ( 1)) was calculated for each scan, which was used to estimate the total propagated error in a DEM of Difference (δ h(DoD) (Equation ( 3)) between two successive scans.Finally, Equation ( 4) provided an estimation of the volumetric error.

Comparison of Topographic and Volumetric Change Methods
Using the Riscan Pro DOD method to compare the DEMs of successive scans, the volume of erosion and deposition were estimated as 190 m 3 and 460 m 3 , respectively for April, 2013-September, 2013 (over an area of 4970 m 2 ), and 250 m 3 and 180 m 3 for September 2013-September 2014 (over an area of 5740 m 2 ). Figure 4 provides a schematic of where erosion and deposition occurred based on the DOD method of comparison within the stream channel for two years following the fire.
This comparative analysis calculated the volume between two surfaces (base mesh, reference mesh) based on the DOD method.The two meshes were converted to raster structures with 0.1 m cell size.Three DEMs were derived from the TLS and point cloud preparation for Williams Downstream reach (April 2013, September 2013, and September 2014; Figure 3), which showed the topographic variations within the stream channels.Volumetric changes in channel sediment between scans were calculated by subtracting the DEM of the post-fire, pre-storm from the post-storm DEM [2,11].We then compared three methods, C2C, C2M, and M3C2, to estimate topographic changes between two successive scans of the same reaches.

Error Analysis
To quantify the error for two TLS scans between April 2013 and September 2013, two point clouds from the regular TLS survey and fixed targets with high resolution [31] were compared.Surveys of fixed targets were not available to evaluate the error for September 2014.Therefore, empirical error was calculated, which was the uncertainty of two point clouds from 12 immovable objects (large boulders) within the overlapping areas represented by adjacent scans (i.e., point offset in any x, y, or h direction) [41,42].This represented the uniform mean error or average overall error in point clouds including registration error for the reach.The errors in x, y, and h directions were propagated in the temporally separate scans and calculated according to Equation (2).
The DEM error was estimated for each scan (April 2013, September 2013 and September 2014) separately by comparing ten random cross sections of the channels before and after triangulation.The RMSE (Equation ( 1)) was calculated for each scan, which was used to estimate the total propagated error in a DEM of Difference ( ( ) (Equation ( 3)) between two successive scans.Finally, Equation ( 4) provided an estimation of the volumetric error.

Comparison of Topographic and Volumetric Change Methods
Using the Riscan Pro DOD method to compare the DEMs of successive scans, the volume of erosion and deposition were estimated as 190 m 3 and 460 m 3 , respectively for April, 2013-September, 2013 (over an area of 4970 m 2 ), and 250 m 3 and 180 m 3 for September 2013-September 2014 (over an area of 5740 m 2 ). Figure 4 provides a schematic of where erosion and deposition occurred based on the DOD method of comparison within the stream channel for two years following the fire.performed with the CloudCompare C2C, C2M, and M3C2 methods (Figure 5). Figure 5A compares channel morphological changes between April 2013 and September 2013 for each method.The C2C method only represents the absolute distance difference and does not show the erosion and deposition separately.All red areas are channel variations including both gain and loss (Figure 5A).Also, since C2C compares point clouds, areas with no points are visible along the channel.The M3C2 shows the distance difference of the points and is similar to C2C method.However, the erosion (blue) and deposition (red) are separated.The C2M model is a triangulated distance difference surface that shows erosion (blue) and deposition (red) along the channel.Figure 5B shows the comparison of channel morphology changes between September 2013 and September 2014 for each method.The C2C model provides point cloud comparison results in red as an absolute value (change is not separated by erosion or deposition), indicating only a total change [38].The M3C2 shows the distance difference result of point clouds comparison and is similar to C2C (based on points comparison).However, M3C2 documents erosion and deposition along the channel separately.The C2M model compares triangulated surfaces and does not have any gaps.The result of distance difference is primarily erosion (values less than 0).Estimated net changes (deposition minus erosion) and accumulated changes (deposition plus erosion) allowed comparison across the methods of distance differencing.C2C represents accumulated changes, whereas e C2M and M3C2 represent net change.The DOD method provides both net and accumulated changes.Table 2 shows the net changes in bed elevation and accumulated changes for Williams Downstream reach.For the DOD method, volume is divided by area.The C2C, C2M, and M3C2 methods estimate a mean distance difference based on a Gauss distribution, a continuous distribution of data as a symmetrical bell-shaped graph.These results provide a basic comparison of different methods for estimating surface variations.The distance differences in All three models in Figure 5A,B show the same erosion and deposition patterns observed in Figure 4. Noticeable "hot spots" or regions with relatively large deposition or erosion (0.9-1.5 m) indicated are evident (Figure 5A,B), which are artifacts of areas with spares points due to shadows from natural features such as boulders and trees, resulting in gaps in the point cloud.Collecting several scans with different perspectives partially removes shadowing effects.However, areas with low point density may still exist.These areas contribute to the error in DEMs and estimation of channel volumetric changes by introducing errors in data point interpolation.
Estimated net changes (deposition minus erosion) and accumulated changes (deposition plus erosion) allowed comparison across the methods of distance differencing.C2C represents accumulated changes, whereas e C2M and M3C2 represent net change.The DOD method provides both net and accumulated changes.Table 2 shows the net changes in bed elevation and accumulated changes for Williams Downstream reach.For the DOD method, volume is divided by area.The C2C, C2M, and M3C2 methods estimate a mean distance difference based on a Gauss distribution, a continuous distribution of data as a symmetrical bell-shaped graph.These results provide a basic comparison of different methods for estimating surface variations.The distance differences in elevation from each method has a relatively small absolute value range between 0.009 and 0.150 (Table 2).The net changes during September 2013-September 2014 for DOD and M3C2 methods are the same (0.054 m), and the maximum difference is 0.020 m.The only significant difference exists in accumulated changes for September 2013-September 2014, where the value for the C2C method (0.140 m) is almost twice the DOD method (0.075 m).Since the C2C technique is based on the measurement of the closest point distance, the distance estimation between two point clouds is sensitive to the point density.It is only appropriate for very dense point clouds and not an accurate distance measurement method [34,38].The mean uncertainty of two point clouds (April 2013 and September 2013) was estimated based on a comparison of two sets of points from the TLS survey and fine scanning (larger number of points collected) of fixed targets (Figure 3).The resolution of the fine scans from the fixed targets in TLS scanning using a Riegl Scanner was 9 to 555 points, whereas the resolution of the fixed targets for the same area with fine scanning mode was 2106 to 319,394 points.Thus, the fine scans from the fixed targets were considered as the reference data and presumably error free [35].However, it should be noted that the resolution of scans can be different based on the TLS instrument and analysis software system used.
The error analysis was performed to evaluate uniform error for the September 2014 scan data with no fixed target by comparing points from stationary rocks in overlapping area (Figure 3).RMSE was estimated for each scan in x, y, and h directions (Table 3).The values range from 0.003 m to 0.008 m for all three scans and demonstrate that the TLS survey methodology has a low overall uniform error regardless of the interpolation technique.The errors in x, y, and h directions were propagated and estimated as 0.015 m for April 2013-September 2013, and 0.016 m for September 2013-September 2014 (from Equation ( 2)).This procedure yields an estimate of minimum LoD.The LoD min of points is important in the C2C and M3C2 methods, which are based on comparing point clouds to estimate the distance differences.Consequently, in the C2C and M3C2 methods, any change estimation should be subtracted from 0.015 m (for April 2013-September 2013) or 0.016 m (for September 2013-September 2014) to incorporate the uncertainty.Also, change estimated below 0.015 m or 0.016 m should be excluded from analysis and results.

Error Analysis after Triangulation (DEM Error)
The interpolation error was estimated for each scan (April 2013, September 2013, and September 2014) by comparing the actual points with its triangulated surface (Table 4).The estimated RMSE of the DEMs is 0.009-0.014m, which is relatively larger compared to the RMSE of the point clouds (0.003-0.008 m; Table 3).This result indicates that triangulation of the point clouds increases the error, specifically in areas with large gaps due to the interpolation of points.The triangulation error can depend on the algorithm used for converting point cloud to DEM.However, the same workflow to estimate errors is encouraged for all processing algorithms.The propagation of DEM error for each scan utilized Equation (3) to calculate the potential error in distance difference of two successive scans, which is the LoD min for DOD and C2M methods.The LoD min for April 2013-September 2013 is 0.018 m, and for September 2013-September 2014 is 0.017 m.Similar to LoD min calculated for points, change estimates for DOD and C2M should be subtracted from 0.018 m (for April 2013-September 2013) or 0.017 m (for September 2013-September 2014) to incorporate uncertainty.Additionally, the DOD and C2M analysis results should exclude any changes estimated below 0.018 m or 0.017 m.  4) to determine the uncertainties found in the volumetric calculations.Considering the entire volume of the channel, the volumetric error (δ υ ) is 89.5 m 3 (0.03%) for April 2013-September 2013 for a total volume of 2.92 × 10 5 m 3 , and 97.6 m 3 (0.032%) for September 2013-September 2014 for a total volume of 3.02 × 10 5 m 3 .The volumetric error of the erosion and deposition was estimated as a proportion of the channel volume and is 190 ± 0.06 m 3 and 460 ± 0.14 m 3 , respectively, for April 2013-September 2013, and 250 ± 0.08 m 3 (erosion) and 180 ± 0.06 m 3 (deposition) for September 2013-September 2014.Since the volumetric error is a proportion of the volume, the error increases with increasing erosion or deposition volume.
The volumetric errors are 0.03%-0.032% of the volumetric changes and attributed to existing errors in point clouds and DEMs, which result from triangulation and potential error from data preparation and processing.In comparison to the erosion and deposition volume, these volumetric errors are negligible.Despite the errors in the volumetric change evaluations, TLS provides volume changes within terrain channels with relatively low uncertainty.

Discussion of Errors
In prior works, researchers have applied Terrestrial Laser Scanning to watersheds after disturbances; however, minimal guidance is available to report errors and uncertainty associated with estimates.Therefore, this step is often left out [3].This section compares error analysis methods developed in this study with methodologies and results from prior works.
Staley et al. [2] used TLS before and after storms to estimate volume of erosion and deposition in a small (0.01 km 2 ) headwater sub-basin of the Arroyo Seco after the 2009 Station Fire in Los Angeles, California (USA).This investigation incorporated 11 control points (locations scanned in situ and assumed to have minimal errors) within the basin to estimate errors associated with TLS data.Topographic changes ranged from −1.34 ± 0.040 m (erosion) to 1.77 ± 0.040 m (deposition).The minimum level of detectable topographic change was 0.011 m and the authors excluded areas with changes below the 0.011 m threshold from the volumetric analysis [2].Rengers et al. [13] collected five successive TLS datasets in the Colorado Front Range (USA) and estimated an average erosion along the main channel of 0.073 m with minimal deposition (less than 1 cm).Following Klapinski et al. [42], Rengers et al. [13] estimated the uncertainty of the maximum offset of each TLS point cloud by comparing stationary objects (i.e., large rocks and boulders) with overlapping areas through a fuzzy inference system (FIS) based on Wheaton [47].An uncertainty between 0.005-0.02m was estimated as the minimum level of detection, which was assigned to each pixel in the DEM and used as a threshold to remove data from the analysis.Another study by Florsheim et al. [1] incorporated multiple scan positions to construct changes in topography for 176 m of stream channel before and after a storm season, which was used to calculated volumetric change following the 2013 Springs Fire in southern California (USA).The authors utilized error analysis methods similar to the current study to quantify uncertainty in the triangulated surfaces compared to the point cloud data.Florsheim et al. [1] estimated a propagation error of 0.011 m, or ~16% of the total volumetric change resulting from deposition, which is larger than the errors estimated in this study (0.030%-0.032%).This may be due to the number of scan positions used and the amount of volumetric change observed within the reach.
In general, studies [1,2,13] that quantified errors in TLS datasets often utilized a variety of methods.However, the estimated uncertainty in elevation change ranged from 0.005-0.02m, which is comparable to the errors quantified in this paper (0.017-0.018 m).Thus, the current study reviewed and discussed the C2C, C2M, M3C2, and DOD methods to provide a workflow for quantifying topographic changes and errors for analyses with multiple datasets.

Conclusions
This paper provides a workflow to guide effective and efficient methods to process and interpret LiDAR information and to quantify associated uncertainty.The application of the workflow to TLS data collected after the 2012 Waldo Canyon Fire provided a demonstration of quantifying error in four methods for estimating topographic and volumetric changes in channel sedimentation.This comparison advances our understanding of the extent of error propagation during data preparation and processing.
To summarize the lessons learned from this analysis: 1.
The error analysis quantifies a minimum level of detection, which is important for considering the uncertainty in both point clouds and DEMs in TLS processing, analysis results, and conclusions about channel alteration.

2.
The error analysis demonstrates rapid propagation of error from individual DEMs into DEMs of difference, highlighting the need for consideration.Estimating errors not only for elevation (h), but also in the horizontal directions (x and y) is necessary in point clouds analysis, especially in complex and rugged terrain.However, errors for each single dataset are significantly smaller in point clouds before triangulation (0.003-0.008 m) than in DEMs (0.009-0.014 m) as converting data points within a LiDAR point cloud to DEMs increases the error by interpolation.The propagated error for comparing multiple datasets recording surfaces over time is similar for point clouds (0.015-0.016 m) and DEMs (0.017-0.018 m).

3.
Approximately 0.030%-0.032% of the estimated erosion and deposition results from errors that exist in the TLS data set, indicating a lower potential for uncertainty for estimating volumetric changes.4.
Triangulation of the bare Earth points in the point cloud increases error, especially in complex topography, where occlusions are typically larger than flat surfaces.Therefore, techniques that estimate topographic changes that are solely based on point comparison (i.e., C2C and M3C2) are more reliable than those that interpolate surfaces (i.e., C2M and DOD).

5.
The C2C method, the absolute distance difference, does not negative and positive values, making data interpretation challenging for studies that are interested in erosion or deposition.In general, C2C compares the closest point to the reference point and provides a simple method for rapid change estimation.6.
The M3C2 method extracts erosion and deposition rates from point cloud comparisons and calculates the most accurate and reliable results for the comparison of two surfaces (of the methods reviewed in this work).The M3C2 method does not use interpolation, which minimizes the errors present, especially in complex terrain.7.
The C2M and DOD methods are based on triangulated surfaces.The C2M method is recommended to estimate distance difference and provide elevation, whereas the DOD method is recommended for volumetric changes.The DOD estimates are similar to the M3C2 method.The findings of this review and analysis suggest DOD for estimating volumetric changes, and M3C2 for estimating distance differences in elevation in point clouds.
TLS applications have advantages for assessing topographic changes, and for providing three-dimensional surfaces to detect changes with relatively high accuracy.Although this paper compared several methods for topographic and volumetric change and error estimation based on TLS, these methods are also applicable to other high-resolution methods such as aerial-based LiDAR or SfM photogrammetry data analyses.The suggested workflow is especially necessary as point cloud data are increasingly incorporated into research that seeks to improve understanding of short-and long-term changes after disturbance.

Figure 1 .
Figure 1.Four-part workflow for (I) LiDAR data acquisition, (II) preparation and processing, (III) change estimation and analysis, and (IV) error analysis.

Figure 1 .
Figure 1.Four-part workflow for (I) LiDAR data acquisition, (II) preparation and processing, (III) change estimation and analysis, and (IV) error analysis.

Figure 2 .
Figure 2. Williams watershed is located northwest of Colorado Springs (located at the diamond in the map inset).Elevation is shown and tributaries are highlighted in black (source: USGS National Map 30 m digital elevation map).The rectangle near the outlet of the watershed shows the approximate location of Williams Downstream reach (see Figure 3 for an enlargement of the study reach).

Figure 2 .
Figure 2. Williams watershed is located northwest of Colorado Springs (located at the diamond in the map inset).Elevation is shown and tributaries are highlighted in black (source: USGS National Map 30 m digital elevation map).The rectangle near the outlet of the watershed shows the approximate location of Williams Downstream reach (see Figure 3 for an enlargement of the study reach).
The first set of TLS datasets collected in April 2013 comprised six separate scans of the Williams Downstream reach.The second and third TLS data sets consisted of five separate scans of the Williams Downstream reach.These scans were collected in September 2013 and September 2014, respectively.Scanner instrument locations are illustrated on Figure3.

Figure 3 .
Figure 3. Scans of Williams Downstream reach for April 2013, September 2013, and September 2014.The locations of the fixed targets and boulders used in the error analysis are represented by solid black circles.Scanner instrument locations are represented by triangles.

Figure 3 .
Figure 3. Scans of Williams Downstream reach for April 2013, September 2013, and September 2014.The locations of the targets and boulders used in the error analysis are represented by solid black circles.Scanner instrument locations are represented by triangles.

Figure 4 .
Figure 4. DEM of Difference (DOD) method used to compare successive scans of Williams Downstream reach.Blue and red denote locations of erosion and deposition, respectively.Distance difference estimation, which is a quantification of the absolute value of elevation changes within the channel caused by both erosion and deposition, for Williams Downstream reach was Remote Sens. 2018, 9, x FOR PEER REVIEW 12 of 19

Figure 5 .
Figure 5.Comparison of successive scans of Williams Downstream reach, showing the distance difference in elevation using the C2C, C2M, and M3C2 methods for channel changes between (A) April 2013-September 2013 and (B) September 2013-September 2014.

Figure 5 .
Figure 5.Comparison of successive scans of Williams Downstream reach, showing the distance difference in elevation using the C2C, C2M, and M3C2 methods for channel changes between (A) April 2013-September 2013 and (B) September 2013-September 2014.

Table 1 .
LiDAR point cloud analysis methods and products used in the current study.

Table 2 .
Net changes (deposition minus erosion) and accumulated changes (deposition plus erosion) for Williams Downstream reach April 2013-September 2013 and September 2013-September 2014.

Table 3 .
TLS point cloud uniform uncertainty for the lateral (x, y) and elevation (h) directions for the TLS scans within Williams Downstream reach.

Table 4 .
Root mean square errors for Williams Downstream DEMs after triangulation.