3D Modeling of Large-Scale Geological Structures by Linear Combinations of Implicit Functions: Application to a Large Banded Iron Formation

Abstract

Implicit methods for modeling geological structures such as stratigraphy and faults have been developed for more than a decade, and they have made automatic model construction feasible. The implicit potential field method is such a method that is capable of incorporating multiple types of data including contact points for geological boundaries and their measured orientations. The implicit potential field method relies on the solution of a co-kriging system. However, applying the method to 3D modeling of large-scale geological structures constrained to dense data and strict geological rules remains challenging. Due to the non-stationary and complex nature of large-scale geological structures, and difficulty in estimating an adequate variogram model, performing global interpolation with all dense data together may create geologically unrealistic artifacts. We propose a framework that uses a divide-and-conquer strategy. The core idea is to create intermediate 3D geological models that match subsets of data and then recombine them into a single large 3D geological model, while maintaining data and geological rule constraints. We also prove that linear combinations of potential fields preserve properties of conditioning. The paper presents an application of the framework in modeling the stratigraphy model of a large banded iron formation in Western Australia with dense boreholes, but scarce orientation measurements.

Appendix

In the appendix, we prove that the conditioning to contact points and orientations is preserved with linear combination of two potential fields. Let $f\left( \varvec{x} \right)\varvec{~}$ and $g\left( \varvec{x} \right)$ denote two potential fields, $h\left( \varvec{x} \right)$ the linearly combined potential field and $w\left( \varvec{x} \right)$ the spatially varying weight.

Statement 1 The linear combination of two potential fields maintains fitting contact points that both potential fields already fit.

Proof 1 At contact points locations $\varvec{x}_{c}$, the linear combination of two potential fields $f\left( {\varvec{x}_{c} } \right)$ and $g\left( {\varvec{x}_{c} } \right)$ is:

$$h\left( {\varvec{x}_{c} } \right) = w\left( {\varvec{x}_{c} } \right)f\left( {\varvec{x}_{c} } \right) + \left( {1 - w\left( {\varvec{x}_{c} } \right)} \right)g\left( {\varvec{x}_{c} } \right)$$

Let $\tilde{f}$ and $\tilde{g}$ denote iso-values that fit $\varvec{x}_{c}$ for $f\left( {\varvec{x}_{c} } \right)\varvec{~}$ and $g\left( {\varvec{x}_{c} } \right)$, respectively. Let $\tilde{w}$ denote the weight at location $\varvec{x}_{c}$. Then, we have:

$$\begin{aligned} \tilde{f} & = f\left( {\varvec{x}_{c} } \right) \\ \tilde{g} & = g\left( {\varvec{x}_{c} } \right) \\ \tilde{w} & = w\left( {\varvec{x}_{c} } \right) \\ \end{aligned}$$

The statement “$h\left( \varvec{x} \right)$ fits contact points $\varvec{x}_{c}$”, entails “$\tilde{h} = h\left( {\varvec{x}_{c} } \right)$”$:\tilde{h}~$ is the iso-value in $h\left( \varvec{x} \right)$ that represents the boundary. We know that

$$\begin{aligned} h\left( {\varvec{x}_{\varvec{c}} } \right) & = w\left( {\varvec{x}_{c} } \right)f\left( {\varvec{x}_{c} } \right) + \left( {1 - w\left( {\varvec{x}_{c} } \right)} \right)g\left( {\varvec{x}_{c} } \right) \\ & = \tilde{w}\tilde{f} + \left( {1 - \tilde{w}} \right)\tilde{g} \\ \end{aligned}$$

so $\tilde{h} = h\left( {\varvec{x}_{c} } \right)$ if and only if:

$$\tilde{h} = \tilde{w}\tilde{f} + \left( {1 - \tilde{w}} \right)\tilde{g}$$

In other words, the iso-value used to represent the boundary in $h\left( x \right)$ should be the linear combination of iso-values in $f\left( \varvec{x} \right)~$ and $g\left( \varvec{x} \right)$, with the same weight used for combining $f\left( \varvec{x} \right)~$ and $g\left( \varvec{x} \right)$. To make sure new boundaries fit contact points, we not only need to do linear combination for two potential fields, but we also need to do a linear combination for the thresholding values and use it to create the categorical model. □

Statement 2 For any contact point, if only one potential field fits it, then the linear combination of two potential fields maintains fitting that contact point if the fitted potential field has a weight of 1 and the unfitted potential field has a weight of 0 at the contact point location.

Proof 2 Let us assume the potential field that fits $\varvec{x}_{c}$ is $f\left( {\varvec{x}_{c} } \right)$. As stated, $w\left( {\varvec{x}_{c} } \right) = 1$, so.

$$h\left( {\varvec{x}_{c} } \right) = f\left( {\varvec{x}_{c} } \right)$$

and

$$\tilde{h} = \tilde{f}$$

Therefore, $h\left( {\varvec{x}_{c} } \right)$ also fits $\varvec{x}_{c}$ as $f\left( {\varvec{x}_{c} } \right)$ does. This means that as long as we assign a weight of 1 to the fitted potential field at the contact point location, and a weight of 0 to the unfitted potential field, then our linearly combined potential field will fit to the contact point. □

Statement 3 The linear combination of two potential fields maintains fitting orientations that both potential fields already fit.

Proof 3 At first, we know the existing linearity of gradient theorem, which states that, for two scalar fields $f\left( \varvec{x} \right)$ and $g\left( \varvec{x} \right)$ at point $\varvec{x}_{o}$, we have.

$$\nabla \left( {\alpha f + \beta g} \right)\left( {\varvec{x}_{o} } \right) = \alpha \nabla f\left( {\varvec{x}_{o} } \right) + \beta \nabla g\left( {\varvec{x}_{o} } \right)$$

It is easy to see that if we replace $\alpha$ with $w\left( {\varvec{x}_{o} } \right)$, and $\beta$ with $1 - w\left( {\varvec{x}_{o} } \right)$, the linearity of gradient still holds.

Then, for orientation observation locations $\varvec{x}_{o}$ where $f\left( {\varvec{x}_{o} } \right)$ and $g\left( {\varvec{x}_{o} } \right)~$ have the same gradients, i.e., $\nabla f\left( {\varvec{x}_{o} } \right) = \nabla g\left( {\varvec{x}_{o} } \right)$, we can obtain the following using the linearity of gradient:

$$\begin{aligned} \nabla h\left( {\varvec{x}_{o} } \right) & = \nabla \left( {wf + \left( {1 - w} \right)g} \right)\left( {\varvec{x}_{o} } \right) \\ & = w\left( {\varvec{x}_{o} } \right)\nabla f\left( {\varvec{x}_{o} } \right) + \left( {1 - w\left( {\varvec{x}_{o} } \right)} \right)\nabla g\left( {\varvec{x}_{o} } \right) \\ & = w\left( {\varvec{x}_{o} } \right)\nabla f\left( {\varvec{x}_{o} } \right) + \left( {1 - w\left( {\varvec{x}_{o} } \right)} \right)\nabla f\left( {\varvec{x}_{o} } \right) \\ & = \nabla f\left( {\varvec{x}_{o} } \right) \\ \end{aligned}$$

In other words, the gradient remains the same after linear combination.

Therefore, as long as we use simulated orientations with the same observed “hard data” to generate $f\left( \varvec{x} \right)$ and $g\left( \varvec{x} \right)$ (i.e., $\nabla f\left( {\varvec{x}_{o} } \right) = \nabla g\left( {\varvec{x}_{o} } \right)$ for all orientation observation locations), all orientation observations will always be matched when doing linear combinations. □