Mapping and validating lineaments

Not all geologic boundaries are necessarily contrast boundaries that produce data edges, however. In addition, the location of interpretable edges is influenced by many factors, one of which is the enhancement method used. As an example of this, Figure 2 of Niccoli (2012) shows a qualitative correlation between occurrences of antimony mineralization in the southern Tuscany mining district and the distance from lineaments derived from the total horizontal derivative. However, Figure 1 of Niccoli (2012) shows that lineaments derived using the maxima of the hyperbolic tilt angle (Cooper and Cowan, 2006) are offset systematically from those derived using the total horizontal derivative. Which should I use?

Using synthetic gravity data, Cooper and Cowan (2006) demonstrate that no single-edge detector method is a perfect geologic-contact mapper. Citing Pilkington and Keating (2004), Cooper and Cowan (2006, p. 1586) conclude that the best approach is to use “collocated solutions from different methods providing increased confidence in the reliability of a given contact location.”

In this tutorial and the accompanying IPython Notebook (at github.com/seg), I show how to (1) create a full suite of derivative-based enhanced gravity maps, (2) extract and refine lineaments to map edges, and (3) combine them into a single collocation map to increase confidence in the edge locations.

I will use gridded gravity data from a survey acquired in the Monti Romani of the southern Tuscany mining district (Niccoli, 2000). The raw measurements from 93 stations were reduced with a standard workflow to Bouguer gravity, from which a regional trend was subtracted to remove the effects of deeper sources. Figure 1 shows the resultant residual anomaly map.

In this area, the complex compressional deformation of the Ap-ennines caused doubling of the Paleozoic metamorphic basement, placing it in lateral contact with the less dense rocks of younger units (the Triassic-Oligocene Tuscany Nappe and the Cretaceous-Eocene Liguride Complex). Figure 1a reflects this. Roughly speaking, there is a high in the southeastern quadrant of about 3.0 mGal corresponding to the location of a large basement outcrop, a northwest-southeast elongated high of about 0.5 mGal in the center bound by lows on both the southwest and northeast, and finally a local high in the northwestern quadrant of about 0.5 mGal.

Using mutually orthogonal forward gravity models, Niccoli (2000) postulates that a system of postorogenic normal faults caused differential sinking of the top of basement in different blocks, leaving an isolated high in the middle, which is consistent with the tectonic history known from the literature. One of the objectives of the study was to help refine the mapping of those faults.

First we load the residual anomaly grid (Figure 1a) and calculate the vertical and horizontal derivatives (Figures 1b, 1c, and 1d). From these, we can generate several enhanced data sets, including the amplitude-stabilized total horizontal derivative, TDXN (Cooper and Cowan, 2006) shown in Figure 2. The full code to generate TDXN and the other enhanced data sets is included in the IPython Notebook.

As an example, the key operations required to load the data and calculate their vertical derivative, dz, are listed in the code block on the next page:

All the above operations use functions from NumPy, a Python library of mathematical functions, and Fatiando a Terra, a Python toolkit for geophysical modeling and inversion by Uieda et al. (2014).

The operations in the next block of code calculate the total horizontal derivative, TDX, and from this, the amplitude-stabilized total horizontal derivative, TDXN, shown in Figure 2a:

Next we threshold TDXN using a scalar value to derive a binary map of TDXN maxima, which we further enhance in several ways using morphological image-filtering operations from the Python image-processing module, scikit-image (van der Walt et al., 2014). For example, closing (a dilation followed by an erosion) removes the smaller dark spots using a disk-shaped structuring element. More details on these morphological operations are included in the IPython Notebook. Figure 2b shows the cleaned binary map.

In the next block of code, we extract lineaments, shown in Figure 2c, using skeletonization, dilate the skeleton so that the lineaments are three pixels wide, and finally sum up these TDXN-derived lineaments with those from two other enhancement methods: THDR, the total horizontal derivative of tilt angle (Verduzco et al., 2004), and NSTD, the normalized standard deviation (Cooper and Cowan, 2008). I chose those three methods from the many available because their solutions are considered nonredundant (Pilkington and Keating, 2010).

Figure 2d shows the resultant collocation (confidence) map with a mask applied so as to display contacts only where at least two methods collocate:

Figure 3 shows a 2.5D forward model along a profile from Niccoli (2000), to which I added in panel (b) a collocation plot extracted along the same profile, yielding increased confidence in the location of geologic contacts in the model shown in panel (c).

In this case, forward modeling had been carried out first and the work on edges and collocation afterward, to be used as validation. The ideal strategy, however, would be to do the work on the edges and collocation first and then extract both along the desired profile and use them to place vertical boundaries between bodies of different density. These boundaries then could be modified based on geologic knowledge of the area and other controls (magnetic or induced-polarization profiles, exploratory drilling, and so forth).

This technique also could be used to collocate lineaments derived from either reduced-to-pole magnetic data or seismic data (e.g., Russell and Ribordy, 2014).

Thank you to Evan Bianco, Gordon Cooper, Tooney Fink, Matt Hall, and Leonardo Uieda for their comments and suggestions and to Michele Di Filippo, my first mentor, for getting me “psyched” on gravity exploration and structural geology.