To invert the attitude of a geological trace

In 2016 I wrote a small utility in C++ and Fortran, geoSurfDEM, that calculates the intersections of a geological surfaces with a topographic surface (in the module "IntersectDEM", see post "geoSurfDEM: a C++ console application for determining intersections between 3D geological surfaces and topography") and that inverts the attitudes of geological traces (module "BestFitGeoplanes", see post "A Linux tool for calculating local best-fit plane attitudes from geological traces"" for the original description).

Provided a set of points of a geological trace, defined by their coordinates x-y-z, the BestFitGeoplanes module derives the local best-fit-planes, using a mathematical technique known as Singular Value Decomposition (SVD).
The used Fortran Lapack subroutine derives the eigenvectors and eigenvalues of the set of the source points expressed as a matrix (m x 3, where m is the number of points and 3 is the dimension of the embedding space, i.e., a 3D Cartesian space).

The eigenvector normal to the plane with the minimum data variance, i.e., the plane that better approximates the data assuming that they approximately subplanar, is the one in which we are interested.

The question is, how reliable is the result? If my data are collinear, no unique solution exists. The same applies when all the points are approximately coincident.

In the original version, in order to characterize the situations related to spurious results, the module calculated the spatial range of the source data along the three reference axes and also the volume of the bounding box (axes-aligned).

To verify whether they can really provided indications of errors (i.e., large differences between expected and calculated values),  a synthetic data set, made up by two geological traces, was inverted and the results analyzed.
The two synthetic geological traces derive from the analysis of the Timpa di San Lorenzo fault structure (Calabria, Southern Italy), where the northern and the southern fault structure segments attitudes can be approximated by two planes, respectively oriented 067.3°/39° and 077.3°/40° (with respect to the UTM 32632 CRS up direction)(see previous post "Along-trace profiles of the Timpa San Lorenzo fault structure (Calabria, Italy)"). The used DEM is a subset of the TINItaly DEM.

 

The intersection of two theoretical attitudes with the DEM were calculated using the 'DEM-plane intersections' of the qgSurf QGIS plugin, with the results saved as two point shapefiles. Points outside the fault spatial range were removed, the two layers were then merged into a shapefile and the data were exported as an x-y-z csv file, to be used as input for the 'BestFitGeoplanes' module.

The inversion results reproduce the theoretical attitudes very well (< 1° of difference), except for 4 values out of the total 311 (1.3 %). The spatial ranges of the source point data for each result do not appear however to be predictors of the result reliability.

To find a better quality indicator, a new code version outputs the values of the three eigenvalues for each result.
Generally the first eigenvalue (S1) has a value in the 10-50 range, the second (S2) is 1 or less, and the third eigenvalue (S3) is very low or zero.

In the analyzed data set, the "erroneous" results are characterized by low S2 values, around e-5. By considering the logs of the ratio between S2 and S1, the results with values lower than -5 are from grossly (>> 1° misfit) to mildly (1° or less) disoriented (Fig. 3). 

Indices comprised between -3.2 and -5.0 are correlated with mildly disoriented results (< 1° disorientation).

So using a threshold of -3.0 for this data set filters out all the low-quality.

Obviously, to check whether this empirical threshold can be considered general requires to consider other test cases. Possibly a collinearity index, to be defined, could be used an additional or alternative quality indicator.

References

 

geoSurfDEM - https://gitlab.com/mauroalberti/geoSurfDEM