Skip to main content

Lyceum 2021 | Together Towards Tomorrow

We will look at the setup and use of different model types for the inversion of TEM/FEM data in Aarhus Workbench.

Most inversions are done with smooth models only, but often the resolution is good enough to sense more abrupt changes in the geology than the smooth model can allow for. This can make the use of subsequent inversions with layered models interesting for complementary insights about the subsurface. We will look at layered models, block models and briefly sharp models, as well as an option to create individual layered models from existing smooth models.



Bjarke Roth
Senior Geophysicist, Aarhus GeoSoftware


17 min

See more Lyceum content

Lyceum 2021

Video Transcript

<v ->Hello, thanks for joining this presentation.</v>

My name is Bjarke Roth and I’m a geophysicist

with Aarhus GeoSoftware.

Today we are going to look at model types

as we use them for inversions in Aarhus Workbench.

Let’s say that we have measured some data,

it could be did TDEM soundings, like this one,

from an airborne time-domain EM system,

or it could be a frequency domain EM systems.

And we have spent some time to remove couplings

and noise from the data

so that it not only reflects the geology.

Before we can pass it on,

we need to turn it into a physical model of the subsurface

that a geologist then can use

to create the more elaborate models that they work on.

Unfortunately, there is no straightforward way to do that.

And in fact,

it turns out to be much easier to go the other way.

We know the physics and with an accurate system description

it’s possible to calculate the forward response

that one would measure for a given model.

So, what we can do is to compare this forward response

with the measured data.

We do this comparison with an objective function,

an example of that is this data residual.

Here we compared the observed data and the forward data

and we normalize it with the data uncertainty.

Really the C here is the diagonal entries

in the covarience matrix,

but it contains the data uncertainty squared.

A value of one

means that the data have fitted within the noise.

Anything below one is, in that sense, equally correct.

A value equal to two

means that data have fitted

with an interval of two times to noise and so on.

By minimizing this function,

we will find the version of our model

that comes closest to giving us the data we have observed.

The simplest 1D model we can have is a half space model,

that is a model with only one resistivity parameter

to describe the subsurface.

So let’s have a quick and somewhat simplified look

at how the inversion process works for the situation.

If we could calculate the data residual

for all allowed residual values it might look like this,

but inversion cannot see this curve.

Instead it will, for each iteration,

calculate the data residual

and then the derivative for the model parameter.

And then it can adjust the resistivity model parameter

based on that derivative taking a jump into subspace.

This repeats and repeats

until it cannot improve the result anymore.

And hopefully it ends up in the global minimum

rather than a local minimum.

This approach is called Levenberg-Marquardt algorithm

or damped least-square.

I have included a few references

for anyone who’s interested in those.

It doesn’t always find a global minimum,

but it is quite robust for where well behaved functions

and reasonable starting parameters.

And in most cases,

the functions we work with will be fairly well behaved.

You can, however, in some cases, end up in situations

where you need to improve your starting model

to not end up in a local minimum.

We, of course, usually want a little more detail

in our 1D models than just the one resistivity parameter.

Conceptually, the simplest model is the layered model

where we have a few layers,

each layer having a resistivity and a thickness.

Unfortunately, it is not easiest model type to work with.

The layers you end up with

will not always correspond directly

to the overall geological layers,

so it takes practice and understanding of subsurface

to get it right.

But we will talk more about why that is so in a bit.


we will turn to the smooth model as our starting point.

Here we have more layers, 20 to 40 depending on data type,

but only the resistivities is estimated.

The thicknesses are kept fixed

with increasing size with depth.

And the resistivities also vertical constrained

to ensure continuity between the layers,

hence to name smooth.

Let us unpack this a little.

There’s a big advantage to this approach.

The starting conditions are much less important

than for the layered models.

We of course have the starting resistivities

of all these layers,

but we almost always use a uniform starting model

to not introduce structure in the model with resistivities.

And as long as we don’t get too far away

from where we need to go with the resistivities,

this will just translate

into slightly longer inversion times

as the inversion needs to take a few extra iterations

to change resistivities and we do have some options

for improving that starting model,

those starting resistivities.

More importantly, because the layers are fixed,

the layer setup is not as important for the final model,

the final smooth model, as it is for a layered model.

All we need to decide for the setup is

the first layer boundary, the last layer boundary,

the number of layers, and the vertical constraints,

all of which we can do relatively objectively

based on the used instrument and data type.

The expected near surface resolution

give us the first layer boundary,

if we make this too thin,

it will just be determined by the starting value

and the constraints.

The deepest expected depth of penetration

give us the last layer boundary,

again, if we make this too deep,

it will just be determined by the starting value

and the constraints.

And we can sort that away

without depth of investigation measure later on.

But if you don’t make it deep enough,

we will not be able to fit the late term data.

The expected resolution give us an idea about the number of layers

to use, this also ties into the vertical constraints.

They are there to stabilize inversion

and reduce overshoot/undershoot effects

we might otherwise see with poorly resolved layers.

This can, however, again, be set relative objectively

based on the expected geological variation.

Then it’s simply a matter of entering the numbers

into a model discretization calculator

and let it distribute the layers

so that the increase thickness discretization

down through the model.

I talked about constraints,

but not really explained what they are.

All our constraints work like factors.

A value of two means that to resistivity of the layer above

or below are allowed one sigma interval on its resistivity

going from the layer of resistivity divided by two

to the layer of resistivity multiplied by two.

So it’s not a hard limit.

It can in principle change quicker than that,

but it’ll be more costly for inversion to do so.

Let’s have a look at an example.

The blue line is the true model

while the red line is estimated smooth model.

The DOI is the depth of investigation

that I briefly mentioned earlier,

which is the depth to which we can trust the results

to mainly depend on the data

rather than the starting values and constraints.

So the smooth model type reproduced the overall structure,

but we’re not exactly getting the correct layer transitions

and layer resistivities here.

Now let us compare that with estimate layered model.

Here we actually have two layered models,

one without a priori

and one where they resistivity of the third layer

was used to help the inversion in the right direction.

So clearly the layered model

can get much better layer transitions and resistivities,

particularly if we help it along.

If we, however, start the inversion with layer thickness

that are too far off,

it can end up placing the layer transitions differently

and wasting layers where we don’t really need them.

This situation can be improved

with the knowledge of the overall structure

that we can get from the smooth model,

but it can take a few attempts to get it right.

One can either use the first

and last layer boundary approach

or edit the layers individually,

but it tends to be more trial and error

than the smooth model set up,

perhaps because the layered model,

more so than a smooth model,

reflects the nativity of your data.

So make sure to focus your setup efforts on the top part

of the model where most of nativity of your data

will be reflected.

So looking at these two model types,

we are left with “all models are wrong, some are useful”.

Neither of these models hold all the answers,

but both models give useful

and somewhat complimentary insights about the subsurface.

Now let us broaden the approach a little

before we look at the last two model types.

We often invert a lot of models at the same time

so we of course need to talk about

how lateral constraints enter into this.

We expect a certain continuity in the geology.

It will depend on the specific geology of course,

but we expect models next to each other

to have some similarities.

We can either have constraints along the flight line,

this is called lateral constrained inversion or LCI.

Or we can have constraints both along the flight lines

and between the flight lines,

this is called a spatially constrained inversion or SCI.

Both of these are still 1D inversions,

the forward calculation here is still done in 1D.

But the models become quasi 2D for LCI and quasi 3D for SCI.

In practice,

the SCI constraints will look something like this.

The red points are the model positions

and the black lines are the constraints.

The constraints are determined with an algorithm,

it ensures that we get a structure like this,

where we always have constraints along the flight lines

as well as to nearby flight lines.

One important detail here is

that unlike with the vertical constraints

distance clearly needs to be considered here.

We don’t want the constraints to be equally strong

between all model positions.

So we set a reference distance for our inversion.

Any constraints between models

that stay within that distance, use the fact as given,

while the models further apart

will have that constraint scaled to become weaker.

If we set this reference distance

equal to the average sounding distance between our models

along the flight line,

the lateral constraint can be set

based on the expected amount of geological variation

in the area, rather than being some arbitrary factor.

The scaling will then take care of the rest of us.

For some data type

this reference distance is calculated automatically,

but, in other cases, it is part of the model setup.

Something like the altitude is also constrained,

but those constraints are kept separate

and only constraint along the flight lines.

Now we’re able to put all the pieces together.

When we talked about inversion earlier,

we used a simple objective function

that only looked at data residual.

The full objective function

most of course also take the constraints

and the a priori into account.

Just as the data residual compares observed and forward data

normalized with the data uncertainty,

the constraint parts

compares the constraint model parameters

normalized with the factor

set for the strength of those constraints.

And a priori part compares model parameters

with a priori for those model parameters

normalized with the factors

set for the strength of those a priori.

The constraint part is the most interesting here

as it highlights a few differences between the layered

and smooth models.

We have vertical constraints on resistivities

but only for smooth models.

We have lateral constraints on resistivities,

this is used by both layered and smooth models.

And we have lateral constraint on the layer thickness

of all layer depths, but only for layered models.

Both the thickness

and depth approach for the layered model works,

but there are, of course, some differences in setup

when we used factors on depth of all the layers above

rather than just the thickness of the individual layers

that is detailed however.

The constraint part is also interesting

because it is the basis for the two last model types.

The idea with both of them is to have a model

that you set up similar to how you set up a smooth model,

without needing to know much

about where the layer transitions should be,

but where the end result

takes on more of the characteristics of the layout result

with the more distinct layers.

The simplest approach is called a blocky model

and it is really surprising simple.

Rather than using the squared difference

it uses the absolute difference for the constraints part.

I have plotted the penalty cost of the sum terms here.

So what we are looking at is the penalty

given the absolute difference of the model parameters

normalized with the facts that we have set

for the constraints.

Clearly there two regions here.

When we are above one

the blocky model has a much lower penalty.

So it can, in that sense, allow much larger model changes

then the smooth model.

When we are below one

the blocky model has a slightly higher penalty,

so it can, in that sense,

not allow as large or small changes as the smooth model.

And that is exactly what we see.

Here we have a smooth model, a blocky model,

and a layered model inversion result

all of the same sounding.

The blocky model changes more,

changes the resistivity more quickly

and it stays more flat once it has made that change

than the smooth model.

It has the exact same setup as the smooth model

but gives us something that is a little closer

to the layered model in result.

In fact, the setup is so close to the smooth model setup

that we can take a copy of this, a smooth model SCI,

and just flip one switch to invert it as a blocky model.

Here, we have them not just for a single model

but over a line with 3.2 kilometers of SkyTEM data.

The layered model is struggling a bit

in a few transitions areas,

but otherwise the results look quite good

for all the model types.

The last approach is called a sharp model.

As you can see, for this it gets a little more complicated,

so I’ll not try to cover all the details.

I have again plotted the penalty costs of the sum terms,

or at least most of it,

as I’ll likely will be ignoring the beta value here.

This time the penalty stays closer to the smooth curve

for small values, but it fans out for larger values.

So once we get about one or so

it effectively becomes the number of variations

rather than the total amount of variation

that is added to the penalty cost.

This approach is called minimum gradient support.

In practice, we get a set of sharpness parameters,

one with a vertical sharpness

and one with a lateral sharpness.

These are related to this beta value I skipped over.

They influence how many blocks you get.

While they’re used for vertical and lateral constraints

affect the amount of variation within the blocks.

The number of blocks is, however, rather data-driven,

so it’s more like you set

how distinct it should make the blocks

than actually how many blocks.

Let us see some results.

Depending on the setup

the result can go from being very similar to a smooth model

to something like this,

where it ends up pretty close to the layered model.

It can change resistivities very quickly,

but, it is, of course, still limited by the resolution

of the fixed layers we have given it.

It is just like the smooth model in that regard,

but it certainly can come closer to the layered model

in result.

It does take some getting used to.

The numerical values for the sharpness parameters

are somewhat different from the constraints,

but there are suggestions for usable values on the F1 Wiki,

you can open from Aarhus Workbench.

Lateral values will make the blocks more distinct.

The usual vertical and lateral constraints

end up requiring smaller numerical values

than we are used to,

and the presets in Aarhus Workbench have been adjusted

to account for this,

but they otherwise behave exactly as expected.

Larger values will allow more variation

within each of the blocks.

The result I have here was found using the default setup,

so don’t be too discouraged

by the amount that this can be tuned.

Here we again have them not just for single model

but over a line with 3.2 kilometers of SkyTEM data.

In some areas,

the sharp model actually appears to be better

than the layered result,

and this was a fairly good layered result at that.

To summarize model types,

start with a smooth model inversion.

If you are looking to lock down a specific layer boundary,

layered models are your best bet, your best option,

as it is the only one that can actually move the layers,

but the more complicated geology gets,

the harder it becomes to do this.

If you are happy with the smooth model

but would like a more discreet version,

you have the blocky and sharp options.

Blocky uses the exact same setup as smooth,

so it’s very easy to try.

Sharp also uses the same base setup as smooth,

except for the added sharpness parameters

and a somewhat different vertical and lateral constraints

that may take a little while to get used to.

That may however be worth the effort

as it creates some very nice looking results.

That is all I have for you today.

Thank you for your time and have a nice day.