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.
Senior Geophysicist, Aarhus GeoSoftware
<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.
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
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.