Mapping and validating lineaments

Data enhancement can help to emphasize edges that correspond to contrasts in acoustic impedance, magnetic susceptibility, or bulk density. Such contrasts often indicate the presence of important geologic boundaries. Emphasizing the edges can help with mapping and interpretation of those boundaries.

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 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 shows that lineaments derived using the maxima of the hyperbolic tilt angle are offset systematically from those derived using the total horizontal derivative. Which should I use?

Using synthetic gravity data, Cooper and Cowan demonstrate that no single-edge detector method is a perfect geologic-contact mapper. Citing Pilkington and Keating, 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.

Background

I will use gridded gravity data from a survey acquired in the Monti Romani of the southern Tuscany mining district. 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.

Figure 1. (a) Residual anomaly data; (b) dz, vertical derivative of the data; (c) dx, horizontal derivative in the east-west direction; (d) dy, horizontal derivative in the north-south direction. Grid node spacing is 100 m.

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.

Collocation mapping method

Figure 2. (a) TDXN; (b) threshold or binary map of TDXN maxima (black = 0, white = 1); (c) lineaments derived from (b); (d) collocation of results from TDXN, THDR, and NSTD, yellow when all three methods collocate, orange when any two methods collocate.

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 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:

>>> import numpy as np

>>> from fatiando.gravmag import transform

>>> data = np.loadtxt(‘continued.txt’)

>>> d = data.ravel()

>>> zderiv = transform.derivz(x, y, d, data.shape)

>>> zderiv2D = np.reshape(zderiv, data.shape)

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.

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:

>>> tdx = np.sqrt(dx**2 + dy**2)

>>> tdxn = np.real(np.arctan(tdx/np.absolute(dz)))

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. 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.

>>> from skimage.morphology import disk, closing

>>> binary_tdxn = color.rgb2gray(tdxn) > 0.75

>>> closed_tdxn = closing(binary_tdxn, disk(2))

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, and NSTD, the normalized standard deviation. I chose those three methods from the many available because their solutions are considered nonredundant.

>>> from skimage.morphology import dilation

>>> from skimage.morphology import skeletonize

>>> skel_tdxn = skeletonize(closed_tdxn)

>>> dil_skel_tdxn = dilation(skel_tdxn, disk(2))

>>> collocation = (dil_skel_tdxn + dil_skel_thdr + dil_skel_nstd)

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

>>> finalmap = np.ma.masked_where(collocation < 2, collocation)

>>> finalmap = finalmap.filled(0)

Discussion

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).

Figure 3. (a) Forward-gravity model across an approximately southwest-northeast profile, shown in red; (b) TDXN and collocation attributes extracted along the profile; (c) 2.5D geologic model. The collocation attribute corresponds well to the positions of the faults.

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.

Acknowledgements

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.

External links

• Corresponding author: Matteo Niccoli matteomycarta.ca, MyCarta.
• Madagascar workflow - reproducible using Madagascar open-source software

find literature about Mapping and validating lineaments