Optimising 4D Surface Change Detection: An Approach for Capturing Rockfall Magnitude-Frequency

. We present a monitoring technique tailored to analysing change from near-continuously collected, high-resolution 3D data. Our aim is to fully characterise geomorphological change typified by an event magnitude frequency relationship that adheres to an inverse power law or similar. While recent advances in monitoring have enabled changes in volume across more than seven orders of magnitude to be captured, event frequency is commonly assumed to be interchangeable with the 10 time-averaged event numbers between successive surveys. Where events coincide, or coalesce, or where the mechanisms driving change are not spatially independent, apparent event frequency must be partially determined by survey interval. The data reported has been obtained from a permanently installed terrestrial laser scanner, which permits an increased frequency of surveys. Surveying from a single position raises challenges, given the single viewpoint onto a complex surface and the need for computational efficiency associated with handling a large time series of 3D data. A workflow is presented 15 that optimises the detection of change by filtering and aligning scans to improve repeatability. An adaptation of the M3C2 algorithm is used to detect 3D change, to overcome data inconsistencies between scans. Individual rockfall geometries are then extracted and the associated volumetric errors modelled. The utility of this approach is demonstrated using a dataset of ~ 9 × 10 3 surveys acquired at ~ 1 hour intervals over 10 months. The magnitude-frequency distribution of rockfall volumes generated is shown to be sensitive to monitoring frequency. Using a 1 h interval between surveys, rather than 30 days, the 20 volume contribution from small (< 0.1 m 3 ) rockfall increases from 67% to 98% of the total, and the number of individual rockfall observed increases by over three orders of magnitude. High frequency monitoring therefore holds considerable implications for magnitude-frequency derivatives, such as hazard return intervals and erosion rates. As such, while high frequency monitoring has potential to describe short-term controls on geomorphological change and more realistic magnitude-frequency relationships, the assessment of longer-term erosion rates may be more suited to less frequent data 25 collection with lower accumulative errors.


Size distribution of geomorphic events
The processes that erode landscapes involve a broad range of event sizes, the distribution of which is commonly characterised using magnitude-frequency curves.Wolman and Miller (1960) proposed that the frequency of events that denude the Earth's surface is log-normally distributed, and that their geomorphic effectiveness (the product of magni-tude and frequency) is greatest for the frequent, moderately sized events.This concept has been widely applied to study both the geomorphic efficacy of rivers (Wolman and Gerson, 1978;Hooke, 1980;Nash, 1994;Gintz et al., 1996) and the characteristics of landslides (Hovius et al., 1997(Hovius et al., , 2000;;Dussauge-Peisser et al., 2002;Turcotte et al., 2002;Dussauge et al., 2003;Malamud et al., 2004;Guthrie and Evans, 2007;Li et al., 2016) using inverse power-law distributions or similar.
Published by Copernicus Publications on behalf of the European Geosciences Union.
The exponent of the inverse power law describes the proportional contribution of increasingly small events, with larger exponents representing an increase in the proportion of small events in the inventory.However, many landslide volume distributions have been characterised by a decrease in the frequency density of the smallest events in log magnitude-log frequency space, known as a "rollover" (Malamud et al., 2004).At this point, the inverse power law breaks down, and so alternative distributions such as the double Pareto (Stark and Hovius, 2001;Guzzetti et al., 2002) or inverse gamma (Malamud et al., 2004;Guzzetti et al., 2005) have been drawn upon to model observations.Explanations for this rollover have been widely considered and include mechanical differences and physically based minimum possible event sizes (Pelletier et al., 1997;Guzzetti et al., 2002;Guthrie and Evans, 2004) or censoring of the smallest events by the resolution or frequency of monitoring (Lim et al., 2010).For rockfalls, Malamud et al. (2004) hypothesised that a rollover may not occur due to rock mass fragmentation.

Temporal resolution of monitoring
The duration of monitoring relative to the return period of all possible event sizes determines the likelihood of detecting changes that are representative of how a landform evolves over longer timescales.The creation of a representative inventory is also a function of the smallest event size that can be detected and the temporal frequency of monitoring compared to the rate at which such small events occur.Abellán et al. (2014) suggested that the spatial resolution of rockfall monitoring should be sufficient to discretise the smallest events in a magnitude-frequency distribution and that the recording frequency should fall below the timescale on which superimposition and coalescence may occur.In practice, defining this timescale a priori is challenging and requires the ability to monitor the rock face over a sustained period in (near) real time.For rockfalls, high-resolution monitoring also shows evolution of failures through time, with event sequences and patterns related to the incremental growth of scars (Rosser et al., 2007(Rosser et al., , 2013;;Stock et al., 2012;Kromer et al., 2015a;Rohmer and Dewez, 2015;Royán et al., 2015).Barlow et al. (2012) showed that a monitoring interval of 19 months underestimated the frequency distribution of small rockfall events, which coalesced into or were superimposed by larger rockfalls.Treating rockfalls as spatially and temporally independent is therefore problematic, as is experienced in other types of landform change.For example, Milan et al. (2007) found an increase in erosion and deposition volumes within a proglacial river channel when monitored using daily terrestrial laser scan (TLS) surveys as opposed to surveys separated by 8 days.This was attributed to the temporal length scales of discharge and sediment supply, with return periods of less than 8 days.The influence of monitoring or sampling interval on measured process rates is more widely considered in the Sadler effect, in which sediment accumu-lation rates observed in stratigraphic sections exhibit a negative power-law dependence on measurement interval (Sadler, 1981;Wilkinson, 2015).Importantly, therefore, in settings that change little but often, the ability to capture true magnitude and frequency without higher-frequency monitoring is subject to an unknown degree of event superimposition and coalescence, as well as temporal coincidence.Our first aim, therefore, is to capture the influence of near-continuous monitoring on the magnitude-frequency distribution of rockfalls from an actively failing rock face, across which spatially contiguous rockfalls have previously been observed (Rosser et al., 2005).This draws upon a database of ∼ 10 3 individual 3-D scans collected at ∼ 1 h intervals, with each comprising > 10 6 points.

Scanning from a fixed position
The improvement in temporal resolution gained by monitoring from (semi-)permanent installations is weighed against a series of compromises in the quality of data generated.These ultimately arise from scanning a complex surface from only a single position, which has not previously required consideration due to the scarcity of near-continuous lidar monitoring in the geosciences.Scanning from a single position has a direct impact on the similarity of point distributions between successive surveys, which is critical in determining the accuracy of subsequent change detections.Scanning from a single position results in occlusion, leaving "holes" in the final point cloud in areas invisible to the scanner.At topographic edges, range measurements may be averaged from multiple returns recorded from separate surfaces, which are intersected within a single footprint.Since laser scanners never measure exactly the same point twice (Hodge et al., 2009); the perimeter of these holes and the position of topographic edges will move between successive point clouds, even despite no movement of the instrument.The likelihood of generating similar point distributions between surveys is also influenced by surface reflectance characteristics, such as moisture and colour, and by surface relief (Clark and Robson, 2004;Bae et al., 2005;Lichti et al., 2007;Kaasalainen et al., 2008Kaasalainen et al., , 2010;;Pesci et al., 2008Pesci et al., , 2011;;Soudarissanane et al., 2011).Scan lines in most laser scanners result in non-uniformly distributed data, with heterogeneity often on a scale and orientation comparable to surface structure, leading to aliasing that is also inconsistent (Lichti and Jamtsho, 2006).The influence of these combined effects is exaggerated if the scanner view is oblique to the surface, which may be common due to logistical constraints when siting a semi-permanent instrument.Our second aim, therefore, is to minimise the errors that arise from near-continuous monitoring in order to reduce the minimum detectable movement.Kromer et al. (2015b) presented a 4-D smoothing technique to reduce the offset between successive point clouds, such as those from near-continuous monitoring.Similarly, the method presented here is optimised for handling large (10 2 -10 4 ) numbers of high-resolution 3-D scans, critically without user intervention, but can also be applied to point clouds from non-continuous monitoring.

Uncertainty in near-continuous monitoring
The minimum detectable movement, or level of detection (LoD), is a fundamental parameter in the delineation and calculation of erosion volumes, here rockfalls.This involves masking regions of change that exceed a hard threshold at the LoD, which is estimated either locally (e.g.Wheaton et al., 2010;Lague et al., 2013) or across the entire point cloud (e.g.Rosser et al., 2005;Abellán et al., 2009).Methods that estimate spatially variable LoDs have enhanced the ability to identify volumetric loss as compared to the application of a single LoD, with the latter set to exceed a significant portion of the modelled uncertainty across the area of interest (AOI).The spatially variable uncertainties described in Sect.1.3 raise the potential for real change to be masked when using a single LoD but, equally, the benefits of using a single LoD are primarily in the consistency in measurement across the AOI.For example, if the purpose of monitoring is to generate a rockfall inventory in which the relative magnitude of events is important, a single LoD ensures consistency in the minimum detectable rockfall across the AOI and minimises the potential for recording erroneous events.The application of a single LoD also becomes increasingly computationally efficient when dealing with a large number of surveys.While the LoD remains constant, however, the accuracy of volume estimates is also contingent upon the ability to accurately identify depth change within boundary cells that occupy the rockfall perimeter.Previous approaches have ignored cells with a depth change below the instrument precision, thereby assuming that erosion events with an aerial extent below the cell size cannot be detected (e.g.Dussauge et al., 2003;Rosser et al., 2005;Abellán et al., 2006).However, these often fail to model volumetric errors that arise from extrapolating measured depth changes across each cell at the rockfall perimeter.In accounting for this, the volumetric uncertainty increases as a proportion of estimated volume for smaller events and specific geometries.Our third and final aim, therefore, is to describe the impact that a changing magnitude-frequency distribution has on the overall uncertainty of eroded volume estimates through time.

Study site and data collection
Data are presented from a monitored coastal cliff in North Yorkshire, UK.Located at East Cliff, Whitby, the rock cliff is near vertical, reaching ∼ 60 m in height, and is actively eroding.The erosion of this coast has previously been monitored and averages ∼ 0.1 m a −1 (Rosser et al., 2005(Rosser et al., , 2007;;Miller, 2007).Typical rockfalls include small-scale joint defined wedges and larger-scale failures released via rock bridge breakage, which can be inferred from the exposed fresh fracture surfaces visible after failure.Rockfalls have been mea- sured up to 2.5 × 10 3 m 3 , but the volume loss is dominated by smaller-scale rockfalls with median volumes approaching 1.0 × 10 −3 m 3 (Rosser et al., 2013).The slope monitoring system presented surveys the cliff using a remotely controlled Riegl VZ-1000 laser scanner, housed inside the former lantern room at the top of East Pier lighthouse (Fig. 1), ∼ 350 m seaward of the cliff face (Fig. 2a).
The viewpoint of the scanner results in some loss of spatial continuity in surface measurement due to occlusion, as a result of surface relief and the high incidence angle of parts of the rock face relative to the scanner (Fig. 2b).The closest point on the cliff is 342 m from the scanner with an incidence angle onto the strike of the cliff of ∼ 25 • , while the furthest monitored point is 533 m at an angle of ∼ 42 • .Range correction for atmospheric effects and precise point cloud alignment is automatically conducted every 3 h using veryhigh-resolution scanning (5.0 × 10 −4 m point spacing) of six fixed 0.25 m 2 control targets, the precise relative positions (±5.0 × 10 −3 m) of which are known from a total station survey.Atmospheric range correction, derived from comparing scanned to surveyed target ranges, typically scales range measurements by ∼ 1.0 × 10 −5 with minimal diurnal fluctuation and is therefore largely inconsequential at this site, but would be significant for locations with greater ranges or areas subject to more extreme atmospheric conditions.The TLS survey was managed using SiteMonitor4D (3D Laser Mapping Ltd), which schedules scans, manages the atmospheric correction, and applies an affine rigid-body rotation matrix to compensate for tilt and yaw in the scanner position based upon the scanned control target positions in real time.The reported dataset has been collected using this setup between 5 March and 31 December 2015, totalling ∼ 9000 scans with ∼ 1 h intervals between surveys.Gaps in the dataset arise from system outages and are excluded from the analysis.

Method: optimising event extraction from near-continuously collected point clouds
Repeatability in change detection is dependent on point clouds that consistently describe the monitored surface.
While the specification of successive scans is identical here, and despite no positional change in either the instrument or the surface, individual point clouds differ due to the inherent uncertainties described in Sect.1.3.Minimising the impact of these differences on the resulting rockfall inventories is the primary objective of the preprocessing phases of our method, involving point cloud filtering and alignment, as well as the change detection phase.Once 3-D change has been calculated, the point clouds of change are interpolated and classified into 2.5-D datasets from which 3-D change geometry can be analysed (Fig. 3).While this paper does not seek to create a real-time system, the approach described is optimised for computational efficiency to allow data to be processed at a rate that is at least as quick as collection.

Point cloud filtering
Optimising a point cloud for change detection involves filtering points on features that cannot be consistently measured, such as edges and vegetation.For large numbers of scans obtained from near-continuous monitoring, an automated approach is required that finds and removes such points consistently for each cloud.Below we describe the application of filtering prior to change detection.Given that the change detection process draws on a neighbourhood of points to describe the surface morphology, filtering points in advance removes the need for change detection of erroneous measurements and ensures that such measurements do not impact the change detection of surrounding neighbours.

AOI extraction
Scan data collected outside of the AOI can provide useful information for scan-to-scan registration, particularly if they cover a wider geographical area, have a distribution that is less planar than the AOI itself, or if a large portion of the AOI is undergoing deformation.Here, given that the control target network provides range estimate control and the scanner remains unmoved, the spatial extent of the AOI is clipped at the beginning of the workflow in order to increase the speed of subsequent processing steps.For a small number of scans (< 20), AOI extraction can be undertaken manually.However, for the ∼ 9000 scans used, points that lie outside a predefined cuboidal bounding box are automatically removed, which typically reduces the raw dataset from ∼ 1.9 × 10 6 to 1.1 × 10 6 points.

Edge and hole filter
Once the AOI has been extracted, morphological (this section) and radiometric (Sect.3.1.3)filters are applied to remove points with high positional uncertainties.These specific filters have not previously been used to our knowledge; however, the reasoning behind their implementation is well documented as described in Sect.1.3.We refer to edges as morphological breaks in gradient on the slope surface and holes as features that surround regions of occlusion in the point cloud, which frequently occur at edges due to high surface inclination.To detect edges, neighbouring points within a fixed radius of each (query) point, q, are identified and the central position of the neighbourhood points, CoG, calculated.The 3-D Euclidean distance, ED, between q and the CoG is then calculated.For a point at an edge, CoG tends away from the query point.The ED is therefore larger for query points that lie closer to an edge.Applying a threshold to ED needs to account for the varying point density across the cloud as in regions of low point density ED will always be larger.The edge-hole value EH assigned to each point is therefore reported as a ratio of the ED to the number of points in a spherical domain centred on each point: where k is the number of neighbouring points.

Waveform deviation filter
A limitation of many laser scanners is the inability to quantify the accuracy of each range measurement.With the common absence of more accurate data, assessing reliability is therefore challenging.In most TLS systems used for rock slope monitoring, range is estimated using the time of flight of a laser, where time is stamped based upon a characteristic of the measured reflection (e.g.intensity gate, maximum intensity amplitude) that varies between scanners.Some systems have the ability to capture the full energy-time distribution of the reflection.The energy of the received laser pulse structure depends on the spatial and temporal energy distribution of the emitted pulse, which are modified by the geometric and reflectance properties of the target surface (Stilla and Jutzi, 2008;Soudarissanane et al., 2011;Hartzell et al., 2015;Telling et al., 2017).This provides a means of estimating the relative quality of recorded measurements as a function of the number of separate reflections from a single pulse, the incidence angle of the laser beam with respect to the surface (the elongation through time of the reflected pulse relative to the emitted pulse), or the reflectance intensity (the integral of the reflection energy-time distribution).Here, the characteristics of the returned signal are used to remove vegetation (multiple returns per pulse) and edges (elongated reflections) in order to increase the consistency between successive point clouds.While the sensitivity of the waveform to target geometry has previously been highlighted (Williams et al., 2013), it has not previously been documented as a method to filter points acquired from terrestrial lidar.The Riegl VZ-1000 TLS, with "full waveform" capture, records the intensity of each returned signal at 2.01 × 10 −9 s intervals, providing 15-70 amplitude measurements per pulse.The "deviation", δ, of the waveform describes the change in shape of the received waveform relative to a modelled (emitted) Gaussian energy-time distribution according to where N is the observations in pulse s i , compared to the reference values, p i .δ = 0 represents identical emitted and received waveforms, as would be expected from a nadiroriented planar specular surface.δ is less sensitive to target range than incidence angle.

Filtering of partially obscured point clouds
In addition to the filtering of individual points considered unreliable, entire surveys required removal from the overall inventory due to inclement weather conditions.In conventional monitoring with no scan schedule automation, scans that are partially or fully obscured due to rain or fog are manually removed and/or repeated.Due to the frequency and duration of near-continuous scanning, some scans will almost certainly be partially or fully obscured.Scans that are entirely obscured can be removed automatically with relative ease.Unobscured areas in partial scans still allow some accurate change detection and rockfall identification, and so may be valuable to retain.Given that change is detected between scan pairs, however, it is critical that these partial scans are removed prior to change detection and remain unused.12:00 and 12:30 (all times are GMT) during adverse weather, which partially obscures the impending scan at 12:30.While some areas of this scan allow accurate change detection of the surface, if the rockfall occurs in an obscured area, it can be omitted from the inventory entirely.However, if surfaces are compared between 12:00 and the following scan at 13:00, with both captured during fair conditions, the rockfall will be observed and included in the inventory.
Given the variability in the persistence of poor weather during a scan, a threshold based on the number of points can be difficult to define.At present, no automated method for detecting partial scans has been developed, though the removal of any scan that coincides with measured rainfall may represent a first step.Here, given that the same scan line patterns are applied during each survey, partially obscured point clouds were identified as those > 1 MB below the average file size.While the maximum possible number of change detections was 8986, this was reduced to 8596 because of complete obscuration and finally to 8270 because of partial scan removal.The point distribution of the remaining scans was manually examined by creating a video of every point cloud prior to reanalysis of the dataset.The reduction in the number of scans has a direct impact on the time interval between scans and hence deformation analysis prior to failures that occur during bad weather, which may result in some of the most active periods of rockfalls.

Precise alignment
While range correction factors automatically scale range estimates in response to atmospheric variation, point clouds collected from a single fixed position still require alignment due to small shifts in scanner inclination that have the potential to propagate to several centimetres over distances of several hundred metres.It is therefore assumed that successive point clouds are approximately but not perfectly aligned, and hence require adjustment.In this section we describe the protocol for scan alignment applied to the nearcontinuous dataset.As discussed by Abellán et al. (2014), aligning point clouds can be undertaken using common surveyed and modelled targets combined with measured global coordinates (Teza et al., 2007;Olsen et al., 2009), featurebased registration based on the planarity and curvature of surfaces (e.g.Besl and Jain, 1988;Belton and Lichti, 2006;Rabbani et al., 2006), and point-to-point and point-to-surface methods, which use iterative closest point (ICP) alignment to progressively reduce the distance between two clouds (Besl and McKay, 1992;Chen and Medioni, 1992;Zhang, 1994).The accuracy of alignment is one of the key sources of error when detecting change between two point clouds (Teza et al., 2007).
Here, once filtered, point clouds are automatically registered to a reference point cloud to improve alignment.Registration of datasets into a global system has a significant impact on data file sizes due to multi-digit coordinates, which becomes problematic with large numbers of large point clouds from near-continuous monitoring.By retaining a local coordinate system in this study, resulting file sizes were halved.ICP registration was applied using MATLAB ® 's pcregrigid function (Besl and McKay, 1992;Chen and Medioni, 1992), which searches for the closest point in the reference scan for each point in the moving scan and estimates the combination of rigid rotation and translation that best aligns them (Mitra et al., 2004).Pottmann and Hofer (2003) showed that for point clouds that are approximately aligned, as here, minimising the point-to-plane distance provided the best estimate of convergence.Given that the raw point cloud data commonly have systematic structure (vertical or horizontal scan lines), and where the AOI is approximately planar (a ∼ 2-D rock face rather than a fully ∼ 3-D scene), the success of the ICP can be dominated by aligning structure in the data, rather than the macro-scale cliff geometry.Point-to-point minimisation in ICP was therefore found to be less effective than point-to-plane alignment in this instance, using point clouds down-sampled to a fixed 0.25 m point spacing.
There is no established protocol for choosing the ideal reference scan for alignment.This reference scan may be the first available scan, a later single scan, or an average of a subset of previous scans.Schürch et al. (2011) aligned a series of three scans to the previous scan in a sequence, rather than to the first of the monitoring campaign, in order to gain more precise change estimates between successive surveys, at the cost of the overall positional accuracy.This procedure is advantageous as it ensures that the shape of a rapidly deforming or changing surface can be matched to the previous survey, rather than one captured considerably earlier.For the nearcontinuous monitoring dataset, however, the series of scans is aligned to the first survey.This is undertaken because ICP alignment minimises the point-to-plane distance of downsampled point clouds, such that individual (small relative to the AOI) rockfall events do not impact upon the overall success of the alignment.Second, even with low alignment errors between scan pairs, the potential for the point clouds to drift over time increases with the number of scans sequentially aligned.While Schürch et al. (2011) assessed pairwise change between scans, high-frequency data allow for change detection to be conducted over multiple intervals that the time series enables, so aligning all scans with respect to each other is important.Third, given that all scans were aligned to the initial reference scan, these could be assigned to a single hierarchical structure created only for the reference scan.The partitioning of 3-D data, here into an octree structure, enables faster searching of point clouds with each point assigned a 3 × n bit code, where n is the maximum octree level (Frisken and Perry, 2002;Girardeau-Montaut et al., 2005;Jaboyedoff et al., 2007Jaboyedoff et al., , 2009;;Elsberg et al., 2011Elsberg et al., , 2013;;Hornung et al., 2013).Partitioning point clouds into a single predefined structure reduces computation time and ensures consistency in subsequent operations.For neighbourhood searches, such as during normal vector estimation, points from both the individual octree cube and the surrounding 26 cubes are used.As highlighted by Girardeau-Montaut et al. (2005), the subdivision level at which normal estimation and change detection is performed therefore influences only the computation time.

Normal vector estimation
The distance between successive clouds is measured along the normal vector of each point in the cloud.Accurate estimation of each normal vector is therefore critical in determining the magnitude and direction of change and should be derived from an appropriately sized neighbourhood of points that adds topological context (Riquelme et al., 2014).In this section, we describe the widely adopted method for normal estimation and gains in efficiency that have been made in this study, including the definition of a single reference map of optimal neighbourhood radii.
In order to calculate the normal direction of each neighbourhood, a tangent plane must be fitted to every point and its neighbours, with each being considered as a potential plane subset.Using eigenvectors calculated from principal component analysis, the eigenvector v 3 with the smallest associated eigenvalue is orthogonal to the plane, and therefore defines the normal (Hoppe et al., 1992).It so follows that the plane minimises the sum of squared distances to the neighbours of query point, p, and passes through the centroid, p, where r is the number of neighbours in the neighbourhood, and p i represents the Cartesian coordinates of each point within the neighbourhood (Pauly et al., 2002).The identification of a local surface normal using the third eigenvector is the equivalent of forming a total least-squares fitting plane.However, in a total least-squares fitting the entities in the covariance matrix are not divided by k, and the smallest eigenvalue is equal to the sum of the residuals squared (Pauly et al., 2002;Belton and Lichti, 2006).
The sign ambiguity of each vector is also corrected (Mitra and Nguyen, 2003;Ioannou et al., 2012).This can typically be resolved using the position of the query point, q, relative to the sensor position, s, by (5) In where × denotes the vector cross product and denotes the Euclidean norm of the cross product.α denotes the angle between the unit-normal vector n at q and the vector between q and s, ŝ.If α > π 2 or α < − π 2 , i.e. if the angle between the direction of the normal vector and the vector between the surface and the sensor is not within ±90 • , the normal direction n is reversed: In order to minimise the computation time required to apply Eqs. ( 5) and ( 6), the axis orthogonal to the surface is introduced as a string of either "X" or "Y ", similar to other point cloud processing packages such as MATLAB ® and CloudCompare.With this information, the relevant component of the unit vector is used to determine whether the vector should be reversed or not.For example, if the approximate range in a near-nadir point cloud is measured along the y axis and the vector n = −u, −v, −w then n is directed into the surface and should be reversed using Eq. ( 7).For each normal, this can provide a ∼ 40-50 % reduction in the time taken for sign correction relative to Eqs. ( 5) and ( 6).
The neighbourhood size strongly determines the direction of surface normals (Mitra and Nguyen, 2003;Lalonde et al., 2005;Bae et al., 2009;Lague et al., 2013;Riquelme et al., 2014).If the size of the neighbourhood is below the scale of surface roughness, the resulting normals will fluctuate in direction and are less likely to be consistent between successive point clouds.Here, by varying the size of the neighbourhood for each point between 0.1 and 2.5 m, the radius that produced the most planar surface is identified, which is ideally suited to normal estimation (Lague et al., 2013).An example of this is shown in Fig. 5a, with Fig. 5b illustrating surface planarity across East Cliff.This shows a clear similarity to the distribution of point density, such that the search radius is increased in regions of low point density.Importantly, identifying the optimum neighbourhood radius for 10 3 -10 4 point clouds adds considerable computational cost in processing.As a compromise, the neighbourhood radius of each point in   this study is made equal to the distance to the closest point in the reference cloud in Fig. 5a.Notably, the normal for each point estimated uses the second cloud, such that change is accurately measured along the normal of a planar, post-failure surface, rather than the yet-to-fail surface (Fig. 6).

Change detection
The distance calculation used is based upon the structure of the M3C2 algorithm, developed by Lague et al. (2013).
The algorithm is described first followed by a modification, which has been incorporated to improve the overall accuracy of change detection and to streamline the workflow when applied to large time series scan datasets.
Once the normal vector is estimated, a bounding cylinder with a user-defined radius is created along the normal running through the query point.In order to enforce the boundaries of this cylinder, the orthogonal distance between every point within the current and neighbouring 26 octree cubes and the normal vector was estimated: where d is a vector that connects each neighbour point p to a point on the normal vector n, such as the query point, q.The projection of each point onto the normal P is therefore and the orthogonal distance is Given that the position of each neighbouring point and its orthogonal distance to the normal vector are known, the cylinder boundaries can be enforced using the user-defined cylinder radius, r, retaining only points at which d orth ≤ r (Fig. 7).Once the points, c, in the cylinder are isolated for both point clouds, the mean point CoG is estimated by where k c is the number of points in c.Both mean points are then projected onto the normal vector using Eqs.( 9) and ( 10).
The mean projected points of each sub-cloud, CP, are subtracted to give a distance vector, v: Earth Surf.Dynam., 6, 101-119, 2018 www.earth-surf-dynam.net/6/101/2018/If the vector of change is along the direction of the normal vector (forward movement), the dot product of both vectors is > 0. If the vector of change is counter to the normal direction (backward movement), the dot product is < 0 and the vector is inverted.
M3C2 imposes a user-defined maximum cylinder length to decrease processing times.Cylinder length is critically important for determining the accuracy of change estimation, particularly at topographic edges within the point cloud.As described, edges are likely to be more prevalent in point clouds collected from single or off-nadir viewpoints.A method to reduce the effect of edge change uncertainty in change detection is therefore required.In Fig. 8, the influence of the choice of cylinder length is illustrated with respect to a jointed rock mass surface.The plots illustrate variation in measured change for a single point.When the cylinder extends 0.25 m in both directions, only points from this surface are included in the cylinder; as such, the centroid positions of each point cloud are both fitted onto that surface.The measured distance for this point, the distance between the two centroids, is +0.0011 m.With a cylinder extending ±0.50 m, points that lie between surfaces are included in the change detection.Given that the distribution of points is rarely consistent between point clouds, the position of the centroid of each neighbourhood differs considerably from the centroids estimated using a shorter cylinder and the resulting change estimate is −0.1460 m.At a length of ±10 m, the cylinder intersects multiple surfaces and the centroid positions are averaged between these surfaces.The inclusion of a greater number of points over a wider area increases the similarity of the mean position in both point clouds, but the resulting vector of change is +0.0938 m, a difference of 0.24 m from the 0.50 m cylinder length and significantly higher than the true change estimate.To address this, a distance along the normal with variable cylinder length (DAN VCL) for each point is used.The approach begins with a cylinder that extends ±0.10 m.If fewer than four points are found, the minimum number to estimate a centroid, the cylinder extends.This process is recursive and accepts a user-defined range of cylinder lengths.

Extracting discrete changes and quantifying volumetric error
The use of an increased number of scans over a consistent timescale increases the potential for error propagation and accumulation within the dataset.Here we aim to model the uncertainty of eroded volume estimates by focussing on the conversion of rockfall depth and area into volume for cells at the perimeter of a rockfall scar.Given that the proportion of rockfall cells at the edge of a scar decreases for a larger footprint, the uncertainty in rockfall volume estimates has the potential to vary with the temporal resolution of monitoring.An LoD was identified between scan pairs in which no rockfalls occurred as 2 standard deviations of the 3-D change, after Abellán et al. (2009).This was of comparable magnitude to the LoD recorded for every scan pair in the dataset; hence, the maximum recorded LoD was applied to all scan pairs in the dataset.Similar to Kromer et al. (2017), these change estimates are assumed to include the registration error, which is reduced here through range correction using finely scanned targets and through ICP.For sites whose geometry creates a highly variable point spacing within a single survey, a spatially variable LoD would be more appropriate even for the purpose of compiling an inventory of geomorphic events, so long as a record of the LoD across www.earth-surf-dynam.net/6/101/2018/Earth Surf.Dynam., 6, 101-119, 2018 the surface is kept.Open-pit highwalls, for example, typically comprise a series of benches to minimise the travel distance of rockfalls downslope.This design generates considerable variation in instrument-object distances across the slope and a spatially variable LoD.More broadly, spatially variable LoDs can be considered better suited to measuring total erosion budgets across a single surface than the relative contribution of individual events of varying sizes.
The resulting LoD was used to threshold 2.5-D rasters of the 3-D change data, created by linear interpolation of change values across the x − z plane.The images produced included consistently located holes (no data) due to occlusion, which were identified and masked for each raster.Pixels that consistently exceeded the LoD within the first 100 point clouds, including in a non-systematic manner (e.g. both forward and backward movement), were also masked from all rasters.This prevented several (predominantly single) pixels of noise from being identified as detachments.
Once a change image is thresholded according to the LoD, the volume of each erosion event, V E , can be calculated as where N is the number of cells that are classified as volume lost, d i is the depth of change in cell i, and A C is the cell area.Previous approaches have ignored cells with a depth change below the instrument precision and assumed that erosion events with an aerial extent < A C cannot be detected but often fail to quantify uncertainty in volume estimates derived using Eq. ( 14) for rockfalls with areas greater than A C (e.g.Dussauge et al., 2003;Rosser et al., 2005;Abellán et al., 2006).Basic assumptions about how uncertainty in aerial extent propagates into volumetric uncertainty are needed, in particular for failures of varying geometry.This is of critical importance considering the relatively low spatial resolution of raster cells (here 0.15 m) relative to the accuracy of the change in depth within pixels recorded by TLS (here 1 in 10 000 to 1 in 100 000).Assuming any cell that lies on the boundary of an area of change can contain any fraction (0-1) of true change, the maximum area of change A E max is In reality, Eq. ( 15) represents the largest possible area because the likelihood that border cells are entirely covered by the true change is small.Conversely, the theoretical minimum area A E min approaches where N b is the number of boundary cells.The maximum range in uncertainty associated with the area of change is then This value can be applied as a threshold to the rockfall inventory, such that failure areas below A maxerror are removed.This threshold, however, represents the maximum possible error associated with the rockfall area.Jahne (2000) defined the variance σ 2 x of the position of a single point in an image (cell), introduced by the cell size dx, as assuming a constant probability density function within the cell area, i.e. all positions are equally probable.The standard deviation σ x is approximately 1 √ 12 ≈ 0.3 times the cell size.Therefore, to accommodate for uncertainty in the position of the area of change within each boundary cell as a function of cell size, 2σ can be used as a threshold as follows: The volumetric error is hence Equation ( 22) shows that the number of border cells relative to the total number of cells within the area of change is critical in determining the net volumetric error.A higher ratio of border cells to the total number of cells results in a greater proportional area (and hence volume) error.While this volumetric error assessment is applied to rasters of 3-D-derived change, its use also extends to extraction of discrete events from DEMs of difference (DoDs).

Processing of near-continuous monitoring data
To test the influence of the filtering phases on the subsequent change detection, comparisons were undertaken between two Earth Surf.Dynam., 6, 101-119, 2018 www.earth-surf-dynam.net/6/101/2018/No appreciable change occurred between these two scans.As the cylinder length increases (from 0.25 to 0.50 to 10 m), the number of surfaces that the cylinder intersects increases (direction equal to the normal vector).All points within a 0.25 m radius would be included as cylinder points (circles) and the distance between their mean positions (squares) calculated.From top to bottom, this distance is 0.0011 to −0.1460 to 0.0938 m.Longer cylinders intersect multiple surfaces and therefore measure the distance between projected centroids that do not accurately represent the surface to which the query point belongs.aligned point clouds following point removal based on various thresholds.The theoretical distance between the two point clouds is assumed to be zero given that no rockfalls were observed between their collection.Any offset therefore represents uncertainty, which is quantified for this purpose as the standard deviation of a 3-D change detection.This uncertainty is plotted against the applied edge and hole threshold in Fig. 9a.As the threshold is lowered, points with higher EH values are retained and the offset between the two point clouds (pre-alignment) increases.The distribution of EH values across the cloud is presented in Figs.9b and S1 in the Supplement.Using a 1 m search radius, an inflection at the 95th percentile of points occurs at 5 × 10 −4 .As a threshold, this value typically removes 5 % of points, which, as depicted by the dashed line in Fig. 9, generally account for uncertainties > 0.5 m.In addition to identifying edge features on the cliff face, this also helps to delineate areas of occlusion.The point density in Eq. ( 1), k, is used to filter spurious "floating points" in the dataset (for example, birds or dust).k values < 4, the minimum number to accurately define a centroid with an associated error, were removed.
Points removed using the waveform deviation filter (Fig. S2) occupied similar locations across the point cloud and, once removed, also resulted in decreased uncertainty in change estimates.In Fig. 10a, the mean absolute distance between points is used to the represent the positional uncertainty between point clouds.This uncertainty is ∼ 0.02-0.03m for points where δ ≤ 25.Conversely, points where δ >∼ 25 exhibit more significant scatter, often approaching 2 to 3 times the level of uncertainty in the entire cloud.Removing points with δ > 25 (Fig. 10b) retains 98 % of points, which accounted for a standard deviation of error between point clouds of 0.18 m prior to removal.The sequence with which the edge-hole and waveform deviation filters are applied appears to have little bearing on the outcome and the subsets of points removed by each have common members.When combined and applied to the dataset in this study, the filters reduced the standard deviation of change measurements between two stable point clouds from 0.078 to 0.055 m, thereby lowering the LoD that could be applied during rockfall or deformation identification by ∼ 30 %.Although the mean offset between unregistered point clouds at these sites was 0.51 m, this was improved using the ICP alignment described in Sect.3.2 (0.0053 m).Over the entire dataset, an average registration error of 0.005 m between point cloud pairs was obtained (n to n + x, where 1 > x > 8987).Despite the alignment of all point clouds to the initial reference scan, no increase in error was observed through time.However, considerable alteration to the surface topography, via the occurrence of a single large event or the continued spalling of material over time, would require new reference scans to be assigned over shorter timescales.
Following filtering and alignment, the DAN VCL method significantly lowered the LoD (0.03 m) compared to that achieved using DoDs (1.04 m).These were created using the same pairs of point clouds rasterised at 0.25 m with each pixel containing > 1 point.In this instance, the LoD also shows a 5-fold improvement relative to the M3C2 algorithm applied using the same normal and fixed cylinder radius (LoD = 0.165 m).The LoD was derived for every sequential scan to ensure that no increase in registration or epistemic errors developed through the monitoring period.This value lay consistently between 0.01 and 0.03 m.The maxi-mum LoD, 0.03 m, was therefore applied to each point cloud to prevent recording erroneous pixels in the resulting rockfall inventory.Combined with a cell size of 0.15 m, this provided a minimum detectable rockfall across the survey area of 6.75 × 10 −4 m 3 .More than 180 000 detachments were detected using the highest frequency of scans (∼ hourly) over the 10-month monitoring period.The spatial and temporal distributions of rockfalls observed are shown in Fig. 11.

Understanding the influence of survey interval in magnitude-frequency relationships
In order to assess the influence of more frequent monitoring on the resultant volume frequency distribution, two inventories were compared.These were analysed over the same monitoring period, using scans separated by different intervals (T int ) of 1 h (hours) and 30 days.An increase in the number of small rockfalls and the proportional contribution of small events to the overall rockfall volume distribution is evident at T int < 1 h (Fig. 12a), with the volume contribution from small (< 0.1 m 3 ) rockfalls increasing from 67 to 98 % of the total.The power-law scaling exponent, β, increases from 1.78 (30 days) to 2.27 (< 1 h).Notably, while a rollover occurs at T int = 30 days, this is not apparent at T int < 1 h.Given that both sets of scans were processed using the same LoD for change detection, the comparison at this site demonstrates that the observed rollover occurs due to superimposition and coalescence of events when longer survey intervals are used.For all rockfalls observed at T int < 1 h, volume error was modelled according to Eq. ( 22) (Fig. 12b).Larger rockfall volumes exhibit a smaller error in proportion to their volume.Importantly, however, as the vast majority of rockfall volumes are between 0.001 m 3 (a minimum of 2 pixels) and 0.01 m 3 (a minimum of 14 pixels), the uncertainty in volume ranges between 80 and 160 % of the estimate.A consequence is that the total volumetric uncertainty over 10 months of the T int < 1 h rockfall inventory is greater than that collected at T int = 30 days.High-frequency monitoring where change is dominated by a high frequency of low-magnitude events is not well suited to accurate measurement of total change through time.The error estimates demonstrate that the uncertainty in volume is greatest for the datasets in which T int is low.For the highest-frequency dataset, the total estimated volume is 110.87 ± 52.44 m 3 (±47 %), while the total estimated volume for the 30-day dataset is 72.37 ± 27.51 m 3 (±38 %; Fig. 13).

Processing techniques for near-continuous surface monitoring
Improvements in near-continuous point cloud acquisition currently outstrip the development of techniques for data treatment and analysis (Eitel et al., 2016;Kromer et al.,  2017).Near-continuous TLS has the potential to generate a considerable number of point clouds (> 10 3 -10 4 ), representing a 2 to 3 order of magnitude increase in data volume over previous terrestrial lidar monitoring campaigns with lower temporal resolution (e.g.Teza et al., 2007;Abellán et al., 2010;Rosser et al., 2013;Royán et al., 2015).Key attributes of the techniques developed to process such datasets therefore relate to computational efficiency, the ability to automate processing, and minimising the accumulation of error between each survey pair.These have necessitated tailored approaches to lidar processing, such as those described here, which often differ from previous applications.The relative gains of each processing step applied are evident in a gradual improvement to the applied LoD.A 30 % improvement with the application of radiometric and morphological filters occurs due to the removal of points considered to be both less accurate and less repeatable.This approach to the removal of unreliable points can also be used in alternative processing techniques, for example, as mem-bers of the fuzzy inference approach developed by Wheaton et al. (2010) to quantify spatially variable DoD uncertainty.The approach to change detection adopted in this study also yielded an improvement to the LoD, highlighting the importance of precise calculation of difference, here as a function of the cylinder length, for scenarios in which multiple surfaces may be intersected by the same normal vector in the resulting point cloud.Such surfaces can be characterised as increasingly three-dimensional on the scale of the applied cylinder length (e.g.Brodu and Lague, 2012) or as rough surfaces that are surveyed at oblique viewing angles (Hodge et al., 2009).This problem is exacerbated for near-continuous monitoring, given that scanning from a fixed position increases the width of occluded zones on surfaces inclined away from the scanner, yielding higher offsets between measured surfaces inclined towards the scanner.Critically, adaptive change detection techniques are necessary to account for the variability in point cloud quality across surfaces surveyed from a fixed position, where the installation location may yield unfavourable target geometries.These may take several forms, including spatially variable LoDs, varying cylinder lengths, or varying cylinder widths.
The cylinder radius determines the degree of spatial averaging over change measurements and, as such, should be informed by the style and scale of movements under investigation.In theory, the smaller the radius, the finer the spatial detail that can be established.However, this comes with a compromise in that the increase in accuracy by accounting for neighbouring points is reduced, the likelihood of intersecting points in the second point cloud is reduced, and the statistical significance of calculations is reduced by only drawing on a small number of points.Future development of a variable cylinder radius may draw upon the same point density estimated during normal estimation, varying in size until ∼ 20 points, the minimum suggested by Lague et al. (2013), are found.Here, cylinder radii between 0.15 and 2.00 m were applied to scan pairs in which rockfalls had, and had not, occurred from East Cliff.For scans in which no change occurred, the standard deviation of differences across the point cloud was greatest for the 0.15 m radius (0.018 m), approaching the point spacing, but decreased and stabilised at a radius of 0.25 m (0.013 m).For scan pairs in which rockfalls had occurred, the size and shape of rockfalls were contrasted with the results of varying cylinder radii and a Hausdorff distance estimation.The Hausdorff distance measure itself is influenced considerably by the scan line spacing and the local point density but, for this purpose, it provides an indication of rockfall geometry with the smallest degree of spatial averaging.As the cylinder radius increases, the difference in shape relative to the Hausdorff distance also increases.While a radius of 0.15 m best approximates the size and shape of the rockfall, this value is too close to the scan line spacing at all but the most proximal regions of the cliff face.A search radius of 0.25 m was therefore selected, providing similar rockfall geometries to the 0.15 m.This emphasises the potential to apply a variable cylinder radius across the point cloud based upon local point density, which could be helpful in future research.
The applied method removes scans undertaken during poor weather conditions (e.g.rain or fog) in order to preserve the accuracy of the resulting rockfall inventory.Given the potential for rockfalls to occur during poor weather conditions, this constitutes an important drawback of nearcontinuous TLS monitoring for rockfall inventory compilation.Techniques that can operate during inclement weather conditions, such as ground-based interferometric syntheticaperture radar (InSAR), are therefore better suited to maintaining temporal consistency in rockfall datasets.However, while the precision and measurement frequency of InSAR monitoring surpasses that of lidar, highly precise change measurements are also spatially averaged across large pixel sizes, resulting in a minimum detectable rockfall volume several orders of magnitude higher (∼ 1 m 3 ) compared to terrestrial lidar datasets.While small magnitude changes that occur across large areas can be accurately characterised at high frequency, detecting the frequency density of small erosion events is therefore compromised.Point cloud generation for slope monitoring has been supplemented in recent years by the development of new photogrammetric techniques, in particular structure from motion (SfM; Niethammer et al., 2011;Westoby et al., 2012;Lucieer et al., 2013;Turner et al., 2015;Carrivick et al., 2016).When imagery is acquired from unmanned aerial vehicles, SfM has the advantage of far lower operational costs than TLS, minimising areas of occlusion that occur from ground-level monitoring and providing highly dense point clouds due to the potentially small distances between the unmanned aerial vehicle and the slope.At present, however, the technique requires further development before it can be deployed for near-continuous monitoring.

Future developments in processing near-continuous TLS collection
Given that small events present the highest aerial, and hence volumetric, uncertainty, the recent development of 3-D volume estimation from point clouds (Carrea et al., 2012;Benjamin et al., 2016) may play an important role in nearcontinuous scanning.The uncertainty of these techniques is determined by the precision of the point cloud, thereby eliminating uncertainty in object aerial extents due to linear interpolation into a fixed grid.However, these techniques also contain uncertainties that arise in part from the meshing approach adopted (Soudarissanane et al., 2011;Hartzell et al., 2015;Telling et al., 2017).Due to the dependence of these techniques on a minimum of four points to create a closed hull, fully 3-D techniques are also limited in their ability Earth Surf.Dynam., 6, 101-119, 2018 www.earth-surf-dynam.net/6/101/2018/ to resolve small, single point detachments.The development of scanners with increasingly small angular step widths and increased rates of point acquisition, however, will decrease the minimum resolvable detachment.At present, the 3-D clustering required to isolate points belonging to geomorphic change, combined with subsequent meshing of these points, comes at a considerable computational cost.These techniques therefore remain to be applied for > 10 scans (e.g.Carrea et al., 2012;Benjamin et al., 2016;van Veen et al., 2017).This study demonstrates the need to adjust the frequency of data collection and processing in accordance with the study objectives.Here, monitoring has been undertaken to detect near instantaneous discrete changes to the slope (rockfalls) where both the spatial and temporal resolution of monitoring are important.Longer-term total change is more prone to error when change accrues from many small events, and big changes can occur in both the short and long term.There is a lack of research into this trade-off in spatial and temporal resolution but approaches that allow this to occur would be helpful in the future.The collection of a high-frequency time series of scan data presents the opportunity to reduce uncertainty by averaging point positions through both time and space, as points are independent in neither space nor in time.This averaging can take the form of averaging the 3-D position of each point, as utilised here and in M3C2 (Lague et al., 2013), or the averaging of differences between points (Abellán et al., 2009;Kromer et al., 2015b).

Implications for rockfall magnitude-frequency
While the precise magnitude-frequency exponent reported is specific to East Cliff, the scale-invariant behaviour of rockfalls is similar to that observed in other rockfall inventories, albeit across a narrower range of magnitudes.Along the same stretch of coastline, previous monthly monitoring has yielded exponents of β = 1.43-1.91(Rosser et al., 2007;Barlow et al., 2012), similar to that identified in this study using the T int = 30 days inventory (β = 1.78).Both the exponent and presence of a rollover show a dependence upon monitoring interval on temporal scales considerably lower than the intervals used by Barlow et al. (2012).For rockfall distributions created from TLS surveys, Young et al. (2011) noted that the ability to resolve small-scale changes should not introduce a rollover because the smallest reported rockfall is larger than the minimum detectable event identifiable in change mapping.From a statistical perspective, this statement holds true as long as the frequency density is not estimated using a moving kernel, which enforces an extrapolation of density that extends one kernel half width beyond the range of the observations both below the minimum and above the maximum, introducing inflections in the frequency density at the tails (Lim et al., 2010).Here, a rollover in the magnitude-frequency distribution is identified for the T int = 30 days inventory.However, this rollover was not present in the T int < 1 h dataset.Only rockfalls larger than 0.03 m (LoD) × 0.15 m × 0.15 m (the area of each cell, which exceeds the minimum point spacing) were analysed, which equates to a volume of 6.75 × 10 −4 m 3 .Critically, this implies that where events coincide, or coalesce, or where the mechanisms driving change are not spatially independent, event frequency is partially determined by survey interval.Given the development of spatially contiguous rockfall scars that has been observed in this setting (Fig. 11) and in other studies (Rosser et al., 2007(Rosser et al., , 2013;;Stock et al., 2012;Kromer et al., 2015a;Rohmer and Dewez, 2015;Royán et al., 2015), the creation of magnitude-frequency distributions from near-continuous monitoring has the potential to generate improved understanding of the underlying mechanisms of rockfall failure.
Monitoring at lower frequencies may provide more accurate estimates of rates of total change over longer periods.This is related to both the longer and hence time-averaged conditions captured but also to the fact that the same level of change measured infrequently has less volumetric error than when measured frequently, particularly when change is accrued by many small, discrete events.A decrease in T int that approaches near-continuous monitoring, 1 h, results in a shift in the exponent of the inverse power law of rockfall volumes from to 1.78 (30 days) to 2.27 (1 h).With a maximum plausible volume error, uncertainty in total rockfall volume ranged from 20 to 160 % of the measured volume.Although critical to measure the full rockfall volume distribution, highfrequency monitoring in this setting is not suitable for measuring net volume loss as a result of large numbers of small events.In summary, magnitude-frequency analysis of rockfall volumes (Fig. 12a) indicates that more frequent scanning detects a greater proportion of smaller rockfall events.Consequently, more frequent scanning also presents increased uncertainty in cumulative volume.Cumulatively, this error can be significant relative to the total flux over the monitoring period given that the size distribution of rockfall volumes adheres to a power law.

Conclusions
The magnitude-frequency distribution of geomorphic change is an important descriptor of the relative efficacy of event sizes and the nature of the hazard that they pose.Improvements in the ability to resolve the magnitude of events have surpassed the ability to constrain event frequency over short time intervals (< days).However, increasing the temporal resolution of monitoring of a changing surface increases the cumulative error over the same monitoring periods, particularly where change is dominated by numerous small events.In this study, we have discussed the practicalities and techniques to reduce this error for near-continuously acquired 3-D monitoring data, using one www.earth-surf-dynam.net/6/101/2018/Earth Surf.Dynam., 6, 101-119, 2018 of the highest-temporal-resolution 3-D datasets collected to date.
The findings of our workflow are distilled here.Both morphological and radiometric filters can be effective in removing unreliable points, such as edges, surfaces of high incidence angle, and vegetation, the effects of which become increasingly prominent when scanning from a fixed position.Applying these filters lowered the standard deviation of change detected between two stable point clouds from 0.078 to 0.055 m.Scans with any degree of occlusion, arising from atmospheric conditions (e.g.rain) should be entirely removed to ensure that no rockfall, which are more probable during these conditions, are likely to be missed in pairwise change detection.The alignment of large numbers of scans (10 2 -10 4 ) to the first scan of the dataset prevents drift without detriment to the overall alignment accuracy, though this may not apply in settings where the overall slope morphology changes through the monitoring period.Threedimensional change detection improves the LoD in relation to 2.5-D change detection techniques (DoDs).However, in its application to an actively failing rock slope, the DAN VCL approach described here yielded a 5-fold decrease in uncertainty relative to M3C2, the effectiveness of which increases with the frequency of holes and edges in the dataset.
By comparing rockfall inventories collected at T int < 1 h and T int = 30 days, it is apparent that more frequent monitoring captures a higher proportion of small rockfalls, represented by a higher magnitude-frequency scaling coefficient.Importantly, both the size and shape of rasterised events determines the ability to accurately quantify their volume, with smaller events and events with a higher proportion of boundary cells producing a higher degree of volume uncertainty.The method proposed to quantify this uncertainty represents an important consideration for volume estimation during future near-continuous monitoring campaigns.Critically, a higher proportion of small rockfalls in an inventory increases volumetric uncertainty, which accumulates because of the increased frequency of these events.Net long-term eroded volume is therefore most accurately quantified between the first and last surveys of a monitoring period.

Figure 1 .
Figure 1.Map of Whitby with the area scanned delineated with red lines.A Riegl VZ-1000 scanner is installed within East Pier lighthouse.The targets installed for the SiteMonitor4D range correction factor estimation are illustrated (T1-T6) in addition to the weather stations.Whitby Abbey lies 180 m from the cliff top.Map produced using shape files from Ordnance Survey© Crown Copyright and Database Right 2016.Ordnance Survey (Digimap Licence).

Figure 2 .
Figure 2. (a) Image of the cliff taken 1 h before high tide on 25 November 2015.Horizontally bedded strata are evident, with upper beds stained orange from downslope wash from glacial till of variable depth.The lower buttress comprises shales and some sandstone, while the near-vertical upper portion of the cliff comprises outcropping sandstone and sandstone interbedded with carbonaceous muds.(b) Slope model of the cliff showing the area covered by the TLS (light grey) draped over a 3-D model of the cliff, surveyed from multiple positions along the foreshore (dark grey).The total area measured is 8561 m 2 , or 89 % of the cliff face area (9592 m 2 ).The cliff is ∼ 210 m across and ∼ 60 m high.

Figure 3 .
Figure 3. Flow diagram representing the stages of the optimised near-continuous monitoring change detection method.All stages following ASCII to MAT conversion were written in MATLAB®, with ICP alignment and rockfall vectorisation using the built-in functions pcregrigid and bwboundaries.Point clouds are initially rotated to become approximately planar across the x − z plane, enabling the removal of points outside a tight bounding cuboid and rasterising of the point clouds of change.This rotation also enables an efficient solution to subsequent normal direction ambiguity, where the y component of each normal vector should always point out of the surface.

Figure 4 .
Figure 4. Conceptual illustration of the significance of removing partial scans.While parts of these scans provide accurate estimates of surface change, if a rockfall occurs in an area of no data, the failure will be missed using pairwise change.These scans should therefore be removed prior to change detection of the scan database.

Figure 5 .
Figure 5. Search radii used for normal estimation across the rock face.(a) The radius for each point on the cliff at which the point clouds are most planar, with a mean value of 1.1 m, used to estimate the normal vector prior to change detection.This point cloud was used as a reference model, such that the normal radius of points in subsequent scans was assigned based on the radius of closest point in this scan.(b) Surface planarity at a radius of 1 m, where higher values indicate a more 3-D neighbourhood.These occur at inflections in slope profile and in areas of high local relief, such as the sandstone beds near the cliff top.Gaps in point cloud are zones of occlusion illustrated in Fig. 2b.

Figure 6 .
Figure 6.Conceptual variation in the distance along the normal for a rockfall.(a) Change detection along a surface-normal direction estimated from the pre-failure surface (Cloud 1).(b) The normal direction estimated using a planar, post-failure surface more accurately represents the direction of change than the post-failure surface vector due to the complexity of the pre-failure surface.With both change detections originating in the same approximate positions on the cliff face, the difference in vector lengths illustrates the sensitivity of the 3-D change measurements to the normal estimation.

Figure 7 .
Figure 7. Approach to change detection used in this study based on empirical rockfall data.(a) A 2 m wide transect taken midway up East Cliff.The black points are taken from Cloud 1 and the grey from Cloud 2, with a 1.75 m high rockfall shown.Points within the cylinder radius, which intersects the two clouds, are shown as red and blue.The cylinder axis, which travels through the query point, is also shown.(b) Area of interest selected from (a) the centroids of each point cloud are determined and their orthogonal projection onto the normal vector (cylinder axis) is estimated (dashed lines).The distance measured in this study is between these projections, along the normal.

Figure 8 .
Figure 8. Inputs used for distance estimation with varying cylinder lengths.No appreciable change occurred between these two scans.As the cylinder length increases (from 0.25 to 0.50 to 10 m), the number of surfaces that the cylinder intersects increases (direction equal to the normal vector).All points within a 0.25 m radius would be included as cylinder points (circles) and the distance between their mean positions (squares) calculated.From top to bottom, this distance is 0.0011 to −0.1460 to 0.0938 m.Longer cylinders intersect multiple surfaces and therefore measure the distance between projected centroids that do not accurately represent the surface to which the query point belongs.

Figure 9 .
Figure 9. Sensitivity analysis for the edge-hole filter.(a) A change detection is undertaken between two point clouds where no observable movement occurred.The standard deviation for a single point cloud therefore indicates the level of noise between the two.When high EH values are retained, the standard deviation of change increases.(b) Cumulative proportion of EH values within an entire point cloud.Each line represents a different neighbourhood radius search, with a 1 m radius selected to ensure a minimum of four points, the minimum needed to estimate the CoG, would be found.An inflection in the number of points retained is used to define the threshold at 95 %, ensuring that artefacts such as holes are not introduced into the point cloud by over-removal of points.While the EH values change with the neighbourhood search radius, their distribution across the point cloud remains the same due to the normalisation by point density.

Figure 10 .
Figure 10.Sensitivity analysis for waveform deviation filter.(a) Mean absolute distance between two point clouds, attributed to points of each waveform deviation, from 1 to 50.Similar to Fig. 9, this indicates the comparison uncertainty between both scans.Error increases from ∼ 0.03 to 0.06 m at values > 25.The variability in error also increases such that the selection of an appropriate threshold > 25 is not possible.(b) Cumulative distribution of waveform deviation values.A threshold of 25 removes 2 % of points.

Figure 11 .
Figure 11.Distribution of rockfalls across East Cliff monitored at sub-hourly intervals between 5 March and 30 December 2015.Rockfalls are distributed across the entire cliff face, in particular in areas of exposed bedrock.Although the high water mark is below the portion of cliff shown in this figure, the largest and most frequent rockfalls occur at the base of the cliff.Accumulation and loss of material in the areas of unexposed bedrock on the cliff buttresses, which runs across the cliff face at ∼ 17 m elevation, were removed.Colours represent the age since 31 December, where red represents the oldest rockfall.

Figure 12 .
Figure 12.(a) Magnitude-frequency distribution for rockfall inventories acquired at varying T int showing that a higher proportion of small events is established by monitoring at high frequencies.Rockfalls used range from 5 March to 30 November 2015 to enable direct comparison between T int < 1 h and nine T int = 30 days change detections.(b) Rockfall volumes from < 1 h rockfall inventory.Percentage volume error is estimated using the LoD, number of internal pixels, and number of edge pixels.Frequency densities (kernel density estimates) are appended to each axis, showing that rockfall volumes < 0.01 m 3 account for the greatest proportion of measured rockfalls (modal volume = 0.0081 m 3 ).As a result, errors range from ∼ 60 to ∼ 140 % for most rockfalls (modal error percentage = 109 %).Cumulative volume estimations using rockfalls of this size may vary by at least the actual volume.

Figure 13 .
Figure13.Cumulative rockfall volumes measured though the monitoring period, using data from all 11 monitoring intervals.The results show that far higher volumes of material, up to twice those recorded by 30 days of monitoring, are measured at sub-daily intervals.The times of pairwise change detections are recorded as the date of the first scan, rather than the second.As a result, although all scan intervals record a significantly increased rate of rockfall activity during November, this appears earlier on the plot for longer scan intervals.The total estimated volumes are not included for comparison as change detections cannot be recorded up to the final day of monitoring for longer time intervals (30 December).