Constraining geophysical inversions with geologic information

Constraining geophysical inversions with geologic information Peter Lelièvre∗ , Douglas Oldenburg and Nicholas Williams, University of British Columbia Geophysical Inversion Facility SUMMARY Earth models used for mineral exploration should be reliable and consistent with all information available. The current focus of the Geophysical Inversion Facility at the University of British Columbia (UBC-GIF) is towards the development of a new generation of geophysical inversion codes and utilities to advance the integration of geologic and geophysical data through appropriate inversion methodologies. This research will provide more functional methods for applying geophysics to general mineral exploration problems. Here we outline some of the available types of geologic information that can be incorporated into UBC-GIF inversions and we provide an example that illustrates some of our methods. INTRODUCTION Geophysical inversion seeks to recover models of the Earth’s physical properties (e.g. density, conductivity) that can adequately reproduce anomalies in geophysical survey data (e.g. gravity, DC resistivity) while being consistent with geologic information. The physical properties are related to mineral structure and, hence, the models are an important source of information for understanding subsurface geology as it applies to mineral exploration. During the last several years the focus of the UBC-GIF group has been towards advancing the integration of geologic data (e.g. lithology, structure), physical property data (e.g. measurements on rock samples) and geophysical survey data through appropriate inversion methodologies. Procedures to incorporate geologic information into UBC-GIF inversions are being designed that exploit and push the boundaries of the existing functionality of the inversion codes. We are also producing some new technical methods that increase the kinds of geologic information that can be incorporated. These developments can help lead the practitioner towards the desired outcome of a common Earth model (Marquis et al., 2003), consistent with all sources of information, which can help answer specific geologic and exploration questions. TYPES OF GEOLOGIC INFORMATION AVAILABLE The geologic information available can come from many sources: surface mapping, drill-holes, hand samples, in situ measurements, preliminary mining or any other manner in which geologic data is collected. We categorize geologic information as either located or non-located. Located information, as the name suggests, is spatially tied. Possible examples of nonlocated information include expected target shape (e.g. an intrusive body should be disk-like or pipe-like) and that a partic- ular rock type is known to lie above another but the location of the contact is unknown. Below we list some of the many types of geologic information that can be, or have the potential to be, incorporated into UBC-GIF inversions. • Located information: – physical property measurements; – lithologic (rock type) observations; – structural orientations; – structural contacts between rock units. • Non-located information: – physical properties change sharply (e.g. across offset faults) or smoothly (e.g across zones of alteration) between rock units; – geostatistical information about the physical property values of each rock unit present; – relative positions of rock units; – physical properties increase or decrease in particular directions (e.g. density often increases with depth); – expected target shape and/or aspect ratio. An added consideration is that we should consider the reliability of any type of information included in an inversion, be it geophysical data or geologic information. For example, direct observations should be considered more reliable than interpolations, which should be considered more reliable than inferences or hypotheses. If the information is qualitative (e.g. expected shape of a target) it may be sound or speculative. Even if the information is quantitative (e.g. physical property measurements on rock samples) there may be a wide range of reliability reflected by the procedures taken to obtain that information. If a relative or absolute value can be placed on the reliability then this information should be incorporated along side the related geologic information. INCORPORATING GEOLOGIC INFORMATION INTO UBC-GIF INVERSIONS In the UBC-GIF numerical inverse solutions, the Earth region of interest is divided into many cells within an orthogonal mesh with the physical property of interest being constant across each cell. A model vector m holds the values in each discrete cell. Due to data uncertainty and other aspects inherent to the inverse problem, there are an infinite number of models that can fit the data to the desired degree. Further information (i.e. geologic information) is essential for a unique solution. The inverse problem is formulated as an optimization that involves minimization of an objective function Φ that combines a misfit measure Φd with a regularization measure Φm : min m Φ(m) = Φd (m) + β Φm (m). (1) Constraining geophysical inversions The tradeoff parameter β controls the relative size of the Φd and Φm measures for the resulting model and allows us to tune the level of data fit as desired. Li and Oldenburg (1996) designed a generic regularization function of the form Z
2 Φm (m) = ws m − mre f dv V  
2 ∂ dv m − mre f ∂x V   Z
2 ∂ + wy dv m − mre f ∂y V   Z
2 ∂ + wz dv m − mre f ∂z V + Z wx (2) which can be readily discretized on the inversion mesh. This regularization function makes the inverse problems tractable (it allows a single model to be recovered) and it allows incorporation of much geologic information into the inversion. The first term in (2) is a smallness or closeness term. In the case of a zero-valued reference model mre f this term encourages the inversion to recover models with low values of the physical property. When a reference model is incorporated, the inversion attempts to match it as closely as possible (while still fitting the data to the desired degree). The latter terms in (2) involving derivatives are smoothness terms that encourage recovery of spatially smooth models. The reference model may or may not be included in the smoothness terms depending on the application. In the absence of more specific, constraining geologic information we choose to seek a smoothly varying model that does not contain unreasonably high values. We recognize that this “minimum structure” inversion is just one of several options, others including “compact” inversions (Last and Kubik, 1983), “focussed” inversions (Portniaguine and Zhdanov, 1999) and non-smooth inversions (Farquharson and Oldenburg, 1998). However, we consider the minimum structure option here due to its computational simplicity. The UBC-GIF codes allow bounds to be placed on the property values within specific model cells. This is accomplished by adding bound constraints to the optimization problem: min Φ(m) = Φd (m) + β Φm (m) (3a) s.t. Li ≤ mi ≤ Ui (3b) m for some or all i where Li and Ui are the lower and upper bounds for the ith model cell. Li and Oldenburg (2003) discuss one possible method of solution to the bounded inverse problem, the logarithmic barrier approach, although more robust and efficient gradientprojection methods are currently being developed. An extension to these simple bound constraints under development is to allow for general linear constraints: min Φ(m) = Φd (m) + β Φm (m) (4a) s.t. Am ≥ b. (4b) m This functionality will allow bounding of both physical property values and spatial model gradients so that one can specify directional changes (trends) in the properties across the model. Much of the located geologic information listed above can be incorporated into UBC-GIF inversions in a natural way using the reference models and weights in the regularization measure in (2) or through the constraints in (3) and (4). Williams (2006) uses a simple real-world-based synthetic example and comprehensively demonstrates how this can be done to provide far superior results than default, geologically unconstrained inversions. Phillips (2001) and Farquharson et al. (2008) used UBC-GIF inversions to incorporate located geologic information into their inversion investigations of the San Nicolás copperzinc deposit and the Voiseys Bay nickel copper-cobalt deposit respectively. Farquharson et al. (2008) state that “incorporating geological information significantly enhanced the fidelity and geological consistency of the constructed models”. Some non-located information can be incorporated into the optimization problem using the same techniques but when this is not possible iterative methods can be developed to introduce this more complicated information. AN ILLUSTRATIVE EXAMPLE There are several commonly encountered inversion scenarios that deal with incorporation of some of the geologic information introduced above. These can be simplified into canonical synthetics such as a block in a half-space, a dipping block, finding the location of a sharp interface between two know units, and finding locations of anomalous material within a confined volume when other parts of the model domain are fixed. In the synthetic example below we combine some of these canonical problems into one to illustrate our methods for incorporating some types of geologic information. A dipping block with cover In this example we create a synthetic density model, calculate the gravity response above it and add a small amount of random noise to this data (shown at top in Figure 1) before inverting. The true Earth model, shown at bottom in Figure 1, contains a dipping target unit with anomalous density +0.6. The target reaches the surface and to either side lies a 30m layer of cover with anomalous density -0.5. The rest of the material contains the background density of 0.0. We assume that we have the following information with which to constrain the inversions: • there are three possible rock units: a bedrock unit (background density), a lower density cover layer, and a higher density target unit; • we have geostatistical information regarding the density distributions of those rock units; • the surface outcrop is known from surface mapping; • we expect the cover unit is horizontal and is no thicker than 50m; • we have some structural measurements taken at the surface that indicate a possible dip of the anomalous high density body. Constraining geophysical inversions indicated and the dip of the target unit apparent. Figure 3: The inversion result with the top layer of cells bounded by surface outcrop information. Incorporating orientation information Figure 1: The true density model (g/cm3 ) and gravity data response on a profile above it. Default inversion A default (unconstrained) inversion result is shown in Figure 2. Without further geologic information, this would be our best guess at the subsurface density distribution. However, the result is inconsistent with much of our geologic information: density values lie across a continuous scale rather than fitting into three separated narrow ranges and the surface of the model does not match the known information. Figure 2: The unconstrained inversion result. Superimposed white lines indicate the rock unit boundaries in the true model. Incorporating surface outcrop information To incorporate the surface outcrop information we use bounds on the top layer of cells in the model. The geostatistical information would provide a density mean and standard deviation for each unit; the bounds can be set around the mean values and the tightness of the bounds determined by the standard deviations. Here, we set the density bounds to [-0.52,-0.48] where the cover unit is observed and to [0.58,0.62] where the high density unit is observed. The recovered model with the surface outcrop information incorporated is shown in Figure 3. The result is immediately improved, with the cover layer now By altering the relative values of the smoothness weights wx , wy and wz in (2) one can cause the recovered models to become smoother in some mesh-orthogonal direction(s) compared to the other(s). The smoothness weights can be homogeneous across the entire mesh or can be set to different values in different regions; we can thereby specify orientations globally or locally. A further generalization of the objective function allows the coordinate axes to be rotated such that an elongation (i.e. a preferred direction of smoothness) can be specified in any (perhaps non-axial) direction; we do not show the form of this further generalization but invite the reader to refer to Li and Oldenburg (2000). In our example, we expect that the cover unit is horizontal and is no thicker than 50m. Hence, we set the wz /wx ratio below the default value of unity (we use 0.2 here) down to 50m to encourage features that are elongated in the x-direction (Easting) and may take larger jumps in the z-direction (vertical). To incorporate the structural orientation information regarding the dip of the high density unit we need to rotate the coordinate system such that the new x-axis (which we will denote as x0 ) lies along the expected dip direction. Hence, in a region immediately around the outcrop of the target unit we specify a dip of 45◦ . Outside of this region we do not rotate the coordinate system. We set the smoothness weights high where this information is most reliable: at the surface. Reliability estimates and knowledge of the scale of geologic information can be used to guide to how far sparse information can be extrapolated to influence its surroundings and thereby populate a larger region of the inversion mesh. Here we are less sure about the dip of the target unit below the surface so within the rotated region we increase the ratio w0z /w0x from 0.2 at the surface to the default value of unity at a depth of 50m. The recovered model with this new information incorporated is shown in Figure 4, which shows further improvement, especially in the depth extents of the units. Specifying three rock units In this final stage, we wish to better constrain the density to lie within three separated narrow ranges. As is evident in the previous results, the use of 2-norms in (2) results in recovered Constraining geophysical inversions cover layer and the target unit have been improved. Figure 4: The inversion result with surface outcrop and orientation information included. models that exhibit a smooth, smeared-out appearance. Sharp boundaries can be generated using Ekblom or Huber norms, for example as implemented by Farquharson and Oldenburg (1998), but here we explore how well we can do with 2-norms by employing an iterative strategy. We first take the model in Figure 4 and set the values to those of their closest rock types: i.e. any cells with density below -0.25 are set to the value of the cover (=-0.5), any cells with density above +0.3 are set to the value of the target (=0.6) and all other cells are set to the value of the background (0.0). The resulting binned model is shown in Figure 5. Figure 5: The binned model derived from the result in Figure 4. This binned model does not fit the data as well as desired so we can not stop here. We have less confidence in the binned rock model near the contacts between rock units. Hence, we continue by using the binned model as a reference model in the smallness term only and we set the smallness weights low near the contacts in the binned model and high away from those contacts. We also set the smoothness weights low near the contacts in the binned model to encourage the inversion to place sharper density gradients (i.e. jumps) at those locations and maintain smoothness elsewhere (within the rock units). Again, we increase the weighting values back to their default values as we move away from the contacts in the binned model. Doing so allows the inversion the flexibility to ignore our estimated contact locations if that is needed to fit the data. This procedure can be carried out in an iterative manner, again to avoid specifying the estimated contact locations too forcefully. The result after three iterations is shown in Figure 6. This final model is a much better representation of the true model than that for the default (geologically unconstrained) inversion in Figure 2: the three rock units are clearly defined and the depth extents of the View publication stats Figure 6: The inversion result after encouraging the density values to lie within separated narrow ranges as expected from the three rock units present. CONCLUSION When not constrained by geologic information, default UBCGIF inversions can generate reasonable results, recovering spatially simple physical property distributions that honour the survey data. However, such first-pass results may not honour the geologic information, which can be incorporated to drastically improve the results, as is evident from the example above and in the work of Williams (2006) and Farquharson et al. (2008). The UBC-GIF inversion codes contain the functionality to incorporate many types of geologic information and this functionality has shown great promise for practical use. However, by leaving this process to the end user some of this functionality has yet to reach its full potential, due in part to a lack of ease of incorporating the related geologic information. We recognize that the development of a software library that provides a suite of reliable and efficient utilities to ease this process would be highly beneficial. The vision of this software is for it to take files containing geologic information in native formats and use that information to build the required reference models, weights and bounds required for input into the UBCGIF inversion codes. Work-flows can also be documented to demonstrate methods and indicate best-practices. Research into incorporating geologic information into UBCGIF inversions is providing more functional methods for applying geophysics to general mineral exploration problems. We are continually discovering more ways to use and improve our inversion methods as they are used in real-world scenarios. ACKNOWLEDGMENTS We thank Richard Lane (Geoscience Australia) and Nigel Phillips (Mira Geoscience, previously UBC-GIF) for their contributions to this research.