Geologist in the Loop: A Hybrid Intelligence Model for Identifying Geological Boundaries from Augmented Ground Penetrating Radar

Abstract

Common industry practice means that geological or stratigraphic boundaries are estimated from exploration drill holes. While exploration holes provide opportunities for accurate data at a high resolution down the hole, their acquisition is cost-intensive, which can result in the number of holes drilled being reduced. In contrast, sampling with ground-penetrating radar (GPR) is cost-effective, non-destructive, and compact, allowing for denser, continuous data acquisition. One challenge with GPR data is the subjectivity and challenges associated with interpretation. This research presents a hybrid model of geologist and machine learning for the identification of geological boundaries in a lateritic deposit. This model allows for an auditable, probabilistic representation of geologists’ interpretations and can feed into exploration planning and optimising drill campaigns in terms of the density and location of holes.

1. Introduction

The accurate modelling of geological boundary transitions is important to ensure mine feasibility, with particular reference to maximising tonnages. Typically, these boundaries are estimated from exploration holes. The incorporation of ground-penetrating radar (GPR) into the modelling of boundary transitions is an attractive idea, as GPR is a non-destructive technology, and GPR data acquisition is relatively cost-effective compared to executing a drilling program for data acquisition. One of the challenges related to GPR data is that it does not inherently inform the geologist of the location of a boundary transition. The inference of transition locations from GPR data is subjective and open to geological interpretation. This is especially true in regions that have nondefinitive mineralogically extensive transition zones, where the change from one assay suite to another is a gradual one.

GPR is a non-destructive technology that can detect subterranean objects through the scattering of electromagnetic waves. It is a versatile technology that can be applied to a variety of fields, such as land mine detection [1,2], archaeology [3,4], the detection of buried people [5], civil engineering and construction [6,7], soil surveying [8], and mining and geology [9,10,11]. An overview of GPR and its applications in civil engineering is presented in reference [12]. A more broad analysis of early (1995–2014) GPR research and developments is presented in reference [13]. Advantages that GPR has over other non-destructive sampling methods such as seismic, acoustic, infrared, and optical are that it has a low survey cost, compact equipment, and an excellent trade-off between resolution and depth of detection [14].

This work introduces the idea of a “geologist in the loop” process that combines subject matter expert (SME) interpretations with Gaussian processes (GPs) [15], a machine learning algorithm, to represent these SME interpretations of the transition from one assay suite to another in a laterally extensive deposit. Specifically, this work is demonstrated in an extensive tropical region, a lateritic bauxite deposit in Cape York, Australia. The use of GPR, particularly in the pisolithic deposits of Cape York, has the potential to overcome the limitations of bauxite resource estimations from sparse exploration holes [16]. To obtain their interpretations, SMEs are provided with GPR data augmented with exploration hole data. The exploration hole data contains ground truth observations, such as the compositional percentages of all present assays. Of particular interest to this work is the percentage of SiO₂. Such data can be useful for pinning down interpretations of geological boundaries and removing potential ambiguities in the GPR domain. It should also be noted that, while the results presented in this paper are acquired from GPR data augmented with exploration hole data, the process can work on GPR data augmented with alternate data sources or simply on just GPR data.

SME interpretations are provided with some measure of uncertainty as an input into the GP model. The GP model then returns a mathematical representation of the location of the geological boundary. A GP model is ideal for this task, as the returned estimate is continuous, nonparametric, and probabilistic [15]. GPs have also been used in other geological tasks, such as the mapping of landslide susceptibility [17] and distinguishing geological units from measure while drilling data [18]. Having a mathematical representation of SME interpretations is attractive, as it allows for the incorporation of these data with other sources, such as exploration holes, gravimetric data, or aeromagnetic data, into further modelling processes. The inclusion of mathematical uncertainty in such a representation is also critical, as it allows for SMEs to highlight regions where they are unsure about where a boundary transition might lie. This is information traditionally not captured in the interpolation of boundary transition locations from GPR data. This mathematical representation of SME interpretations with uncertainty can influence mine planning decisions. One such example is the optimisation of drill hole locations, allowing mine site planners to best choose the location of exploration holes given the objective criteria and constraints related to their particular mine site.

While automated machine learning models have been applied to modelling from GPR data, such as the classification of a geological structure [9] or the estimation of the water content in soil [19], this paper presents a hybrid intelligence model for the identification of geological boundaries. Hybrid intelligence models employ a combination of human and machine intelligence, allowing for the strengths of each to be exploited [20]. The employment of a hybrid intelligence model in geological boundary identification is appealing, as it allows for the SME interpretations to be auditable and properly scrutinized. The review of a hybrid intelligence boundary estimate allows for an independent reviewer to not only understand the final boundary estimation but, also, to understand the thought process of the original geologist when they were generating the boundary estimate. This is especially true if both the SME inputs and boundary estimate outputs are represented with some uncertainty of measure.

The remainder of the paper is structured as follows: Section 2 explains the bauxite domain that we are using for the model. Section 3 provides a quick overview of how GPR works, with Section 4 following on how the data used in this paper was collected. Section 5 then describes the modelling process and how an interpretation of the boundary transition is generated. The paper then follows with Section 6 presenting the results and discussion, Section 7 covering the limitations and future works, and concludes in Section 8.

2. Geological Background and Profile of the Test Site

The bauxite province surrounding Weipa is found on the west coast of Cape York in Queensland, Australia, and the location of interest is noted (Figure 1). The laterally expansive bauxite region is roughly 14,000 km² and estimated to be around 50 Ma ([21,22,23]). Weipa bauxite deposits are noted to be substantially different to other bauxite deposits, so far as they mostly comprise “free” pisoliths, which are therefore not matrix-bound [21]. Further to this, there are also areas of pisoliths hosted in kaolin-rich horizons and areas of cementation [24].

Characteristic pisoliths (pisoid—singular) are described as the following: round to subrounded concreted particles of concentric mineralogies, created from the surrounding bauxite materials. All pisoliths vary in mineralogy and size, and therefore, no two are the same. Generally speaking, pisoliths are >2 mm in diameter, and particles below this size are known as oolites [23]. Figure 2 provides an example of pisoliths from the Weipa region.

The following are general minerals of bauxite deposits but may not be known at all bauxite deposits and may comprise many of the profiles described: böehmite (AlO(OH)), diaspore (AlO(OH)), gibbsite (Al(OH)₃), goethite (α-Fe³⁺O(OH)), hematite (Fe₂O₃), and kaolinite (Al₂(Si₂O₅)(OH)₄) [25]. The Weipa deposits range between 49% and 53% Al₂O₃ and are known to be some of the world’s highest grades [26].

Depth Profiling

The depth of the bauxite deposits is known to be around 15 m from the top to the ironstone basement, though this varies within the site [24]. The image shown in Figure 3 demonstrates a typical bauxite deposit regolith horizon in the area of interest for the GPR studies. Although individual horizons may vary in thickness over this particular bauxite deposit, the horizons are present throughout.

The lateral profiles forming boundaries within the bauxite deposits form via a chemical weathering process known as lateritisation and take place in tropical zones. Lateritisation is a leaching process [28] that takes place in high rainfall locations [23]. Tropical environments provide the catalyst for oxidation and the dissolution of minerals within an existing parent rock. Un-leached elements (immobile) remain as the constituents of soil. The lateritisation process leads to the oxidation of any mineral sulphides that may have been present in the parent rock material. In the right conditions, the relic soil left behind can reach thicknesses in excess of 10 m. The lateritisation process leads to the formation of relic soils rich in aluminium hydroxides, iron oxides, clays (kaolinite), and immobile elements such as titanium and zircon [23]. It is the high proportion of aluminium hydroxides that make this material sought after as a raw material for aluminium production. This process has taken place since the upper Cretaceous (approximately 100 million years b.p.) when the hot, humid climate resulted in extensive weathering of the Rolling Downs Group, forming deep lateritic profiles with gibbsite nodules and pisoliths in a kaolinitic matrix in the upper weathering horizon and oxyhydroxides and kaolinitic peds in the lower mottled zone. Uplift and erosion of the Rolling Downs Group resulted in the deposition of the Bulimba Formation, which continued to develop a new lateritic profile in situ [23].

3. Ground-Penetrating Radar

GPR is a geophysical technique that provides a reflectance from below the surface, which is displayed as a greyscale image to the end user. GPR uses propagating electromagnetic waves that respond to changes in the subsurface electromagnetic properties of the materials—in this case, bauxite. The dielectric permittivity, conductivity, and magnetic permeability of the shallow subsurface are the so-called responses [29,30,31]. The schematic “diagram” in Figure 4 provides a simple illustration of how GPR technology works.

Electromagnetic waves are radiated into the ground from a transmitting antenna at the surface. When propagating electromagnetic waves encounter any contrast in the electrical properties, they reflect off the interface of subsurface materials. The reflected signal is then received by a receiving antenna and recorded as a function of time [33].

4. Data Collection

There are two data sources used in this paper: exploration holes and GPR data. The methods for acquiring data relating to each of these sources is described in the remainder of this section.

4.1. Exploration Holes

The exploration holes are drilled using the Rotary Air Blast (RAB) technique and are drilled on a regular 50-m by 100-m grid. RAB uses compressed air to remove rock samples from the drill hole in a controlled manner. These samples are then sent to a laboratory for analysis, resulting in chemical and mineralogical information. A detailed explanation on drill cores in surface mines is provided in reference [34] for the interested reader. While several outputs of the data are supplied from the laboratory (assay compositions, hardness, etc.), expert geologists have indicated that SiO₂ is the primary dimension of interest for the identification of the bauxite transition. The incorporation of the other data outputs into this model is scoped out for future works.

4.2. Ground Penetrating Radar (GPR)

For this work, GPR data was collected using a frequency of 45 MHz. The choice of GPR frequency is dependent on the domain being modelled and the use of the data—for instance, a frequency of 50 MHz in a limestone deposit [35] or a lower frequency of 25 MHz for better subsurface resolution in a lateritic bauxite deposit [36]. The analysis of alternate frequencies at this Weipa site is for future works beyond the scope of this paper. The GPR data was collected along the same lines as the exploration holes at 50-m spacing between GPR lines. The lines are run in a north–south direction. Due to some environmental or heritage constraints, the GPR lines may deviate to accommodate this but will revert to north–south and in parallel with the adjacent GPR lines. To acquire the GPR images used in this work (and shown in this paper), raw GPR data was passed through a high-pass dewowing filter and amplified with an attenuation compensation gain [37]. Dewowing is used to eliminate low-frequency and DC bias in the data [38], while attenuation gain is used to help better visualize deeper targets [39].

5. Modelling Process

There are three main stages in determining a probabilistic boundary estimates from geologist-picked point data: GPR data representation, geologist point picking, and the training of the regression model from the geologist data. In the data representation step, a potentially augmented GPR image is presented to the geologist. It is this image that the geologist uses to identify boundary transitions. Note that non-GPR data is optional, but it allows the geologist to construct a more informed mental understanding of the transition. Exploration holes, or even another GPR image collected at a different frequency, are two examples of how an initial GPR image might be augmented. Once the geologist has selected the points that they have identified to sit on the boundary transition, the regression model can then be trained. This process of the geologist selecting a set of points and then training the model is an iterative one and can be repeated until the geologist is satisfied that the resulting model accurately represents their SME understanding of the boundary transition in a mathematical format. The employment of these three steps in this paper are described in the following subsections.

5.1. Data Representation

The GPR data for a “single line” is presented to the geologist as a two-dimensional image. This is a greyscale image that shows the radar response from different depths at different positions. We augmented this GPR with data from exploration holes that were “near” the GPR track. In this instance, we chose “near” to mean exploration holes that are within 2 m of the GPR track. An example image is shown in Figure 5. The exploration holes were coloured by relative SiO₂ values (Please note that, while these are the images that were shown to geologists, numerical information was abstracted away for commercial sensitivity reasons. Likewise, the spatial information was similarly altered).

5.2. Geologist Point Picking

From the provided image (Figure 5), the geologist selects a set of points that they believe lie on the boundary transition. For each of the chosen points, the geologist also provides a vertical measure of uncertainty, representing their understanding of where the vertical boundary transition might lie for the given horizontal coordinate. An example of this is shown in Figure 6. The geologist’s chosen points are presented with the augmented GPR image. Each of the geologist’s chosen points has an associated uncertainty, with points near the exploration holes having significantly less uncertainty due to the included exploration hole data. With these points and the associated uncertainties defined, a model can now be trained.

5.3. Machine Learning Model

Given the supplied geologist points, GPs are used to estimate the boundary transition. GPs are a Bayesian machine learning model that define a distribution over a function space, with inferences occurring in this function space [15]. Inferred functions from a GP model are continuous and use observation data to help describe the function performance in the “nearby” region. These model properties allow for GPs to model arbitrary, continuous spatial functions in a probabilistic manner, making them ideal for modelling geological boundary transitions.

A GP model requires a covariance function, a function that partially describes the relationship between the model inputs and outputs. The choice of a covariance function is defined by the user and defines the hyperparameters of the model. When training a model, the aim is to optimize these hyperparameter values to generate a model that best represents the training data. In this work, a squared exponential covariance function was used to spatially relate the input coordinates $x$ and $x^{\prime }$:

$$k\left(x,x^{\prime }\right)={\sigma }_f^2\ast \exp (-\frac{{\left(x-x^{\prime }\right)}^2}{l}),$$

(1)

where ${\sigma }_f^2$ is the variance, and $l$ is the characteristic length scale. Given this covariance function and a set of observation or training data, $T=\left(\left\{x_i\right\},\left\{y_i\right\}\right), i=1:N$ of $N$ input points, where $y_i{}_{}$ is the height of the observed boundary transition at location $x_i$, and the function of interest, $f\left(x\right)$, can be sampled at a point of interest, $x_{\ast }$, to infer the corresponding predictive distribution, $f\left(x_{\ast }\right)$.

To achieve a predictive distribution, the GP model maps the input space to the output space by placing a multivariate Gaussian distribution over the function variable space. The GP model can then be defined by its mean function $m\left(x\right)$ and covariance function $k\left(x,x^{\prime }\right)$: $f\left(x\right) ~ GP\left(m\left(x\right),k\left(x,x^{\prime }\right)\right)$. By denoting inference points of interest as $T=\left(\left\{x_{\ast i}\right\},\left\{y_{\ast i}\right\}\right), i=1:N_{\ast }$, the joint Gaussian distribution with a zero mean is defined as

$$\left[\begin{array}{c}f\\ f_{\ast }\end{array}\right]~ \mathcal{N}\left(0,\left[\begin{array}{cc}K\left(X,X\right)& K\left(X,X_{\ast }\right)\\ K\left(X_{\ast },X\right)& K\left(X_{\ast },X_{\ast }\right)\end{array}\right]\right),$$

(2)

where $f$ and $f_{\ast }$ are the noise-free values of the function. $\mathcal{N}\left(\mu ,\mathsf{\Sigma }\right)$ is a multivariate Gaussian distribution with a mean $\mu$ and covariance $\mathsf{\Sigma }$, and $K$ is the covariance matrix computed between all points in the set. For example, the covariance matrix for $K\left(X,X\right)$ is

$$K\left(X,X\right)=\left[\begin{array}{ccc}k\left(x_1,x_1\right)& \cdots & k\left(x_1,x_N\right)\\ ⋮& \ddots & ⋮\\ k\left(x_N,x_1\right)& \cdots & k\left(x_N,x_N\right)\end{array}\right].$$

(3)

By incorporating point-based uncertainty into the model with variance ${\sigma }^2$, the joint distribution of Equation (2) becomes

$$\left[\begin{array}{c}y\\ f_{\ast }\end{array}\right]~ \mathcal{N}\left(0,\left[\begin{array}{cc}K\left(X,X\right)+\left(\eta +{\sigma }^2\right)I& K\left(X,X_{\ast }\right)\\ K\left(X_{\ast },X\right)& K\left(X_{\ast },X_{\ast }\right)\end{array}\right]\right),$$

(4)

where $y$ is the set of noisy observations, and $\eta$ is a column vector of individual observation offsets. In this work, $\eta$ is the set of uncertainties that geologists supply with their observations.

The predictive distribution for the points being estimated can then be defined as

$$p\left(f_{\ast }|\left(X_{\ast },X\right),y\right)=\mathcal{N}\left({\mu }_{\ast },{\mathsf{\Sigma }}_{\ast }\right),$$

(5)

where the predictive mean and covariance are, respectively, calculated as

$${\mu }_{\ast }= K\left(X_{\ast },X\right){[K\left(X,X\right)+\left(\eta + {\sigma }^2\right)I]}^{-1}y$$

(6)

$${\mathsf{\Sigma }}_{\ast }=K\left(X_{\ast },X_{\ast }\right)-K\left(X_{\ast },X\right){[K\left(X,X\right)+\left(\eta + {\sigma }^2\right)I]}^{-1}K\left(X,X_{\ast }\right).$$

(7)

To learn the GP hyperparameters values that fit the given geologist observations, the log of the marginal likelihood with respect to $\eta$ is maximized:

$$\mathit{log}p\left(y|X,\theta \right)=-\frac{1}{2}{y}^T{\left[K\left(X,X\right)+\left(\eta + {\sigma }^2\right)I\right]}^{-1}\mathrm{y}-\frac{1}{2}\mathit{log}\left|K\left(X,X\right)+{\sigma }^{2}I\right|-\frac{N}{2}\mathit{log}2\pi .$$

(8)

From left to right, the terms in Equation (8) represent how well the model fits the data, a complexity penalty, and a normalization constant. Strategies for determining $\theta$ must be able to handle the nonconvex nature of the log of the marginal likelihood. Example optimization strategies to achieve this include: simulated annealing [40,41], stochastic gradient descent [42,43], and the limited memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) algorithm [44]. The LBFGS algorithm was used in this work. For further information on GPs, the interested reader can refer to references [15,45,46].

Figure 7 shows the resulting GP model estimate of the boundary location given the geologist input data shown in Figure 6. The black line represents the mean location of the boundary transition, with the shaded green region representing the model uncertainty. Once this model is generated, the process of choosing/adjusting geologist inputs and generating a boundary estimate can be iterated until the geologist is satisfied with the mathematical representation of their interpretation of the boundary transition.

6. Results and Discussion

In this section, we present multiple interpretations of two different sets of augmented GPR data. These examples demonstrate how interpretations can change from geologist to geologist, how different levels of SME can alter GPR interpretations, and how the final interpretations from each geologist can be audited by an independent source.

6.1. Location 1

The augmented GPR data for this location can be seen in Figure 5 (in the previous section). The points and associated uncertainties that each geologist chose are shown in Figure 6, Figure 8a and Figure 9a, with the respective resulting regression estimates of the boundary shown in Figure 7, Figure 8b and Figure 9b. A comparison of all regression estimates is shown in Figure 10.

It can be seen from these results that two geologists had similar interpretations, while the third interpretation varied towards the right of the GPR data. In review of the results, the third geologist’s interpretation of the transition on the far right side turned out to be cemented bauxite. Cemented bauxite occurs where lenses of pisoliths are fused together and typically requires a large amount of force to break up [24]. Cemented bauxite has a higher silica value than non-cemented bauxite, which results in the spike shown on the exploration hole. The change in silica values, along with the trend shown in the GPR data, can result in a completely automated boundary detection process choosing an erroneous boundary. This demonstrates how the incorporation of SME knowledge can be leveraged to improve model robustness. Further to this, the auditability of the hybrid modelling process allows for a thorough review of the third geologist’s interpretations, providing the opportunity to treat the mistake as a learning opportunity related to the presentation of cemented bauxite in GPR and exploration hole data.

6.2. Location 2

The GPR and exploration hole data for this location is shown in Figure 11. The observations made by the three geologists and their corresponding transition interpolations are, respectively, shown in Figure 12, Figure 13, Figure 14 and Figure 15, showing the three resulting geologist transition estimates overlaid. It is interesting to observe, in this instance, that the geologists have adopted differing observation strategies. The first geologist (Figure 12) has opted to place minimal observations on the right side, only marking near the exploration holes. This allows for the uncertainty of the model to “naturally” increase as it moves away from data. The second and third geologists (Figure 13 and Figure 14) generated observations that were equi-spaced across the entirety of the GPR data. Despite the differences in the modelling strategies, we can see that the two geologists generated similar interpretations that were within each other’s uncertainty profiles.

Another geologist interpretation is shown in Figure 16. Whilst the exploration holes were included for reference, they were not used in the generation of the geologist interpretation; only GPR data was used. The interpretation followed the stronger GPR signal and dipped under the exploration holes on the right. Comparing this result with the previous interpretations (Figure 15), we can see that the augmentation of the GPR data with the exploration hole data allowed for deeper SME insight, which resulted in a more informed transition interpretation. Similar to the previous example, it is likely that a completely automated model would have generated an inappropriate interpretation of the transition given that it would have lacked this depth of SME knowledge to draw upon.

7. Limitations and Future Works

While this work looks at the augmentation of GPR image data with exploration hole data, there are several related questions to explore. While the GPs were largely chosen in this work due to their ability to return probabilistic model estimates, it is of interest to see how they perform in other modelling methods based on seismic [47,48] and radar data [49,50]. Extended use of this work is going to be dependent on user requirements, and these requirements could help to dictate suitable alternate modelling techniques. Geologist interpretations of the data (and how accurate the geologist is) are going to be heavily influenced based on the data presented to them. We would therefore like to investigate how geologist responses change when different GPR data-processing methodologies are applied and when GPR data captured at different frequencies are used. It would also be of interest to present a wider variety of exploration hole-based information to the geologist, such as the full suite of assay compositions, or even measure-while-drilling data to see how this influences their interactions with the model. The consideration of several GPR samples, especially when they intersect with each other, is also something to investigate in future works.

We would also like to introduce this concept into 3D model generation and update the model with interpolation and extrapolation techniques generally associated with the geological visualisation of site-based data (e.g., reference [51]). GPR data and little (or no) exploration data could be used to generate an initial 3D model that could then be used to infer ideal places for future exploration holes to be drilled and then update the 3D model. This process could be repeated several times, with the specifics dependent on the time, cost, and resources available. The resulting process of exploration hole placement could be compared with alternate sampling patterns (e.g., reference [52]), and the final 3D estimate could also be compared to the estimates derived from other models (e.g., reference [53]).

A current limitation of this work is that the modelling is currently performed in individual 2D slices of GPR data. The consideration of nearby or intersecting GPR tracks is something that we would like to further explore in the future. The extension of the hybrid intelligence algorithm to non-lateritic deposits is also something to explore in the future.

8. Conclusions

In this paper, we presented a hybrid intelligence model that can generate probabilistic interpretations of geological transitions from GPR data augmented with exploration hole data. This method allowed for the uncertainties of geologists to be encapsulated into a mathematical model, resulting in objective and more well-informed transition interpretations. These interpretations can then be audited, allowing for a level of governance in the generation of these interpretations, and, if adopted, would present a change in the standard operating procedure for the industry. Finally, as these hybrids intelligence models are mathematical, they can be fed into consequential mine strategies such as exploration planning, mining optimisation, and the optimisation of drill campaigns in terms of the density and hole location.