«LANDSLIDE DEFORMATION CHARACTER INFERRED FROM TERRESTRIAL LASER SCANNER DATA A DISSERTATION SUBMITTED TO THE GRADUATE DIVISION OF THE UNIVERSITY OF ...»
McCoy, S.W., Kean, J.W., Coe, J.A., Staley, D.M., Wasklewicz, T.A., and Tucker, G.E., 2010, Evolution of a
natural debris flow: In situ measurements of flow dynamics, video imagery, and terrestrial laser scanning:
Geology, v. 38, p. 735-738.
DETERMINING GROUND DISPLACEMENT FIELDS OF SMALL SPATIAL EXTENT
USING TERRESTRIAL LASER SCANNER DATA: A COMPARISON OF 3D METHODS
APPLIED TO LANDSLIDE MONITORING
3.1 Introduction As discussed in the first two chapters of this dissertation, terrestrial laser scanning (TLS) is an emerging technique for detecting surface displacements of small spatial extent (meters to subkilometers) more accurately (sub-cms to cms) at higher spatial and temporal resolutions (e.g., Oldow and Singleton, 2008; Stewart et al., 2009; McCoy et al., 2010; Aryal et al., 2012).
Geologic examples of these types of small spatial extent surface displacements include land subsidence, active faults and volcanoes, glacier movement, and landslides (Figure 3.1).
However, significant difficulties may arise when deriving 3D displacement fields using TLS data primarily because the data are not necessarily from the exact same location of the reflector between the observational epochs due to changes in the scanner’s orientation and/or changes in the reflective surface (Figure 3.2). Furthermore, returns from vegetation that change over space and time complicate the use of TLS data for displacement analysis. Therefore, matching of a pattern or surface is needed to analyze the TLS data.
There are various approaches to derive surface displacement fields using TLS data but each is limited, and there is no accepted best practice for automated analysis. DEM differencing is one of the most common approaches in the literature, but it is a scalar measurement of displacement along a single, vertical axis (Baldo et al., 2009; Prokop and Panholzer, 2009; McCoy et al., 2010). Therefore, this technique is appropriate for the rare case when expected motion is only in one direction (vertical in most cases but it can be any direction). In general, landslide surface displacement fields cannot be determined using vertical DEM differencing. Another approach uses manual feature tracking to estimate displacements of identifiable features including treetrunks or user-installed reflectors such as large spheres in each scan (Collins et al., 2009;
Wilkinson et al., 2010). This approach can be quite precise, particularly if there are adequate identifiable features and measurement scatter from the reflective object is damped using geometric modeling of the feature. This technique is not automated and therefore it is time consuming with results that are user dependent. Least squares 3D surface matching (Gruen and Akca, 2005) has been applied to TLS data (Monserrat and Crosetto, 2008), but this technique seems to work well only for tracking features with regular shapes. Two of the most promising approaches to estimate 3D displacement using point cloud data are particle image velocimetry (PIV) (Aryal et al., 2012; Aryal et al., 2013) and iterative closest point (ICP) (Teza et al., 2007;
Nissen et al., 2012), but the strengths and limitations of applying these methods to TLS data analysis have yet to be explored.
In this chapter, we test the efficacy of using PIV and ICP methods to estimate 3D surface displacement from TLS data. In particular, we compare the PIV and ICP methods applying synthetic signals to TLS data from the slow moving Cleveland Corral landslide (CCL) in California to compare the performance of PIV and ICP. Then we apply both methods to TLS data from the CCL and compare the results with independently measured GPS and feature tracking displacements. We also use the TLS-derived displacement fields to compute strain fields and characterize the surface deformation pattern of the toe part of the landslide in space and time.
Movement rate and spatial extent of the most common geologic features. Spatial extent of most landslides is too small to use space-based InSAR to detect variations in surface displacement.
TLS point cloud data of a stationary building from two temporally different acquisitions (red and blue dots). The building is identifiable in both scans, but there is neither one to one relationship between data points nor is there an equal number of data points in the two scans.
3.2 Displacement Estimation Methods
Two of the most promising approaches in the literature to estimate 3D displacement using point cloud data are PIV (Aryal et al., 2012; Aryal et al., 2013) and ICP (Teza et al., 2007; Nissen et al., 2012). Here, we compare these two methods by applying a synthetic signal to a TLS scan and also by using the series of TLS scans from the slow moving CCL in California.
3.2.1 Particle Image Velocimetry
The PIV method has been used for decades to derive the velocity of fluid flows seeded with particles from time series photography (Keane and Adrian, 1992; Westerweel, 1997; Meunier and Leweke, 2003; Raffel et al., 2007). Fundamentally, PIV estimates a velocity field in a plane by cross-correlating a subset of raster images from a series of observational epochs. The PIV method has also been applied to geologic studies with close-range photography (White et al., 2003). Recently, Aryal et al. (2012) adapted the PIV method for 2D TLS data displacement estimates and Aryal et al. (2013) extended the method to 3D. To apply PIV to TLS data and
estimate 3D displacement field, we perform the following steps:
1. Grid the aligned or referenced 3D (x,y,z) point-cloud data with grid size, GR, in the horizontal plane to acquire images I1(i,j) and I2(i,j) where each grid-value contains average z-values from the corresponding TLS data set. Generally, smaller GR is better, as the correlation can potentially introduce estimation error of +/- 0.5 GR, but coarse GR allows faster cross-correlation and gridding. Gridding should be done without data extrapolation as it can introduce significant errors. In this study, we use a GR of 0.04 m.
2. Cross-correlate a window size of WC from the image I1 with an interrogation window of size WI from the image I2 for each grid shift (is, js) to acquire the normalized crosscorrelation function ( rN ) given by
where µ and are mean and standard deviation of z values of respective images indicated by the subscripts I1 or I2. The horizontal components of displacements are then a distance to the peak in the cross-correlation matrix from its origin (no shift position). To acquire the displacement at sub-pixel accuracy, a Gaussian function is fitted to the crosscorrelation matrix and the peak of the Gaussian function is located.
Any non-uniform displacement or displacement gradient within a correlation window can influence the classical PIV results causing the peak in the correlation matrix to be broad or even have multiple peaks. This can cause the estimated displacements to be inaccurate.
The iterative deformations of the correlation window (Huang et al., 1993; Meunier and Leweke, 2003) overcome this problem by applying the cross-correlation in larger windows at the first step followed by the correlation in smaller windows in the second step. Selecting the correct size of WC and I1 can be specific to a data set, the expected displacement, and the displacement gradient. As a rule of thumb: WC 2*dmax and I1 3*dmax in the first run where dmax is the expected maximum displacement in the window.
In the second run, the estimation is less sensitive to the window sizes and both WC from I1 needs to be smaller than in the first run. We refer reader to Aryal et al. 2012 for detailed parameter selection criteria.
3. Obtain the vertical component of the displacement by translating the elevation map according to the 2D horizontal displacements from step 2 and then differencing the elevation maps (Figure 3.3). This provides an approximate 3D displacement field.
To perform the second step above, we adapt the freely available DPIVsoft tool (Meunier and Leweke, 2003; Meunier et al., 2004) that has also been applied for TLS data (Aryal et al., 2012).
Conceptual sketch showing components of landslide displacement at surface of the sliding block. Ground surface displacement of G to G' consists of horizontal components ux and uy and a vertical component uz. The vertical component uz is the difference in elevation from G to G'. Elevation is available almost everywhere from TLS DEMs. Displaced location G' of vertical grid G is located using PIV-TLS estimated ux and uy.
3.2.2 Iterative Closest Point The ICP algorithm is one of the more commonly used methods for matching 3D point cloud data (Besl and McKay, 1992; Chen and Medioni, 1992; Zhang, 1994). Several commercial and research tools use ICP (e.g. Polyworks software by Innovmetric Inc.) to align TLS data.
Although there are different variants of ICP, its main goal is to refine the matching between two point cloud datasets (often referred to in the literature as model and data or destination and source) by estimating the best transformation (rotation and translation) parameters based on iteratively minimizing the distance between data points from two scans (Figure 3.4).
Let Mi=(m1, m2, … mn) and Di=(d1, d2, … dn) be two TLS scans where mi and di are composed of x,y, and z locations on the ground surface. In order to find the displacement, the goal is to find a rigid body transformation composed of a rotation matrix R and translation vector T so that M (model) and D (data) have the best alignment. The best alignment results when the sum of squared distance from points in one cloud to their nearest neighbors in the other point-cloud is
minimized, often referred as the error metric E such that:
E i 1 ( Di T M i ) 2 n (2) The error metric in equation (2) is the sum of the squared distances between corresponding points in M and D, often referred to as point-to-point minimization (Besl and McKay, 1992).
Finding corresponding points, however, is not trivial and can be computationally challenging.
Therefore, very often point-to-plane minimization (Chen and Medioni, 1992) is performed which sums the perpendicular distances of the data points to tangent planes containing the matched model points (Figure 3.4) and minimizes the error metric iteratively. Mathematically, the error
metric E is:
E i 1 ( Di T M i ) 2. ni n (3) where ni is the normal plane at the ith point in the reference point cloud. In equation (3), R is a function of nonlinear trigonometric functions but when the direction cosines in the x, y and z directions ( , , and are small, equation (3) can be solved using a linear least-square approximation (Low, 2004).
Sketch showing the point-to-plane distance matching in ICP between model and data that minimizes the sum of the Euclidian distance (dotted black lines) iteratively. Points represented by mi and di are from two data sets.
In this study, we perform ICP point-to-plane distance minimization using the commercially available Polyworks 10.1 software. First we align all the scans masking the potentially moving area and then import the initial scan as a model (reference data) image that results in a surface to which we fit data points from the second scan. We then divide the second scan into different subsets using square grids (e.g., 5x5 m). We use a sequence of operations in the Polyworks program, as described in Teza et al. (2007) to estimate the transformation matrix for each subset grid. Teza et al. (2007) contains a detailed description of the estimation scheme.
3.2.3 Synthetic Tests
To compare the performance of ICP versus PIV and to understand the effect of TLS data acquisition parameters, such as data density, we performed a test applying a synthetic displacement signal to TLS point cloud data from the toe portion of the CCL. ICP was originally developed for point cloud data and therefore its performance for purely synthetic data is well documented (e.g., Besl and McKay, 1992; Chetverikov et al., 2005; Minguez et al., 2006).
Similarly, use of PIV for purely synthetic data and its performance has been discussed in Aryal et al. (2012). Therefore for comparative purposes we applied a synthetic signal to actual TLS data instead of a purely synthetic example. We introduced a synthetic displacement pattern into the January 2010 point cloud (Aryal et al., 2012) with a maximum value of -1 m in y and -0.4 m in z-directions and recovered the signal using both PIV and ICP methods.