Comparing Two Methods of Surface Change Detection on an Evolving Thermokarst Using High-Temporal-Frequency Terrestrial Laser Scanning, Selawik River, Alaska

Abstract

Terrestrial laser scanners (TLS) allow large and complex landforms to be rapidly surveyed at previously unattainable point densities. Many change detection methods have been employed to make use of these rich data sets, including cloud to mesh (C2M) comparisons and Multiscale Model to Model Cloud Comparison (M3C2). Rather than use simulated point cloud data, we utilized a 58 scan TLS survey data set of the Selawik retrogressive thaw slump (RTS) to compare C2M and M3C2. The Selawik RTS is a rapidly evolving permafrost degradation feature in northwest Alaska that presents challenging survey conditions and a unique opportunity to compare change detection methods in a difficult surveying environment. Additionally, this study considers several error analysis techniques, investigates the spatial variability of topographic change across the feature and explores visualization techniques that enable the analysis of this spatiotemporal data set. C2M reports a higher magnitude of topographic change over short periods of time (∼12 h) and reports a lower magnitude of topographic change over long periods of time (∼four weeks) when compared to M3C2. We found that M3C2 provides a better accounting of the sources of uncertainty in TLS change detection than C2M, because it considers the uncertainty due to surface roughness and scan registration. We also found that localized areas of the RTS do not always approximate the overall retreat of the feature and show considerable spatial variability during inclement weather; however, when averaged together, the spatial subsets approximate the retreat of the entire feature. New data visualization techniques are explored to leverage temporally and spatially continuous data sets. Spatially binning the data into vertical strips along the headwall reduced the spatial complexity of the data and revealed spatiotemporal patterns of change.

1. Introduction

In the last decade, techniques to rapidly acquire high density topographic data have proliferated; however, the tools to adequately analyze change within these complex data sets have only begun to emerge. Topographic change detection has evolved from techniques that relied on visual comparisons (e.g., Webb et al.[1], Williams [2], Tape et al.[3], Bierman et al.[4], Kaab et al.[5], Nichol and Wong [6] and Frankl et al.[7]) to survey-based techniques (e.g., Mackay [8], Ruhlman and Nutter [9] and Wheaton et al.[10]). Currently, airborne laser scanners (ALS) and terrestrial laser scanners (TLS) enable the acquisition of topographic data sets with postings between 1 m and <5 cm. ALS systems allow for the rapid collection of regional topographic data, while TLS systems allow for high density, local topographic surveys of complex terrain. Additionally, the advent of GIS has allowed the analysis of spatially extensive and dense topographic data sets (e.g., Wheaton et al.[10]). The increase in topographic complexity captured at fine spatial scales with ALS and TLS surveys underscores the need for more sophisticated point cloud processing and change detection methods.

Both ALS and TLS instrumentation support the filtering of laser returns to remove vegetation via waveform or multiple discrete return processing. All returns can be classified by height or physical geometric parameters, which is useful both for isolating points and identifying topographic form and change (e.g., Lane et al.[11] and Charlton et al.[12] for ALS examples and Gordon et al.[13], Rosser et al.[14], Girardeau-Montaut et al.[15], Heritage and Hetherington [16], Abellán et al.[17,18] Brodu and Lague [19] and Kociuba et al.[20] for TLS examples). Many early ALS and TLS topographic change detection techniques relied on the comparison of digital elevation models (DEMs) for use in change detection analyses using a GIS (e.g., Charlton et al.[12] and Krieger [21]). Subsequent change detection methods, such as direct cloud to cloud comparisons (C2C) and cloud to mesh (C2M) techniques, have been developed to compare topographic data sets. C2C directly compares two point clouds, while C2M compares a point cloud to a triangulated model of a feature. Recently, Lague et al.[22] has proposed the Multiscale Model to Model Cloud Comparison (M3C2) algorithm, which can carry out topographic change detection using repeat TLS point clouds with minimal cleaning and processing of the scan data. This study seeks to compare and contrast two TLS topographic change detection methods, C2M and M3C2, in a complex topographic environment.

1.1. Change Detection via TLS

Change detection analyses using TLS data are difficult because of (1) the complexity and richness of the point clouds representing the natural environment and (2) because TLS point clouds from different temporal epochs or physical scan positions may not overlap closely enough for accurate point to point comparison. Additionally, the software packages used to acquire, process and analyze TLS data are often optimized for the built environment and lack analytical tools relevant to natural landscapes. Many solutions have been proposed to account for these shortcomings, such as gridding point clouds into DEMs. Gridded data sets can be compared to produce a DEM of difference (DoD), which highlights areas of loss or accumulation [10,15,23–26]. Other methods include C2M, where surface models are created from point cloud data sets via meshing or triangulation and compared with subsequently gathered point data sets (e.g., Abellán et al.[17]; Abellán et al.[18]; and Olsen et al.[27]). However, isolating regions of interest and creating DEMs or surface models is time-consuming and requires interpolation [22]. Furthermore, a topographic data reduction occurs when point cloud data is reduced from 3D to 2.5D. C2M comparisons allow non-horizontal surfaces to be analyzed without the rotation required in DoD analyses; however, surfaces with significantly different orientations should be isolated before analysis, so that they may be treated separately.

Cloud to cloud methods have also been used to directly compare point clouds (e.g., Girardeau-Montaut et al.[15]), but they often do not return signed displacements, which are required for many applications. Although both C2M and rotated DoD techniques enable comparisons of complex point cloud scenes, they require a high degree of data processing and isolation to work effectively. M3C2 directly compares two point clouds and conducts change detection on entire scenes of points with minimal manual processing [22]. The way each change detection method incorporates uncertainty to differentiate true from perceived topographic change is also an important factor when comparing these methods. In order to assess the relative merits of C2M and M3C2, we apply both to a TLS data set collected at a challenging and complex field site in northwest Alaska.

To clearly explain both techniques, the following nomenclature will be used to reference TLS products. TLS measurement epochs refer to individual point clouds collected sequentially through time and are denoted C₁, C₂,...,C_n. Surface meshes/models derived from point clouds are denoted similarly as M₁, M₂,...,M_n. When the C2M and M3C2 methods are described, data sets consisting of sequentially gathered point cloud pairs are described as before (C_b) and after (C_a), and before and after surface meshes/models are referred to as M_b and M_a, respectively. M3C2 makes use of core points [19], which are a user-defined subset of a full point cloud used to reduce computational intensity. Core points are denoted with the addition of c to an epoch’s subscript. C_1c is then a lower density subset of C₁. When computations are performed for C_1c, data are returned for every point in C_1c; however, all data in C₁ in the vicinity of the points contained within C_1c are used to inform the calculations. This allows the number of points requiring analysis to decrease, while still leveraging the full density of C₁. For example, a core point file of one million points could be generated from a scan with six million points and could be compared to a similarly dense data set from a second scan epoch using M3C2. The analysis run using the core point file calculates displacements for each point in the core point file and takes the full density data sets of the before and after scans into account. This results in the calculation of one million displacements rather than the six million, had every point in the full data set been fully analyzed. Lague et al.[22] provide a detailed account of how core points are used in M3C2.

1.1.1. C2M Analyses

C2M comparisons have been used by many workers to detect a change between sequentially gathered point clouds [17,18,27]. Consider two sequential point clouds of the same feature, C_b and C_a (Figure 1(A)). C_b is triangulated to create M_b, a reference mesh that models C_b (Figure 2(B)). Triangulation methods range from Delaunay triangulation, where long skinny triangles are minimized, to more customized methods, where a vantage location in the point cloud scene and a look point in C_b are selected to specify how C_b is viewed by the meshing algorithm during triangulation. Additionally, a maximum allowable triangle edge length (MATEL) may be specified to generate a mesh of complex vertical or overhanging surfaces. Setting a MATEL also allows C2M-type comparisons to ignore areas in C_b that are topographically shadowed or have sparse point density during the creation of M_b; however, if the width of topographic shadows or gaps in C_b are below the MATEL, these areas will be interpolated as a continuous surface by the facets of the resulting mesh.

For every point, i, in C_a, the distance, L_C2M, to the nearest vertex or facet in M_b is computed and stored as an attribute of i (Figures 1(A) and 2(B)). A maximum allowable distance can be specified to eliminate unrealistic results when the point clouds being compared have different extents or when data sets are spatially discontinuous. This method has been used with success to investigate structural and surface deformation [27,28]. Collins et al.[29] have used similar techniques to investigate land surface change in the Grand Canyon, AZ. It should also be noted that Wheaton et al.[10] and Glennie et al.[30] used triangulated irregular networks, a meshing method, as an intermediate step when creating DEMs of their study area.

Lague et al.[22] note that C2M-type comparisons provide good approximations of surface change, yet they do not take into account the local orientation of the surfaces represented by the point clouds. In this study, surface orientation was forced towards the center of the RTS feature to remove ambiguity. Not taking point cloud surface orientation into account causes the displacement calculated between point clouds to vary depending on point density, topographic shading and changes in surface roughness [22,31]. If the change between C_b and C_a is spatially uniform and the exact same scan position is used for both epochs, these shortcomings are minimized; however, this is rarely the case. Furthermore, Lague et al.[22] note that C2M and DoD change detection methods are labor- and time-intensive, because areas of interest must be carefully isolated and treated in the TLS data processing workflow to create reference meshes and DEMs (Figure 2(A,B)).

1.1.2. M3C2 Analyses

The M3C2 technique allows rapid analysis of large point clouds with complex surfaces that span a range of surface orientations. Lague et al.[22] combine the M3C2 technique with the point cloud classification and vegetation removal technique presented in Brodu and Lague [19] to assess change on a variety of surfaces. A short synopsis of the technique is provided here.

The M3C2 algorithm can be broken into two steps: point normal estimation and difference computation (Figure 2(C)). Users may specify if local point normals are calculated or if normals are fixed in either the horizontal or vertical direction. Horizontal point normals allow true horizontal erosion rates to be calculated from M3C2 analyses, while vertical normals allow M3C2 data to be used for strictly vertical erosion and aggradation measurements. Point normals are computed for C_bc; however, C_b in its entirety can be used as C_bc if a greater spatial resolution of normal calculation is desired. For every point, i, in C_bc, the neighboring points within a radius of D/2 are used to calculate a best fit plane (Figure 1(B)). D is a user-specified scale appropriate for the scene (see

Table 1. Data used with Equation (1) and the resulting RMS error estimate. Registration errors in 2012 necessitated the use of a surface matching algorithm to re-georeference a portion of the point cloud data, which resulted in a larger georeferencing error. The 2012 GPS error is also larger than the 2011 GPS error. M3C2 spatially variable confidence interval (SVCI) values are the averages of the confidence intervals computed for all displacements by the algorithm.

Year	GPS (2σ) (m)	Instrumental * (1σ) (m)	Georeferencing (2σ) (m)	RMS Error (2σ) (m)	Empirical Error (m)	M3C2 SVCI (2σ) (m)
2011	0.008	0.01	0.015	0.02	0.028	0.052
2012	0.014	0.01	0.022	0.03	0.019	0.069

^*Instrumental error provided by the manufacturer lists the error as 0.005 m at 100 m [55]; however, the slump headwall is, on average, 200 m from the instrument, so the published value was doubled for this study.

in Lague et al.[22]). The pole of the best fit plane is the normal for i. Lague et al.[22] note that the appropriate normal estimation scale, D, is based on point cloud density and local roughness.

Following point normal calculation, a cylinder with its axis aligned with the normal for i is projected toward C_a (Figure 1(B)). The cylinder has a radius of d/2, where d is the scale specified by the user. The maximum length of the cylinder is specified to reduce computational intensity. All of the points in C_b and C_a that reside in the cylinder are spatially averaged to calculate mean surface positions, i₁ and i₂, respectively. The distance, L_M3C2, between i₁ and i₂ is recorded as an attribute of i. Additionally, the local and apparent roughness, σ₁ and σ₂, of C_b in the vicinity of i₁ and C_a in the vicinity of i₂ are calculated along the normal vector. These values are used in conjunction with the error associated with registering C_b and C_a to a common coordinate system to determine the change detection confidence interval for i. The calculation and treatment of the M3C2 change detection confidence interval is discussed in greater detail in Section 2.3.1.

1.2. Objective

Detecting topographic change in diverse natural environments is a persistent problem for earth scientists that has increased in complexity as new survey methods enable the rapid acquisition of complex topographic data. The magnitude of this problem increases further at sites where it is difficult to maintain georeferencing control. Previous work to compare the C2M and M3C2 techniques have relied on simulated point clouds to test differences in the methods [22]. This study builds on the foundation emplaced by Lague et al.[22] and seeks to (1) explore the limits of topographic change detection by comparing C2M and M3C2 using TLS scans gathered in an environment problematic for high-resolution topographic change detection; (2) compare TLS error analysis methods, which assist in differentiating real topographic change from noise in survey data sets; and (3) leverage an unusually long repeat TLS data set to investigate the spatial and temporal variability of topographic change.

1.3. Site Description and Significance

In the Arctic, topographic change detection has long been employed to assess landscape change and landform evolution. Mackay [8] resurveyed a retrogressive thaw slump (RTS) to investigate its growth; Lewkowicz [32] employed specialized ablatometers to measure the retreat rate of a similar RTS feature; and, more recently, Matsuoka [33] delineated a variety of detection methods to monitor periglacial processes. Krieger [21] used both ALS and TLS methods to assess thermokarst form and evolution on Alaska’s north slope.

RTS are found throughout the Arctic in a variety of environments, such as coastlines, river banks and lake shores [34–36]. These features mobilize large quantities of sediment, soil carbon and nutrients to downstream environments [36–38]. Large permafrost disturbances, such as RTS, also alter the permafrost CO₂ flux and change the surface albedo of the tundra [38]. The frequency and magnitude of these features are anticipated to increase as the Arctic experiences enhanced climate change via a variety of positive feedback mechanisms [37,39]. Over the last century, the Arctic has warmed 0.09 °C/decade, which is 0.03 °C/decade greater than the Northern Hemisphere [40]. A better understanding of the mechanisms and environmental factors that drive RTS growth are needed, because RTS features have a high potential to impact downstream ecosystems and communities.

Although RTS features can be quite large, the individual processes that drive their growth are small. As the vertical RTS headwall thaws, material is destabilized and fails. Failures are usually water-rich from a combination of active layer surface overflow and the melt of interstitial ice within the headwall. In localized areas of the headwall, material may dry as ice melts and water leaves. The dry material then crumbles. Failures are typically <10 cm thick and can range in area from <0.1 m² to 4 m². Sapping processes are not possible on the headwall, because the permafrost is impermeable. TLS scans of an RTS headwall are able to integrate the spatial complexity of these ablation processes. Failed material accumulates at the base of the headwall and is evacuated via episodic mudflows or continuous streams that traverse the RTS floor. The transport of water and sediment away from the headwall prevents burial and stabilization. Even a small build up of material over the permafrost headwall would insulate it and halt ablation [32]. Barnhart [41] provides a more detailed discussion of the mass loss processes through which the Selawik RTS retreats and the hydrometeorological drivers that affect the retreat rate and form of the feature.

RTS are difficult to survey for a variety of reasons. (1) Tundra surface properties (dynamic frost table, extensive soft ground cover, absence of stable ground) destabilize reflective TLS target tripods used for georeferencing control; (2) shifting tripods require frequent TLS target resurveys to minimize georeferencing error; (3) remote Arctic locations require GPS baselines on the order of 225 km to process position solutions, which increases position uncertainty; and (4) inclement weather and the wet surfaces on RTS features can cause poor and/or spatially discontinuous TLS returns, which lowers confidence in the displacements calculated between point clouds. These factors, combined with a rapidly evolving feature, make the Selawik RTS an ideal natural laboratory to explore the limits of TLS change detection in an environment with a high degree of survey uncertainty.

Our study deals exclusively with topographic change detection and error analysis using a 58 scan repeat TLS data set of the Selawik RTS headwall in northwest Alaska (Figure 3). The Selawik RTS is a south-facing feature located at 66.5° N, about 100 km inland and upriver of Selawik Village, AK, at the boarder between the continuous and discontinuous permafrost zones [42]. The region is largely tussock tundra with some stands of willow and alder, as well as black spruce in the sheltered areas near the river and heath communities in the uplands [43]. The feature initiated in 2004 on a bluff on the north bank of the Selawik River [44] and has retreated headward at a rate of approximately 20 m/yr. In planform, the Selawik RTS is subcircular with an arcuate headwall approximately 250 m wide and 15 m tall (Figure 3(C)). The river bluff and headwall are composed of glacial diamict [45] with a silty upper loess layer [46]. Ice wedges, lenses and epigenetic ice are visible in the RTS headwall. The water and sediment from headwall failures not deposited on the slump floor drain through a narrow canyon near the position of the original river bluff. The efflux from the RTS has built a large fan into the Selawik River and, at times, has even dammed the river (Figure 3(B)). The Selawik RTS is active from late spring through early fall, when the maximum daily air temperature is above 0 °C. During this period, even at midsummer, the site experiences significant day-to-night variability in air temperature and insolation.

2. Methods

2.1. Data Acquisition and Processing

We used a Riegl VZ-400 TLS to collect point clouds of the south-facing Selawik RTS headwall at a range of 200–250 m over two weeks during the summer of 2011 and four weeks during the summer of 2012. Scan point spacing in 2011 and 2012 was 0.05 m and 0.045 m at a range of 200 m, respectively. Scans of the RTS headwall occurred at 9:00 and 21:00 to capture nighttime and daytime change, respectively. We used the same scan position in 2011 and 2012 for all data collection to minimize differences in topographic shading between the two data sets (triangle in Figure 3(C)). We applied an atmospheric correction to the instrument prior to data collection to account for variations in air temperature, relative humidity and air pressure. We used both the C2M and the M3C2 methods to measure spatially distributed headwall retreat rates (Figure 2).

TLS scans were georeferenced (Figure 2(A)) using a reflective target array deployed around the Selawik RTS (Figure 3(C)). We resurveyed the targets daily and frequently leveled them to minimize georeferencing errors due to shifts in the reflector array. This study used a Trimble 5700 post-process GPS receiver and antenna (PP-GPS) to measure the position of each reflector. We installed a continuously running GPS base station on a gravel bar approximately 0.75 km from the feature to provide a local static correction for the GPS data. In 2011 and 2012, we used daily occupations of 20 and 30 min at each reflector, respectively. Increasing the GPS occupation length between 2011 and 2012 resulted in a 0.003 m (2σ) decrease in the scan georeferencing error.

For each field campaign, we located the GPS base station using averaged position solutions from OPUS [47]. We used the OPUS coordinates to correct the local GPS base station observations with a position calculated using several published GPS stations in the region. We then used the corrected local base station position to post-process GPS solutions for the reflector target array using Trimble Geomatics Office and Trimble Business Center in 2011 and 2012, respectively [48,49]. It should be noted that long PP-GPS occupations in an unstable setting, such as the Selawik RTS, could reach a point of diminishing returns in precision with exceedingly long occupations. During very long occupations, sensitive target tripods may shift, while the GPS antenna is still attached and collecting data, resulting in GPS position processing errors or less precise PP-GPS position solutions. At more stable sites, it may still be possible that longer PP-GPS occupations would further increase georeferencing precision. Additionally, Oldow [50] found that constraining the geometry of the TLS target array with a total station and geolocating it with GPS increases georeferencing precision.

RiScan Pro [51] was used to georeference each scan into an Earth-centered Earth-fixed (ECEF) coordinate system (Figure 2(A)). Table 1 shows the GPS and georeferencing errors from 2011 and 2012. Georeferenced point clouds were exported from RiScan Pro at full density and converted to Universal Transverse Mercator (UTM) coordinates using a series of MATLAB functions [52] before processing and analysis using I-Site Studio [53] and M3C2. The error introduced by this conversion is not considered, because all scans received the same treatment. In 2011, scans were initially georeferenced directly into UTM coordinates; however, it was found that the ECEF georeferencing process decreased the 2011 georeferencing error by 0.028 m (2σ). In 2012, scans that showed significant georeferencing errors required re-georeferencing with the Surface Registration tool in I-Site Studio (Figure 2(A)) [54]. This tool worked best when we used only reflector faces, because they are the most stable features in each point cloud. The necessity of this type of secondary georeferencing lead to a larger total georeferencing error in 2012 than in 2011 (Table 1).

Georeferenced point clouds of the entire RTS were cropped to a long swath along the RTS headwall (white polygon in Figure 3(C)). This was done to remove areas with overhangs and the bottom of the headwall, which is variably covered with debris, depending on the level of material on the RTS floor. Cropping also ensures that the same region of the headwall is analyzed throughout each season of data. Each point cloud was then checked for overhanging vegetation and spurious points.

2.2. Change Detection Method I: C2M

We used I-Site Studio’s Spherical Triangulation tool with a MATEL of 0.25 m and fixed vantage and look points to mesh processed point clouds (Figure 2(B)). We chose these parameters to minimize interpolation in areas of low point density and to preserve holes in each data set, due to either poor returns or topographic shadowing. Specifying these stringent meshing parameters minimized the amount of topographic change masked by interpolation. Change detection analyses used I-Site Studio’s Color by Distance tool (CbD) with a maximum allowable distance of 0.3 m. The CbD tool uses the C2M technique described in Section 1.1.1. CbD accepts both mesh and point data; however, at least one mesh must be supplied for signed distance calculations. When two meshes are supplied, the vertices of M_a are treated as points resulting in a C2M-type analysis. When CbD is used, points in C_a are colored based on the distance between M_b and C_a. CbD returns a text file of the vertices in C_a with a signed distance attributed to each point in the file (Figure 4(A)). If the maximum distance is exceeded, a no data value is assigned to that point. Assigning a maximum distance for this analysis keeps CbD from calculating distances between the two clouds for regions that have holes or poor return density. Change detection analyses were implemented sequentially, such that when M₁ and M₂ were compared, M₁ became M_b and M₂ became M_a. For the subsequent comparisons, M₂ then becomes M_b and M₃ becomes M_a, and so forth. Two error analysis techniques were used to calculate a change detection threshold to be used in conjunction with CbD to differentiate noise from measurable topographic change.

2.2.1. Additive RMS Error Analysis

Collins et al.[29] proposed an additive RMS error analysis technique to estimate the accumulated errors associated with (1) the GPS solutions used for georeferencing; (2) the instrument; (3) the georeferencing itself; and (4) the registration of multiple point clouds. This error estimate can then be used to assess whether measured changes between two point clouds are significant. For this study, Collins et al.[29] equation was modified, because only one scan position was used to generate all point clouds, which allows the registration error to be dropped from the equation such that:

$$E_{\mathit{total}}=\sqrt{{E_{\mathit{gps}}}^2+{E_{\mathit{inst}}}^2+{E_{\mathit{georef}}}^2}$$

(1)

where E_gps is the error associated with the GPS solutions used to georeference the point cloud, E_inst is the error associated with the instrument used and E_georef is the error associated with georeferencing the point clouds (Table 1). The form and components of this equation are similar to the one used by Brasington et al.[23] for propagating uncertainty to change detection work with DoD analysis; however, the equation of Collins et al.[29] deals more explicitly with the sources of error involved in point cloud change detection. The uncertainty of the OPUS positions used in 2011 and 2012 are not considered in this error analysis, because inter-year comparisons of the feature are not undertaken in this study.

2.2.2. Empirical Error Analysis

The error associated with C2M was also assessed using a single vertical subset of the RTS headwall that only experiences retreat (yellow box in Figure 3(C)). The region was isolated in the point cloud data sets and analyzed using the C2M method described in Section 2.2. Any positive differences calculated between subsequent measurement epochs correspond to errors in the GPS-instrument-georeferencing system. Abellán et al.[17] used a similar technique to estimate the error in their system for monitoring rock fall events. The mean of the positive errors from the 2011 and 2012 (Table 1) agree with the results of the additive RMS error analysis, suggesting that both error estimates do a good job of aggregating the total error of the system.

2.2.3. Application of the Error Analysis Results

A MATLAB script was used to uniformly change measurements in the CbD output that were less in magnitude than the additive RMS error (Table 1) to zero, to represent no detectable change, because they fall below the calculated limit of detection (Figure 4(B)). The mean of the remaining significant differences was used in conjunction with the elapsed time between M_b and C_a to calculate the spatially averaged retreat rate of the headwall (m/d) during that period. Changing insignificant displacements to zero removes displacements due to noise from the system, which decreases the magnitude of the average change observed between scan epochs.

2.3. Change Detection Method II: M3C2

The processed point clouds (Figure 2(A)) were exported as xyz text files. Each point cloud was run sequentially using the batch file provided with the M3C2 algorithm using the following parameters: 0.3 m for D, 0.5 m for d, a maximum displacement distance of 0.3 m and total registration errors of 0.02 and 0.03 m for 2011 and 2012, respectively. For this study, C_b was used as C_bc, because the point clouds generated for this study are not very dense and horizontal point normals were specified, so that a true horizontal retreat rate could be calculated for the Selawik RTS headwall. We found no significant difference when we performed M3C2 analyses with free point normals, because the RTS headwall is vertical. The M3C2 algorithm (Figure 2(C)) returns a text file for each analysis containing the core points used, their normal vectors, the displacement, L_M3C2, for each point in C_bc (Figures 1(B) and 4(A)), the confidence interval and a logical attribute stating whether the displacement is significant or not.

2.3.1. Error Analysis: Spatially Variable Confidence Interval

In conjunction with their change detection algorithm, Lague et al.[22] also proposed a 2σ, spatially variable confidence interval (SVCI) to account for how the local point density and roughness of a point cloud data set affect the precision of the distances computed by M3C2. See Lague et al.[22] for a thorough discussion of the error analysis techniques employed by the algorithm. The local and apparent roughness, the standard deviation of the points in the vicinity of i₁ and i₂, respectively, of the two point clouds being compared can decrease the certainty of differences between the two clouds (Figure 1(B)), leading to a larger local confidence interval [22,31]. For example, if the measured displacement, L_M3C2, between i₁ and i₂ is not very great and the local roughness, σ₁, around i₁ and the apparent roughness, σ₂, around i₂ are very great, the magnitudes of σ₁ and σ₂ could mask the measured displacement between i₁ and i₂ (Figure 1(B)). The SVCI would then consider this displacement insignificant. Surface roughness can also cause local shadowing in point clouds, which also increases surface position uncertainty when comparing point clouds, because a larger search radii must be used to find corresponding points in each point cloud for analysis [22,56]. In addition to roughness considerations, the registration error associated with assigning scan data to the same coordinate system also affects the SVCI computed for every point in E_bc. For this study, the additive RMS error (Table 1) was used as the registration error for M3C2 analyses. Topographic change detection using M3C2 is, therefore, inherently more conservative than C2M, because the M3C2 SVCI estimation accounts for both the RMS error and roughness aspects of the point clouds analyzed, rather than just the RMS error used with C2M.

2.3.2. Application of the M3C2 SVCI

A MATLAB script was used to change M3C2 distance calculations marked as insignificant to zero, because they fall below the SVCI for that point (Figure 4(D)). For each epoch analyzed, the mean of the corrected differences between C_b and C_a was used in conjunction with the elapsed time between C_b and C_a to calculate the retreat rate (m/d) of the RTS headwall. To compare the overall M3C2 level of detection to the C2M detection threshold, all of the SVCI estimates were also extracted from the M3C2 analyses and averaged (Table 1).

3. Results and Discussion

3.1. Comparison of Error Analysis Methods

Using the methods described above, we are able to consistently detect surface changes down to 2–7 cm, depending on the error analysis technique used and whether surface roughness is used to calculate the change detection threshold (Table 1). The average M3C2 change detection thresholds are both larger than the corresponding C2M RMS change detection thresholds. The difference in the change detection thresholds is due to the fact that M3C2 considers the RMS error and the uncertainty due to local and apparent roughness when determining if a calculated displacement is significant. This results in the M3C2 algorithm converting more displacements to zero than C2M (95% vs. 61%, Figures 4 and 5(B)). This difference in the techniques also causes M3C2 to report a larger detection threshold than C2M when the two processes are compared (Table 1). Furthermore, this discrepancy in the techniques allows any increase in the average M3C2 confidence interval over the RMS error to be attributed to uncertainty, due to local and apparent surface roughness in C_b and C_a. In 2011 and 2012, the average uncertainty contributed by roughness is 0.0316 m and 0.0385 m, respectively. This is consistent with Lague et al.[22], who found that most of the uncertainty in direct cloud to cloud comparisons can be attributed to surface roughness, given their very precise scan registration methods. However, at the Selawik RTS, where such precise scan georeferencing and registration techniques are not possible, the contributions of surface roughness and scan registration error to the average SVCI are similar in magnitude. The differences in the uncertainty contributed by roughness in 2011 and 2012 can be attributed to changes in point cloud density between the two data sets, because the scan positions for 2011 and 2012 were nearly identical.

3.2. Comparison of Change Detection Methods

The RTS headwall retreat rate time series derived via M3C2 and C2M are similar in form; however, there are times when the data sets differ significantly. In this portion of the discussion, different portions of the time series in Figure 6 are indexed using line segments for clarity. For the 2011 data set (Figure 6(A)), the C2M and M3C2 analyses show significant differences during segment I and then follow the same pattern through segment II. For the 2012 data set (Figure 6(B)), the C2M and M3C2 analyses have the same general pattern as the 2011 data set, with the two techniques reporting rates of change with different magnitudes in segment III, similar rates of change in segment IV and distinctly different rates of change in segment V. In both 2011 and 2012, C2M consistently reports higher magnitude rates of change than the M3C2 algorithm during periods of comparatively low magnitude change (segments I, III and V) and sporadically reports lower magnitude rates of change during periods when the rate of change is highly variable (segments II and IV). Only rarely does the M3C2 algorithm report greater magnitude rates of change than the C2M technique. These patterns suggest that C2M consistently over-reports the magnitude of the rate of topographic change or that the M3C2 algorithm consistently under-reports the magnitude of the rate of topographic change. It remains unclear which technique best represents the actual topographic change occurring on the RTS headwall.

This knowledge gap could be addressed by building a irregular surface that resembles the RTS headwall surface on a track in a laboratory setting and moving it progressively away from a TLS scanner in known increments. This would simulate a more natural surface, while also removing the uncertainty due to georeferencing. Unfortunately, because of the difficulty in maintaining geodetic control in a permafrost dominated environment, this type of experiment could not be undertaken at this field site, and the RTS headwall could not be physically measured to obtain an actual retreat rate, as it is very difficult and dangerous to access.

The difference between the two data sets likely stems from how the techniques apply their respective change detection thresholds (Figure 4(B,D)) to the raw data produced by the analyses (Figure 4(A,C)). The M3C2 SVCI is clearly more conservative than the blanket C2M change detection threshold. More M3C2 topographic change measurements are converted to zero than C2M topographic change measurements, which shifts the mean M3C2 distance between C_b and C_a closer to zero than the mean C2M measured distance (Figure 4(B,D)). For example, given the displacements calculated between 2012 scans number 3 and 4, the mean raw C2M and M3C2 displacements are −0.023 m and −0.026 m, respectively. This suggests that both methods calculate similar mean displacements before insignificant displacements have been treated. After insignificant displacements have been changed to zeros, the mean C2M and M3C2 displacements are −0.019 m and −0.005 m, respectively. The large difference between these two values is caused because 95% of the calculated displacements in M3C2 are considered insignificant (e.g., Figure 4(D)) by the algorithm, while only 61% of the C2M displacements are considered insignificant (e.g., Figure 4(B)). By considering more values as insignificant, M3C2 converts more displacements to zero than C2M, which results in a lower magnitude mean filtered M3C2 displacement when compared to the mean filtered C2M displacement. This distinction between the two techniques is especially pronounced during periods of low magnitude change, because the M3C2 algorithm interprets more change detection measurements as insignificant, due to the uncertainty added when local and apparent surface roughness are considered. During periods of high magnitude change, this difference is less pronounced, because fewer M3C2 measurements are masked by roughness. How each technique handles insignificant displacements explains why both techniques are better at measuring large changes; however, M3C2 is better at resolving small, significant changes because of its more nuanced treatment of the uncertainty due to registration and roughness. Interestingly, many histograms of raw M3C2 data exhibit the bimodal distribution illustrated in Figure 4(C). While M3C2 has deemed these displacements insignificant (Figure 4(D)), this pattern in the displacement distribution suggests that the mass loss processes through which the RTS headwall retreats control the position of the larger peak, either via large sheet type spallation processes or via the amount of elapsed time between scan epochs.

3.3. Comparison of Cumulative Change

To investigate which change detection method most accurately approximates the cumulative observed change of the headwall, topographic change analyses were conducted for both the C2M and M3C2 techniques using the first and last measurement epochs from the 2012 field season (C₁ and C₄₀). The 2012 data set was used, because it spans a longer record of incremental change than the 2011 data set. The results from this analysis were treated the same as the incremental change analyses used above and are compared to the mean summed retreat of the incremental change analyses above. The summed incremental change for C2M and M3C2 is −1.546 m and −1.003 m, respectively, and the total change for C2M and M3C2 is −1.338 m and −1.629 m, respectively. We expected that both the C2M and M3C2 summed incremental change would be smaller than the the total change between the first and last measurement epochs, because insignificant measurements were removed throughout the incremental analyses. However, this is not the case with the C2M method, because the cumulative measurement path summed across all incremental epochs is longer than the measurement path when only the first and last epochs are considered. It is likely that this problem is present in the M3C2 analyses; however, because the path direction at each increment is specified by the local surface orientation of C_b, the problem is sufficiently minimized.

The difference in summed incremental change between the C2M and M3C2 techniques is due to how each method handles insignificant change. The M3C2 method is inherently more conservative than the C2M algorithm (see above and Section 2.3.1), which causes more displacement measurements to be converted to zeros (Figure 4(D)) than the C2M method (Figure 4(B)). The difference between the total C2M and M3C2 change can be explained by how each technique calculates displacements. C2M calculates a lower magnitude total change, because it calculates displacements based on the closest mesh elements and points in M_b and C_a (Figure 1(A)), while M3C2 uses the local surface orientation of C_b when calculating displacements (Figure 1(B)), which does not bias it towards close points in each topographic change analysis.

The difference between C2M incremental and total change values are less than the differences between M3C2 incremental and total change values (discrepancies of 0.208 m and 0.626 m, respectively). However, because the C2M technique does not take the local surface orientation of C_b into account when computing measurements between measurement epochs, it is biased by point proximity. This bias leads to C2M total change estimates less in magnitude than those computed by the M3C2 algorithm. Since M3C2 takes the local surface orientation of C_b into account when calculating topographic change, we consider it to be a more true measure of actual change than C2M, which is consistent with the findings of Lague et al.[22].

The transparency of the calculations behind the M3C2 algorithm and its ability to calculate a true horizontal retreat rate, by fixing local surface orientation normals to horizontal, make it a well-suited analytical tool to rapidly calculate incremental differences between successive TLS surveys. Additionally, the SVCI used by the M3C2 algorithm takes both the registration uncertainty of C_b and C_a, as well as the local and apparent roughness of the regions being analyzed into account. However, there are scenarios under which C2M may be a better tool for detecting topographic change (Figure 7). For example, M3C2 may be a better topographic change detection tool when the displacement between two surfaces is non-uniform (Figure 7(A)). Under this scenario, because some points in C_a are much closer to C_b, C2M is biased towards calculating displacements between these points. M3C2, because it takes the local orientation of C_b into account, is able to calculate a more true distribution of surface displacements. Conversely, if an undulatory surface changes a relatively large amount in comparison to the wavelength of the undulations, M3C2 may be prone to calculate displacements across undulations (Figure 7(B)). C2M is still biased towards close points in this scenario; however, this bias may be minor compared to the error introduced by M3C2. It should be noted that M3C2 will only behave in this way if D, the scale at which point normals are calculated, is smaller than the undulations in the surface. Increasing D to smooth out these undulations during point normal calculation would correct M3C2 in this scenario.

The authors have found that M3C2 is the most nuanced treatment of TLS point cloud uncertainties for topographic change detection available, because it calculates displacements along local surface orientations and takes into account the most sources of uncertainty when calculating topographic change and the certainty of that change. Conversely, C2M does not provide as nuanced an approach to statistically significant topographic change detection, because it is biased towards calculating displacements between closest points, utilizes only a blanket level of detection based only on systematic errors and ignores local effects, such as local and apparent roughness, on detecting significant topographic change. Furthermore, C2M requires surface models to be created for those data being analyzed, which requires additional processing and interpolation, which makes assumptions about a surface where there are no data. For these reasons, we only discuss the M3C2 results for the remainder of this section. However, it should be noted that implementing a C2M spatially variable error analysis technique similar to the SVCI in M3C2 would not be difficult and has the potential to greatly improve C2M.

3.4. Capturing Spatial and Temporal Variability in RTS Headwall Retreat Rate

The retreat rate of the Selawik RTS derived from the M3C2 analysis results is both spatially (Figure 5(A)) and temporally (Figure 6) variable. The spatial variability decreases, but still persists when only portions of the headwall that have experienced significant change are considered (Figure 5(B)). This suggests that underlying drivers of RTS headwall retreat may cause some of this remaining, statistically significant, spatial variability. It should be noted that most of the RTS headwall does not experience significant change during the 12 h between measurement epochs (Figure 5(B)). Using a data set of this length to examine the temporal variability of RTS headwall retreat is accomplished simply by taking the mean of the significant differences calculated for each analysis, normalizing by the elapsed time between measurement epochs, and treating the resulting data set as a time series (Figure 6). However, fully visualizing and animating a point cloud data set that has been analyzed in this way is problematic computationally, as well as practically, because of the volume and spatial complexity of the information presented. Visualizing spatial subsets of these types of data can highlight both spatial and temporal patterns that recur in a data set, while not overwhelming the analyst. We present two methods: spatial binning and localization studies to investigate these areas.

Spatial subsets of a complex data set can be created by spatially averaging the data created from either the C2M or M3C2 point cloud comparison techniques. For this study, this type of simplification is defined as spatial binning (Figure 3(E)). The RTS headwall is partitioned into 0.5 m-wide vertical sections by UTM easting (∼200 sections per scan), and all of the filtered displacements computed by the M3C2 algorithm in each subset are averaged and normalized by the elapsed time. This technique works well for characterizing lateral variation at the expense of vertical variation. Each spatially binned point cloud can then be plotted as a profile of the RTS headwall viewed in planform (Figure 8(A)). The form of the headwall is represented by dots whose easting and northing correspond to the spatial bin segment centers and the average northing of all the points contained in each spatial bin. The dots are then colored based on the average rate of change calculated for each spatial bin. This way of visualizing complex vertical surfaces simplifies the system just enough to draw interpretations from such spatially and temporally robust data sets. For example, the spatially binned analyses (Figure 8(A)) show the same diel pattern as the time series of RTS headwall retreat rates (Figure 6(A) segment II). This suggests that diurnal hydrometeorological variables, such as insolation or air temperature, are driving the retreat of the RTS headwall. When individual daytime and nighttime change analyses are displayed; the spatial pattern of daytime and nighttime RTS headwall retreat show a possible aspect control on retreat rate (Figure 8(B,C)). Generally, over night change is greatest on west-facing aspects of the headwall (Figure 8(B)), and daytime change is greatest on east-facing aspects of the headwall (Figure 8(C)). This spatiotemporal pattern, previously masked by the spatial complexity of this data set, suggests that insolation, over air temperature, may be a key driver of RTS headwall retreat.

The spatial variability of RTS headwall retreat (Figure 5) also changes through time. To investigate this, 11 local portions of the RTS headwall (Figure 3(D)) were isolated and analyzed with the M3C2 algorithm to create additional time series of the RTS headwall retreat rate (gray lines in Figure 6). We define these analyses as localization studies. These analyses suggest that, at times, small portions of the headwall accurately represent the averaged retreat rate of the entire headwall. While no one locality on the headwall always approximates the retreat of the entire headwall, when averaged together, all of the localization studies closely represent the retreat of the entire headwall (thick black line in Figure 6). This suggests that average RTS retreat rates could be accurately measured with an array of laser rangefinders or a laser profiler instead of a TLS. Plots, such as Figure 6, also provide information about when the RTS headwall is behaving uniformly and when it is highly spatially variable. Epochs with high spatial variability (arrows in Figure 6(B)) coincide with misty and rainy meteorologic conditions, which suggests that these factors may play a role in the mass loss processes of the RTS headwall or that they affect the ability of the TLS instrument to accurately acquire topographic data or both.

3.5. Unconstrained Sources of Uncertainty

Equation (1) accounts for the readily quantifiable sources of error in the scanning and georeferencing system; however, there are other sources of uncertainty that may also decrease the level of topographic change detectable at the Selawik RTS. The heterogeneous surface properties of the RTS headwall could create a complex spatial noise pattern, depending on the spatial distribution of different surface types (e.g., gravel, silt, cobbles, wet, dry, peat or ice). Additionally, the presence of rain or mist in between the instrument and the headwall could cause small systematic shifts in the position of the headwall. Rain and mist could also wet the headwall surface via rainfall or increased tundra surface overflow, which could also wet the headwall. Wet areas of the headwall were observed to have fewer laser returns than dry surfaces, and returns from wet surfaces could possibly have increased uncertainty. The combination of these unaccounted factors could contribute both to an increase in the spatial variability of RTS headwall retreat, either because of an increase in spatially discontinuous mass loss processes or an increase in the level of noise in the change detection system. Furthermore, although an atmospheric correction was conducted prior to each scan, the long duration of the scans used for this study, approximately 30 min each, coupled with the rapidly changing meteorologic conditions often experienced at the field site, could cause measurement epochs where portions of the headwall were recorded under relatively clear atmospheric conditions, while other portions, within the same measurement epoch, were recorded with mist and rain partially obscuring the feature. Not only are these considerations not described by Equation (1), but the magnitude of their potential effects also varies in space and time as the headwall erodes and weather patterns shift at the field site, further complicating the tabulation of errors associated with point cloud measurements.

4. Conclusions

This study presents a comparison of two point cloud change detection techniques, cloud to mesh (C2M) and Multiscale Model to Model Cloud Comparison (M3C2) [22], which were applied to an extensive TLS data set gathered under difficult survey conditions at a permafrost dominated field site. We also leverage this extensive TLS data set, gathered at the Selawik retrogressive thaw slump (RTS), to compare error analysis methods, as well as to investigate techniques to visualize the spatial and temporal variability of RTS headwall retreat.

Two-to-seven centimeter topographic changes can be resolved under challenging survey conditions; however, daily resurveying of TLS georeferencing targets and two-sigma change detection thresholds were needed to reach this level of precision. We found M3C2 to be the most robust point cloud analysis method available, because it integrates a spatially variable confidence interval (SVCI) as a change detection threshold with displacement calculation. Furthermore, M3C2 can calculate true horizontal displacements, which allow the extraction of horizontal erosion rates from point cloud comparisons. Conversely, C2M uses a blanket change detection threshold and cannot calculate true horizontal displacements. A comparison of first and last measurement epochs using both methods illustrates that C2M reports a higher magnitude of topographic change over short periods of time (∼12 h) and reports a lower magnitude of topographic change over long periods of time (∼four weeks) when compared to M3C2. The SVCI used with M3C2 may classify too many displacements as insignificant, due to the uncertainty caused by surface roughness and scan registration; however, because we lacked true measurements of RTS headwall retreat, we were unable to evaluate this.

The spatial variability of the Selawik RTS headwall retreat changes through time, as some measurement epochs show more spatial variability of retreat than others. However, when averaged together, localized measurements of headwall retreat approximate the average retreat rate of the entire headwall, which suggests that low point density laser profilers could be used to monitor the retreat of RTS type landforms. This study also employed a spatial binning technique to visualize both temporal and spatial changes in the RTS headwall retreat rate. Visualizing data in this way leads to insights, such as a possible aspect and insolation control on headwall retreat (Figure 8(B,C)).

Although we employed a spatially variable change detection threshold in this study, the systematic errors associated with a TLS data set of this temporal length also vary in time. Future studies of this kind could employ temporally variable additive RMS estimates, which could provide a more nuanced error analysis treatment than done in this study. This would likely decrease the overall error observed in the change detection analyses and would limit the error introduced by uncertain measurement epochs to measurements before and after the uncertain epoch. Furthermore, the time between measurement epochs could be increased to reduce the noise to displacement ratio for each measurement epoch. Now that automated point cloud classification and analysis algorithms are available, as well as TLS instruments with greater laser pulse energy and sensor sensitivity, higher density point clouds could be generated from RTS features, and the benefits of spatially averaging points could be used to reduce noise in the point cloud data set. Additional work to better understand how adverse atmospheric conditions affect point cloud precision, as well as how heterogeneous and wet surfaces contribute to noise on natural surfaces represented by point clouds could also contribute to our knowledge of the unconstrained sources of uncertainty mentioned in this work.