Inversion of potential fields by interactive optimization of 3D subsurface models using a spring-based space warping and evolution strategy

Investigation of the Earth’s subsurface is important for sustained evaluation of natural resources and plays a key role in building digital subsurface models. We show how measured potential fields – gravity and magnetics – can be used to model subsurface properties in an elegant interactive way based on numerical optimization.

Abstract

The rapid development of dynamic data-integrative modeling of geologic processes and subsurface structures is an important factor for the sustainable utilization of natural resources. One of the current gaps in data-integrative modeling is the interactive construction of realistic 3D models of the earth’s underground. We have developed here a methodology of 3D interactive processing of potential fields — gravity and magnetics — as well as their potentials and derivatives, combining forward and inverse modeling. Forward computations are based on the approximation of geologic subsurface structures by polyhedra with triangulated surfaces. Inverse computations for the model geometry are performed by means of the covariance matrix adaptation evolution strategy (CMA-ES), which is proved to be efficient in the case of a strongly nonlinear problem and high-dimensional parameter space, as in potential field models. The main disadvantage is related to triangulation, as certain algorithmic constraints must be applied during the approximation of geologic body shapes. To avoid topology distortions, we use a concept of warping the space containing the model, rather than the model itself. However, the optimized lengths of grid sides are dependent on each other, which degrades the self-adaptation of the CMA-ES. The elegant solution is to introduce a system of virtual elastic springs connecting the grid nodes. We develop a numerical formulation of this system and provide a proof of the concept rather than an overview of the theoretical concepts of the inversion schemes. The new workflow is tested on the 3D SEG Advanced Modeling Program model and applied to a real case study of salt dome modeling in the Northwest German Basin. In the context of the inversion procedure described here, we determine how an interpreter can visually control and influence the quality of the inversion on a timeline.

Full text

The full text is unfortunately not open-access, but it is available for institutions through GoeScienceWorld using the following DOI: 10.1190/geo2022-0222.1

Background

Here we present some background ideas behind the idea of geometry optimization using potential fields.

First of all, detailed 3-D subsurface models are very important for understanding geohazards as well as for sustainable use of georesources. Creation of a 3-D model requires integration of the whole spectrum of input structural, geological and geophysical data, like geology maps, seismic images, well data, thickness maps and so on, including other 3-D models.

A good example of such a data-constrained 3-D model is the Upper Rhine Graben model:

The key concept to solve arising inconsistencies in the data is 3-D gravity modelling. For that purpose we recommend to utilize IGMAS+, as it can integrate various data in a comfortable way. It allows adjusting the geometries and physical properties of modelled subsurface bodies interactively, while the corresponding calculated and measured potential field components are visualized together with other independent observations.

However, even gravity-constrained 3-D models contain uncertainties which are a caused by gaps in the input data. Apart from that, manual geometry modifications are very time consuming, and model changes made this way are not reproducible.

In order to improve the modelling consistency and accuracy, we propose a semi-automatic optimization concept utilizing potential fields.

Concept

Non-linearity

First, as we consider optimization of geometrical parameters of a 3-D model, and it is a highly non-linear problem, this requires a suitable nonlinear optimization method.

Evolution Strategy

For optimization we have chosen the Evolution Strategy with Covariance Matrix Adaptation, or CMA-ES (Hansen and Ostermeier, 1997; Hansen, 2016; Weng, 2019).

The performance of the evolution strategy has been numerically compared with other optimization techniques and has been verified to be applicable in the context of an interactive modeling process, as shown by Alvers (1998):

The CMA-ES works well in high-dimensional space, demonstrates fast convergence and is able to overcome local minima.

In this simple problem on the figure above, the population concentrates over the global optimum within a few generations, but for more complicated problems it is not obvious.

One of the difficulties in CMA-ES application is that one should take care of the validity of generated models in terms of topology. To solve this we propose to use space warping complemented with a system of virtual elastic springs.

Space warping

The main idea behind the space warping is to distort the space containing the triangulated model, but not the model itself.

In other words, we propose to put a grid with a certain resolution over the model and then change the grid in such a way that the model becomes optimal, without dealing with the model vertices.

Optimality here is understood in the sense of a quality function represented by a misfit between the measured and the modeled potential fields:

Variation of the grid is done by optimization of the grid legs. If the grid legs are kept bigger than zero, we are able to preserve the topology of the original model. Derivation of the positions of model vertices is done by linear interpolation with respect to the distorted grid edges:

The main disadvantage here is that it is still necessary to check whether the topology is preserved so that the evolution strategy hasn’t introduced non-physical intersections of the triangulated model interfaces.

In the end, we need a space warping method that does not have these disadvantages.

Elastic springs

To avoid disadvantages of space warping using optimization of grid legs, and to keep the model topology, we suggest an elegant solution: to replace the rigid connections in the space grid with elastic springs.

A spring is characterized by two parameters: spring constant $k$ – a factor describing the stiffness of the spring, and its rest length $x_0$, i.e. the length of the spring before any deformation.

Their relation are explained by the Hooke’s law: $$F = k \Delta x$$ which states that the force ($F$) needed to extend or compress a spring by some distance ($\Delta x = x_1 - x_0$) scales linearly with respect to that distance. $x_1$ is in that case the length of the spring after the deformation.

It turns out that if we update the spring constants instead of the grid legs during the mutation step of the evolution strategy, we are able to derive a linear system to update the grid node positions, at the same time keeping the topology model consistent.

Visualization

The workflow presented in the paper is tested on a synthetic model (SEAM) and a real case study of a salt dome in the Northwest German Basin. Both cases are visualized below.

3-D space warping of the salt overhang in the SEAM model

The video shows space warping for the SEAM Phase I project model in 3D:

The 3-D view shows the salt overhang with a dynamically warped spatial grid. The map view on the right side shows the color-scaled residual field of the measured and modeled vertical gravity gradient changing with each iteration.

The SEAM (SEG Advanced Modeling) corporation created a model for deepwater Gulf of Mexico in a region containing salt and conducting variable-density acoustic seismic simulations on the model (see more in the SEAM Phase I documentation).

Optimization of the Büsum salt dome

The video shows the optimization of the Büsum salt dome situated on the North Sea coast of Schleswig-Holstein, Northern Germany.

The model, originally based on the representations in the Geotectonic Atlas of Northern Germany (Baldschuhn, 2001), represents results of the 3-D modeling of the three salt domes: Oldenswort (left), Büsum (center), and Meldorf (right). The upper panel in the video shows the 3-D model and the gravity field in a perspective view (left) and a vertical cross section along the gray plane (right). The solid red curve marks the measured gravity field, and the stippled curve depicts the calculated field which is used in optimization. Lower panel display the optimization status with three maps: the measured gravity (left), the modeled gravity (center), and the residual field (right).

Implementation in IGMAS+

The plugin for the spring-based geometry optimization is currently in the testing phase and is not yet available for the users. We plan to include it into one of our next official releases, together with the detailed documentation for it.

References

• Baldschuhn (2001). Geotektonischer atlas von Nordwest-Deutschland und dem Deutschen Nordsee-sektor Strukturen, Strukturentwicklung, Paläogeographie: Schweizerbart. (Bilingual digital atlas of the geotectonic structure of Northwestern Germany and the German North-Sea sector). ISBN 978-3-510-95881-8.
• Hansen, N., and A. Ostermeier (1997). Convergence properties of evolution strategies with the derandomized covariance matrix adaptation: The (μ=μ−l,λ)-CMA-ES. Proceedings of the Fifth European Congress on Intelligent Techniques and Soft Computing, 650–654.
• Hansen, N. (2016). The CMA evolution strategy: A tutorial. arXiv preprint. DOI: 10.48550/arXiv.1604.00772.
• Weng, L. (2019). Evolution strategies, https://lilianweng.github.io/posts/2019-09-05-evolution-strategies

Materials

• Inversion demo videos