Subsurface Boundary Geometry Modeling: Applying Computational Physics, Computer Vision and Signal Processing Techniques to Geoscience

This paper is distinguished from prior work through its attempt in fully automating a chain of processes required to reconstruct subterranean surfaces using sparse labeled data. These processes include boundary extraction, spatial correspondence and contour evolution, as outlined in Fig. 1. Although the contours for various geological domains (henceforth referred as geozones) can be specified interactively, our model basically requires only a set of unordered points, sampled non-uniformly across an orebody in multiple cross-sections as input.

In [5], Sprague and de Kemp presented a partially automated tool which uses Bézier and NURBS (non-uniform rational B-spline) curves to model non-planar 3D surfaces. Fitting under the guidance of a control frame, it makes use of 2D interpreted plan views provided by mine geologists, created using local slices of projected drill hole data at semi-regular spaced depths. The authors emphasized that in moving from 2D polygonal lines to 3D surface construction, the continuation of common features requires delicate spatial synchronization across neighboring plan-sections, and manual correspondence performed by the geologist was critical to the integrity of the modeled structure and preventing self-intersections. The term map trace was defined as a “geological interpretation delineating the intersection of a geologic boundary as it breaks through the surface…or as it intersects a given elevation plane”. This definition closely resembles our notion of active contour extracted boundaries which represent segmented regions in a geozone. Their technique were refined by imposing positional and orientation constraints using structural ribbons (expert knowledge), thus it may be categorized as a semi-supervised technique.

In [6], Mitášová and Hofierka applied differential geometry, more specifically, regularized spline with tension to topographic analysis of a watershed. This utilized convex/concave sections with smooth curvature for interpolation. However, the focus was to model the top surface (as opposed to sub-surfaces) using dense elevation data obtained via remote sensing (as opposed to sparse data).

In [7], Kaufmann and Martin built 3D subsurface models using a variety of sources (drill holes, cross-sections and geological maps) with different motivations. Their goal was to further understand subsoil characteristics such as hydrogeologic or geothermic properties of the geological bodies. Their surfaces were modeled using DSI (discrete smooth interpolation) which computes the location of nodes by balancing roughness and misfit constraints (see Mallet [8] and Frank [2]). This is representative of a class of information-rich GIS fusion approaches which utilize topographic, geological and structural data [9]. In contrast, our approach makes the most of the limited data in an information-poor environment where only blasthole locations and geozone labels are available.

In [10], Caumon et al. presented general guidelines for creating a 3D structural model made of faults and horizons using sparse field data. Their focus was natural resource evaluation and hazard assessment; triangulated surfaces with variable resolution was also discussed. The authors offered many insights, one notable comment is that “3D subsurface modeling is generally not an end, but a means of improving data interpretation through visualization…to generate support for numerical simulations of complex phenomena (i.e., earthquakes, fluid transport) in which structures [11] play an important role.” This is also true of subsurface models serving as decision support tools in mining and geological exploration.

In [12], Dirstein et al. demonstrated that automated surface extraction and segmentation of peak and troughs from seismic survey can provide insight into structural and sedimentary morphology. In particular, differential geometry was used to select objects of concave and convex curvature, these features can help identify subtle cues for fluid flow events that are perhaps over-looked by conventional interpretation methods. For an in-depth survey of past efforts and current interests in 3D geological mapping, readers are referred to [13], [14], [15] and [16]. Relevant techniques from computational physics and computer vision are presented in Appendixes A and B. As a quick overview, the major themes and related works are highlighted in Table 1.

In both Dirstein’s and our work, curvature preserving spatial representations are used to good effect albeit the objectives and modalities are different. One aspect of Dirstein’s work is waveform analysis, where genetic fitness (similarity) relative to a genotype (waveform signature) can reveal relative stability of neighboring segments above and below a particular surface and confer understanding about its structure and stratigraphy. In our work, we obtain an essentially continuous representation for the slope of geozone subsurfaces from sparse boundary patterns through contour extraction and metamorphosis; this provides useful directions for subsequent drilling and exploration. Our motivation stems from an interest in integrating techniques from pattern analysis, computational physics, computer vision and signal processing, bringing these to bear on a subsurface reconstruction problem, solving it in unconventional ways.

The target application is the modeling of subsurface boundaries in an open-pit mine of sedimentary iron ore deposits.²²2The most common minerals are hematite $\text{Fe}^{3+}_{2}\text{O}_{3}$ and goethite Fe³⁺O(OH). These deposits are believed to have formed as chemical precipitates on the floor of shallow marine basins in a highly oxidizing environment during the Proterozoic eon, circa 1.9–2.4 billion years ago. Here, the sparse spatial patterns derive from blasthole samples where each is labeled with a geozone. The geozone labels provide a classification of geological domains based on mineralization and stratigraphic units [29]. From a mining perspective, these are used to segregate high-grade minerals from low grade or unprofitable materials. With some effort, the geozones can be turned into coherent spatial clusters. However, as unorganized points, or even with a discretized block-based representation, there is no easy way of obtaining a coherent motion field, one that describes the direction a boundary moves in consistently without artifacts or discontinuities. This is one major reason for using a curve-based, deformable surface model.

At the outset, it is important to distinguish this work from 3D segmentation approaches [30, 31] that use energy functional for minimum path optimization given contours at specific depths which is a common scenario in computer tomography. These approaches often require 3D partial derivatives to be computed, our data simply do not have the requisite (vertical) resolution to facilitate that. Furthermore, as our data are literally mined through manual processes, they only become available progressively at periodic intervals. In our application, slope trends need to be estimated in a causal manner using cross-section samples from as little as two successive benches (at two particular depths). Our problem also differs from grade estimation or material classification for which a variety of probabilistic inference and machine learning approaches are known [32], [33], [34], [12], [29], [35], [36], [37]. The focus here is explicit spatial representation which encompasses boundary extraction, region correspondence and capturing changes to subsurfaces using differential geometry.

The proposed system takes input in one of two forms. The contours which describe geospatial regions are either given or derived. Fig. 2(a) depicts a common mining scenario where stratigraphic information are obtained by examining core specimen extracted from the drilled holes. Alternatively, geozone transitions (marker bands) are established from geochemical assays or geophysical measurements taken while drilling. The main characteristics of this data is high resolution in z and sparse sampling in the x-y plane. Geologists typically interpolate the boundary at locations (dotted lines) between the drilled holes based on an understanding of the geological setting. The time and effort required to “join up” these vertical slices to form a preliminary 3D surface model can be substantial. Fig 2(b) depicts a second scenario where horizontal cross-sections are taken. This scenario differs from the previous through denser sampling in the x-y plane and having relatively low z-resolution; blastholes are plentiful albeit non-uniformly sampled. In this case, the contours are not given — only the coordinates and geozone designations (classification labels) of the blastholes are known. This second scenario poses additional challenges with respect to boundary extraction which will be further considered in this paper. These pictures highlight the main objective of this work which is to model how a volumetric region evolves through the cross-sections using contours and estimate the direction in which the boundary moves in.

The proposed boundary geometry modeling framework (henceforth, abbreviated as BGM) consists of four sub-systems. The processing pipeline is shown in Fig. 1. The role of each subsystem is briefly described below.

• •

  Boundary extraction finds connected regions (components) from labeled blastholes. It locates boundary samples on each component and converts them into a closed contour.

• •

  Spatial correspondence automates the process of region association and component alignment. Having identified one or more components on each cross-section, the next goal is to match-up and motion-compensate contours, i.e., find the optimal translation which brings corresponding regions into alignment. By convention, contours found in the top and bottom of any two successive slices are referred as ‘source’ and ‘target’ components respectively.

• •

  Contour metamorphosis uses particle trajectories to facilitate slope estimation. The goal is to model residual differences after source–target components are aligned in a common frame. Contour boundaries are embedded in a 2D level-set function, by tracking the movement of the zero-interface, one establishes pathways for morphing the source region(s) into the target region(s) as the level-set evolves under the dynamics of a PDE.

• •

  Contour prediction reconstructs a 3D volume through interpolation or extrapolation of normalized trajectories.

The main contributions are as follows.

• •

  For boundary extraction, an edge map synthesis procedure is described which converts sparse non-uniform data points to an image representation. Gradient vector field and active contours are used to extract shapes (regularized contours) from unordered edge pixels.

• •

  For spatial correspondence, region association and component alignment problems are solved using FFT cross-correlation under spatial constraints. Obstacle avoidance is approached from a resource allocation view point, multiple-source multiple-target component relationships are explored using intersection graphs.

• •

  For contour metamorphosis, slope estimation is achieved using particle trajectories. A curvelet backtracking algorithm has been devised to overcome tracking failure caused by branching; a phenomenon that leads to boundary distortion. This occurs when there is a significant mismatch in shape (local curvature) between the corresponding boundary segments.

Table 2 provides an overview of the techniques employed in each area. In the ensuing sections, the reasons for choosing these techniques will be elaborated as specific challenges are described. We discuss how standard approaches are modified to address particular needs. The solutions sought for different parts of the problem draw inspiration from different fields.

The modification involves using a kD-tree for nearest neighbor search and adopts a region growing approach. All samples from the same geozone and same cross-section are numbered from 1 to $n$, initially each belonging to the cluster of ‘self’. Neighboring samples within a radius of $r$ are merged with the current sample and labeled with the minimum cluster index amongst the group. Cluster membership information is propagated iteratively until no further changes occur and $S$ connected components remain.

For each connected component, boundary samples are identified by thresholding the local entropy which is significantly non-zero at geozone transition points. Suppose a sample $n$ has $N_{n}$ neighbors within a radius of $r$ and the fraction of samples belonging to geozones $g_{1}$ and $g_{2}$ are $p_{n,1}$ and $p_{n,2}$. The local entropy is computed as $h_{n}=-\sum_{i}p_{n,i}\log_{2}(p_{n,i}+\varepsilon)$. Sample $n$ is marked as a boundary sample if $h_{n}\geq\max\{T_{\text{entropy}},h^{\text{(median)}}_{n}\}$ where $T_{\text{entropy}}=0.5$ and $h^{\text{(median)}}_{n}$ [the median entropy in $n$’s neighborhood] is used to suppress “non-maximum” responses.

One limitation of this detector is the entropy tends to zero when the geozone labels no longer vary. This creates a problem as boundaries remain essentially open at the frontiers of surveyed regions which are not bordered by a different geozone. To remedy this situation, we perform orientation analysis to close these gaps. The objective is to recognize samples on the outskirts of a geozone as edges. The direction of the $K_{\text{orient}}$ closest neighbors from $n$ are computed and sorted in ascending order. If a gap larger than $T_{\text{orient}}$ radian is found, sample $n$ is deemed to be on an open edge that needs to be closed. Considering the blastholes are often sampled on a hexagonal lattice, we set $T_{\text{orient}}=\frac{2}{3}\pi$ and $K_{\text{orient}}=(4\times 2\pi/T_{\text{orient}})=12$. Some intermediate results are shown in Fig. 3.

One observation from Fig. 3(b) is that there is often a large gap between [what humans perceive as] adjacent samples along the [as yet undetermined] component boundaries. Indeed, the variable gap size renders useless boundary tracing algorithms like ‘marching squares’ that rely on fixed separating distances. A more robust approach to boundary extraction is to formulate it as an energy functional minimization problem that provides a tradeoff between smoothness and fidelity. For this, projection of the boundary onto an image grid constitutes the first step.

The main objective is to bridge the gap between adjacent samples located on the boundary. Each sample connects with its $K$ nearest neighbors by forming an edge between them. This edge structure comprising of line segments is then densely sampled and transferred over to the image plane. Specifically, edge energy is accumulated by pixels within some margin of each line segment.³³3Morphological closing is subsequently applied to give the edges adequate thickness. This procedure is further described in Algorithm 2.⁴⁴4Mapping real-world coordinates to the image domain (and vice-versa) involves scale changes. Uniform quantization is used implicitly throughout. For illustration, a synthesized edge map is shown in Fig. 4(a). Although it contains some false edges, the active-contour segmentation approach (to be described next) can tolerate small imperfection.

The final goal is to obtain a contour (polygon of regularized boundary points) given the synthesized edge structure. This is achieved by performing region segmentation using Active Contour Evolution (ACE) in a Gradient Vector Field (GVF). The GVF is obtained from the synthesized edge structure following the procedure given in Appendix A. It may be visualized as a field of arrows (an external force) that pushes any given point in space toward the boundary (see Fig. 4(b)). Accordingly, an active contour which has been initialized as the component bounding box will evolve over time under the action of the GVF and be drawn inward until it converges at the region’s boundary. This is illustrated in Fig. 4(c) and the short video (see multimedia content). The computational aspects of GVF active contour evolution are described in Appendix A. This highlights the first connection with computational physics.

Given two cross-sections with $n_{\text{s}}$ source components and $n_{\text{t}}$ target components, the goal is compute the association matrix $A\in\mathbb{R}^{n_{\text{s}}\times n_{\text{t}}}$ which describes the relationships between $s$ and $t$. In the preliminary assessment, a ‘1’ in $A(s,t)$ implies there is potential interaction between source $s$ and target region $t$.

Having obtained the association matrix $A(s,t)$, source–target relations can be represented by a graph. This gives rise to several possibilities: a 1-to-1 mapping, many-to-one, one-to-many, and many-to-many correspondence between the source and target components. These scenarios are shown in Fig. 6. For complex scenarios (e.g., multiple sources associated with multiple targets), the intersection graph is decomposed into subtrees (see Fig. 6 inset). This section is concerned with devising suitable strategies for each scenario, i.e., finding source translation estimates which maximize target alignment or component coverage. Graphically, each scenario is described by a subtree rooted at an s-node and its expansion includes all connected t-nodes and their children (any s-nodes connected to the t-nodes).

It is important to note that these level-set zero interfaces $\Omega_{t}$ are functions of time. Ultimately, we need to obtain iso-contours at different elevations $z$. However, no single $\Omega_{t}$ (other than $\Omega_{0}$) currently corresponds to a fixed elevation.

As a corollary, not all points on the advancing wavefront will reach the target contour at the same time since the wave propagates at different speeds depending on location.⁹⁹9The propagation speed depends on the geometric distance between the interface and target contour in accordance with (19). Therefore, particle trajectories (see Fig. 10(d)) are used to track the direction in which the zero-interface moves to facilitate slope estimation at a specific height. The relevant procedure called “particle advection” is described in Appendix B wherein (26) represents the key equation for particle update.

The particle speed (spacing) along each trajectory will eventually be adjusted to ensure each advances at a uniform rate and reaches the target boundary at the same time. This allows the zero interfaces to be converted into iso-contours of constant elevation. These ideas capture the essence of contour metamorphosis; Fig. 10(e) also serves as a prelude of the end goal.

In principle, Appendix B provides all the necessary tools for capturing information on local spatial deformation. In practice, when the shape of the source and target regions are dissimilar, the scheme may fail. Specifically, particle trajectories may fail to track the wavefront in areas where complex topological changes occur especially when they coincide with large differences in local curvature. This can be seen in Fig. 11(a) where some target boundary segments have zero particle coverage. Fig. 11(b) highlights regions which end up being unmodeled despite being reached by the wavefront. This usually occurs due to undersampling when the cross-section spacing is large. The goal in this section is to determine its root cause and contribute a solution to this problem.

In each time step, wavefront segments which the particle trajectories had failed to track are archived. This produces a collection of time-stamped pixels (called curvelets) which identify areas of innovation, viz., locations where the wavefront and head of the particle trajectories have diverged. The window of interest in Fig. 11(c) is magnified in (d) to reveal its structure. It can be seen the waves emanate from two separate sources and merge in the middle, this describes a region splitting situation which can be reasonably expected based on the topology shown in (a). For other unmodeled regions shown in (c), the situation is less complex and the behavior can be attributed to significant differences in local curvature. The common theme is that the missing particle trajectories appear to emerge from a single branching point.

The color strands in Fig. 11(e) show the trajectory flow for new particles sitting on the curvelets. The arrows indicate the direction the backtracking algorithm proceeds in. The aim is to establish pathways from the target back toward the branching point. In Fig. 11(f), new pathways deduced from the curvelets are slotted into the unmodeled regions. Evidently, this recovers directional information between the source and target contours in areas where information had previously been lost.

1. (a)

   Curvelet endpoint identification
   Locate the start and endpoints of prospective curvelets within $\mathcal{P}^{t}$ — the set of unordered pixels on the wavefront which have drifted away from the particle trajectories at time $t$.

2. (b)

   Curvelet segmentation (clustering and ordering)
   Extract individual curvelets $\zeta^{t}_{c}$ by tracing the connecting pixels between relevant endpoints.

3. (c)

   Curvelet regularization
   Fit pixels with a smoothing spline to remove jitter. Resample segment as $N_{c}$ point regularized curvelet $\kappa^{t}_{c}$.

4. (d)

   Curvelet track initialization / extension
   The displacement of curvelet points between successive time steps $[\tau^{t}_{c}(n),\tau^{t-1}_{c}(n),\ldots]$ describe possible pathways for the missing trajectories in unmodeled regions.

5. (e)

   Curvelet association
   Given regularized points $\kappa^{t+1}_{c}(m)$ for curvelet $c$ at time $t+1$, find a matching curvelet $k$ at time $t$.¹⁰¹⁰10Nearest neighbor search is employed to find a surjective mapping from $t$ to $t+1$. Each point $\kappa^{t+1}_{c}(m)$ has at least one match at $t$, but not every pixel on $\zeta^{t}_{k}$ is matched with a point at $t+1$. For this reason, as illustrated in Fig. 12(e), this is used only for curvelet association and not for point-wise correspondence. This is because inconsistent (jittery) movement has a detrimental effect on trajectory estimation. To eliminate the cumulative effects of drifting which would otherwise skew the tracks over time, the associated point set is constrained to lie on a regularized curvelet as shown in Fig. 12(f). A case of interest is depicted in Fig. 12(e) where the segmented pixels $\zeta^{t}_{k}(n)$ from two curvelets $k=A,B$ are merged from $t$ to $t+1$. This is handled by considering interval mappings $\mathbf{k}^{t+1}_{c,p}$ for multiple curvelet parts ($p=1,2$).

6. (f)

   Curvelet correspondence
   Once matching pixels $\zeta^{t}_{k}$ are found, the curvelet at time $t$ is smoothed and resampled evenly with $N_{c}$ points. This curvelet $\kappa^{t}_{k}$ is then collectively mapped onto $\kappa^{t+1}_{c}$, the current curvelet.¹¹¹¹11This provides an interesting contrast to the wavefront tracking approach proposed by Tomek et al. [49] which represents two sets of time-stamped pixels with a complete oriented bipartite graph, poses the task of finding the best tracking arrows as a network flow problem, formulates and solves it as an integer program using branch and bound. Similarly, Singh et al.  treated this as an optimization problem and employed dynamic programming in [50]. For a given point $n$, the sequence [$\kappa^{t+1}_{c}(n),\kappa^{t}_{c}(n),\kappa^{t-1}_{c}(n),\ldots$] subsequently defines a track which traverses the unmodeled region.

These steps capture the essence of the backtracking algorithm which produces particle tracks, see color strands in Fig. 11(e) and (f), for the unmodeled regions. The remaining tasks involve branching point identification which enables existing particle trajectories [the trunk portion extending from the source contour to the branching point] to be fused with the newly discovered tracks [from the branching point to the target contour] to complete the missing pieces. The details are described in Algorithms 11 and 12 and the notations are clarified in Fig. 13.

Spatial correspondence and metamorphosis were applied to the segmented regions; culminating in the contours seen in Fig. 14 (b and e). These contours are triangulated to produce the surfaces shown in Fig. 14 (c and f). These results demonstrate that differential geometry can model subterranean boundaries, handling topological changes (e.g., region merging) satisfactorily. The proposed techniques essentially transform irregularly-spaced points of low resolution into continuously deformable surfaces that represent boundaries. Fig. 15 reveals the structure of the two geozones. The two vantage points show that $g_{1}$ lies above $g_{2}$, such bedded planes (layered geological formations) are common for iron ore deposits at the mine site.

To assess the utility of slope information conveyed by the model, precision and recall rates are computed for predicted contours and compared against ground truth boundaries. The terrain extends from 690m down to 610m and is partitioned into disjoint 10m intervals (called benches) and processed top-down. Predictions are made using only information above the current bench floor. Therefore, predicted contours represent extrapolations from the known interval. For each comparison, data at the predicted depths from the bench below are withheld (not used to inform the prediction model) and these contours are used only for validation purpose at a given bench. Precision and recall rates are defined using set notations. The general setup is shown in Fig. 16. For predicted region $P_{i}$, precision is defined as $p(i)=\bigcup_{j\mid G_{j}\cap P_{i}}\left|P_{i}\cap G_{j}\right|/\left|P_{i}\right|$. For ground truth region $G_{i}$, recall is defined as $r(i)=\bigcup_{k\mid G_{i}\cap P_{k}}\left|G_{i}\cap P_{k}\right|/\left|G_{i}\right|$. The overall precision (resp., recall rate) at depth $d$ is the area-weighted average precision (resp., recall rate) of all contours at the specified depth. The depth varies from 0 to 10m in steps of 1.25m. Precision and recall statistics are computed under two conditions: (i) not using slope information [prediction employs zero-order hold], (ii) using the gradient field from the proposed model to obtain displacement estimates, i.e., combining the translation with the local deformation component estimated during spatial correspondence and contour metamorphosis, respectively.