Robust best-fit planes from geospatial data

Geological and geomorphological surfaces observed in nature are rarely perfectly planar. Nevertheless, in some situations, it can be useful to approximate surfaces as planes, for example, to improve computational efficiency (e.g., in discrete fracture network modeling), or when analyzing the spatial distribution of deviation from planarity (e.g., Jones et al., 2009).

To define the orientation of a plane requires three noncollinear points; this is the basis of the “three-point problem” used widely in geological mapping (e.g., Threet, 1956; Berger et al., 1992) and taught in undergraduate classes and textbooks (e.g., Lisle, 1988; Boulter, 1989). Most geospatial methods typically allow many more than three points (10¹–10⁷ points or more) to be rapidly measured on an exposed surface. Relevant technologies include reflectorless total stations (e.g., Feng et al., 2001), laser range finders (e.g., Xu et al., 2001), global navigation satellite systems (including global positioning system [GPS]; e.g., McCaffrey et al., 2005), and satellite imagery (e.g., Banerjee and Mitra, 2004), as well as virtual outcrop methods such as terrestrial and airborne light detection and ranging (lidar; e.g., Ahlgren et al., 2002; Hasbargen, 2012; Cracknell et al., 2013), digital photogrammetry (e.g., Pringle et al., 2004; Bemis et al., 2014), and combined methods (Bistacchi et al., 2011; Jacquemyn et al., 2015).

The purpose of this paper is to review different regression methods used to derive a best-fit plane from geospatial data, and to assess the criteria used to determine the quality of the regression. We propose new criteria to replace those described by Fernández (2005).

One common method of plane fitting is based on ordinary (or “standard”) least squares regression, initially developed by Legendre (1806) and Gauss (1809). The approach is easiest to illustrate in terms of two variables (Fig. 1A), where x is the independent variable, and y is dependent (or “measured”), though it is equally applicable in three (or more) dimensions. This method is the basis for the planar regression method used in section 2 of Fernández (2005).

An important limitation of the ordinary least squares method is that it assumes that all error lies in the dependent dimension, and therefore the regression will minimize the residual errors in that variable (Fig. 1A). For most types of geospatial data, this assumption is invalid, since there are measurement errors in the x, y, and z dimensions. Consequently, in many situations, a more appropriate strategy is to use “total least squares” (or “orthogonal,” or “errors-in-variables”) regression, in which the residual errors are minimized in a direction perpendicular to the best-fit plane (Fig. 1B). The total least squares approach is not new (Adcock, 1877, 1878; Pearson, 1901), and according to Anderson (1984), it has been rediscovered many times (e.g., Sturgul and Aiken, 1970; see also subsequent discussion in Gower [1973, and references therein]). Total least squares is the underlying basis of principal component analysis (Hotelling, 1933; Jolliffe, 2002), and the moment of inertia approach described in section 3 of Fernández (2005).

For many types of three-dimensional (3-D) geospatial data, using total least squares regression, in which the likelihood of error is equal in the x, y, and z dimensions, will be a valid approach to plane fitting. For example, in virtual outcrops derived from terrestrial lidar and/or digital photogrammetry data, we usually consider there to be no systematically greater error in any one dimension relative to the others. For types of data where errors are unequal (such as GPS data, where measurement precision is poorer in z than in x and y), a weighted total least squares approach can be used (e.g., Markovsky and Van Huffel, 2007). Other modifications to the total least squares approach include ways to identify and disregard outlying points (e.g., Li and De Moor, 2002).

A significant advantage of using orthogonal regression is that the magnitudes of the residual errors are independent of the orientation of the best-fit plane with respect to any of the reference coordinate axes. In contrast, with ordinary least squares, the magnitude of the residuals varies with the orientation of the plane relative to the chosen coordinate axes (see section 4.1 of Fernández, 2005).

(where xyzPoints is the input matrix of geospatial data points, and V and D are the eigenvectors and eigenvalues of the covariance matrix of the input data). The best-fit plane passes through the centroid of the input data. The maximum and intermediate eigenvectors lie within the best-fit plane, and the minimum eigenvector defines the pole to the plane.

As discussed by Fernández (2005), there are two aspects in determining the quality of a best-fit solution: the suitability of the sampled data points to define a plane, and the goodness of fit between the input data and the resultant regression plane. The suitability of the input points depends upon their spatial distribution in three dimensions, while the goodness of fit is related to the magnitude of the residual errors. The eigenvalues of the covariance matrix (λ₁ > λ₂ > λ₃) derived from total least squares regression are useful in analyzing both these aspects of the best-fit results.

One way to conceptualize the robustness of a best-fit plane is to consider the dimensionality of the sample data (cf. Jones et al. 2008). Since points within a plane are intrinsically two-dimensional (2-D) in extent, best-fit solutions for data points that sample a space that tends toward one dimension (1-D) or 3-D will have poorer quality (Fig. 2). If the data points are nearly collinear (Fig. 2A), such that λ₂ and λ₃ are small, then the sampled points are unsuitable for defining a best-fit plane, irrespective of the magnitude of the residual errors (in this situation, the covariance matrix is close to a rank-one matrix, and the plane is degenerate). In contrast, if the data are noncollinear (λ₃ << λ₂), and the total residual errors are low (λ₃ is small; the covariance matrix is close to a rank-two matrix), then the input data are nearly coincident with a 2-D plane (Fig. 2B). Finally, where the residual errors are larger (λ₃ >> 0), then the data points are markedly non-coplanar (i.e., the dimensionality of the data is >2), and the goodness of fit is poor (Fig. 2C).

As with all regression methods, for orthogonal regression there is no universally applicable measure of what constitutes an acceptable fit, so the quality of the best-fit solution is arbitrary and depends on the nature of the input data and the purpose to which the interpreted best-fit planes will be used. We illustrate this with analysis of fractured outcrops based on high-resolution outcrop data acquired using terrestrial laser-scanning (lidar) or multiview photogrammetry.

Here, we describe an approach that we typically use to assess best-fit planes derived from fracture surfaces measured from lidar and/or photogrammetry point-cloud data, as part of a workflow to produce deterministic discrete fracture network (DFN) models from outcrop analogues (Fig. 3). For illustrations of the lidar equipment used and further details of the method of field acquisition, see Labourdette and Jones (2007), White and Jones (2008), and Jones et al. (2009). Multiview photogrammetry was described by James and Robson (2012).

Following field acquisition, lidar data and/or digital photographs from multiple tripod positions are coregistered, georeferenced using differential GPS control points, and processed to give a point cloud that represents a digital copy of the outcrop surface (a “virtual outcrop model”; Fig. 3B). Representative points from each exposed fracture surface or intersection trace are picked (for examples, see Fig. 3C) and form the input data for the calculation of the best-fit plane for that fracture (Fig. 3D). For further examples of fractures picked from virtual outcrop models and analyzed using the methods described here, see Jones et al. (2009), Jacquemyn et al. (2012), and Pless et al. (2015). Other examples using alternate workflows include Olariu et al. (2008), Wilson et al. (2011), and Seers and Hodgetts (2014).

A critical step in the analysis is to decide two appropriate thresholds for the quality of the plane fitting. First, to ensure that the sampled data points are sufficiently noncollinear (cf. Fig. 2A), the key requirement is that the extent of the data in a direction parallel to the intermediate eigenvalue (λ₂) is large compared with the precision of the input data. Therefore, it is important to be able to quantify the typical measurement precision of the method used; this depends on the specific type of lidar equipment, or the errors involved in photogrammetric processing (James and Robson, 2012). For example, high-precision medium-range scanners such as the Leica ScanStation C10 or Leica HDS 3000 (see fig. 1A inWhite and Jones, 2008) have low noise levels and, according to the manufacturer, have an accuracy of ∼2–6 mm at 50 m range (see also Lindenbergh et al., 2005). With these types of scanners, our empirical tests have shown that a low λ₂ threshold can be used; e.g., input data with λ₂ > 0.001 typically give well-defined planes throughout the scanners’ measurement range of ∼200 m. In comparison, when using a scanner prone to slightly more measurement noise, such as the Riegl LMS-Z420i (fig. 3 inJones et al., 2009), our tests indicate that a λ₂ threshold of 0.005 or 0.01 is often more appropriate. According to Riegl, under test conditions, this scanner can achieve an accuracy of ∼10 mm at 50 m (see also Renard et al., 2006). In situations where merged lidar data are smoothed (decimated), higher λ₂ threshold values may need to be used.

However, for the statistical analysis of hundreds or thousands of fractures in relation to DFN modeling, we generally find it more suitable to allow larger fracture planes to have higher variance, and therefore for this application, we usually choose to define the threshold in terms of the ratio of λ₂ to λ₃ (e.g., data points with λ₂/λ₃ < 10 are rejected).

In our experience, the criteria used herein to assess the quality of best-fit planes consistently give robust results. This is based on the analysis of many virtual outcrop data sets derived from lidar and/or digital photogrammetry, in a range of lithologies including carbonates, sandstones, mudstones, and crystalline rocks. Other criteria might also be useful and may be equally valid for different purposes.

Figure 4 also illustrates how the K parameter is particularly inappropriate, because planes with high λ₁ are likely to be rejected as unreliable because the resultant K values will be high, although these planes are actually often very well defined in virtual outcrop data. It is common for outcrops to contain fractures (or other types of planes) that can be reliably delineated by picking along the trace of the intersection of the plane with the irregular outcrop surface (e.g., fractures “b” and “d” in Fig. 3C). These intersection traces are often very long relative to their width, but as long as the outcrop surface has sufficient rugosity along the trace, then the definition of the plane is robust. For example, the intersection of the fracture in Figure 4A with the outcrop surface is very elongate. This results in a high λ₁ value, and therefore a K value that is high enough to be deemed unreliable. In reality, however, the only critical factors are that λ₃ should be small, and that there is sufficient relief on the outcrop surface such that λ₂ is large in absolute terms in comparison with the spatial precision of the measurements; the relative size of λ₂ in comparison with λ₁ is irrelevant. This is further demonstrated in Figures 4B and 4C. The fracture shown in Figure 4B is defined by an equant distribution of points so that its K value is low (provided that λ₃ is small). Assume that we are able to extend the analysis by picking an additional point situated on the same planar fracture (beyond the area of fracture initially used), depicted by the red square in Figure 4C. (Assume also that λ₃ remains small.) Counterintuitively, the addition of an extra point, which should help to improve the definition of the fracture, actually increases its K value significantly, such that the best-fit plane might become assessed as unreliable. Adding further additional points (magenta stars) continues to cause the K value to change, but it always remains higher (i.e., “less reliable”) than in Figure 4B.

Consequently, our opinion is that the M and K criteria of Fernández (2005) are generally inappropriate, and to assess the robustness of best-fit fracture planes in relation to DFN modeling, we use the following parameters, xL and either xP₁ or xP₂:

Parameter xL (precision-dependent test for collinearity): Reject the best-fit plane if λ₂ < ρ², (where ρ ∝ spatial precision of the data).

Parameter xP₁ (rigorous test for planarity): Reject the best-fit plane if λ₃ > Tv, (where Tv is the user-defined threshold for variance; i.e., reject if the mean of the squared residuals is deemed too high).

Parameter xP₂ (alternate, less rigorous scale-dependent test for planarity): Reject the best-fit plane if λ₂/λ₃ < C (where lower values of C allow more curvature or irregularity of the fracture surface).

These parameters are applicable to many plane-fitting problems involving 3-D data, including analysis of digital outcrop models. The specific threshold values to be used for both parameters depend on the scale of analysis, the precision of the geospatial data, and the desired stringency when assessing goodness of fit.

Total least squares regression, in which errors are minimized in the direction orthogonal to the regression plane, provides a robust way to define best-fit planes from many types of geospatial data. Eigenvectors and eigenvalues from the regression analysis are used to define the orientation of the best-fit plane and to assess the quality of the results.

The M and K criteria proposed by Fernández (2005), based on adaptation of a method to quantify clustering in orientation data, are inappropriate and are unable to identify reliable best-fit planes in a consistent way. Not only are the criteria scale independent (and hence do not take into account the spatial precision of the input data), but more fundamentally, the criteria do not represent independent measures of planarity: They are a measure of the three-dimensional aspect ratio of the sampled data points. The criteria are only appropriate in the particular theoretical case when there is no measurement error in the data, and when all the planes in a study are represented by sample points with the same length to height ratio (such that the influence of the maximum eigenvalue, λ₁, is removed). Fundamentally, the method of Woodcock (1977) is a widely used, robust, and elegant way to characterize shape fabrics from a set of orientation measurements (i.e., vectors with no defined spatial position), but it is poorly suited to analyze the spatial relationship of georeferenced point data.

Instead, we propose a pair of criteria that in combination allow robust best-fit planes to be derived. First, the intermediate eigenvalue from the orthogonal regression must be sufficiently large in comparison with the spatial precision of the input data to ensure that the best-fit plane has sufficient breadth to be statistically meaningful (i.e., to test that the input data are noncollinear). Second, the minimum eigenvalue needs to be sufficiently small to ensure that the input points are satisfactorily close to being coplanar. This can be assessed using the standard deviation of the residual errors from the regression, although for some purposes, a less stringent approach based on the ratio of the intermediate and minimum eigenvalues gives appropriate results.

We are very grateful to Sabine Van Huffel for her input on total least squares regression. Eni S.p.A. funded field work, including lidar acquisition and data analysis, during Jacquemyn’s Ph.D. program, supervised by Rudy Swennen. Thomas Seers provided a careful, insightful, and helpful review, and we also thank an anonymous reviewer.