Question about structural constrain inversion technique structural weighting adjustment

[p56075607]

Question about structural constrain inversion technique structural weighting adjustment

Problem description

I conducted a synthetic test of a slope model with a curved slip surface case with the structural constrain technique. Based on the example of the website about the structural constrains. I created a model with two layers: the regolith layer and the bedrock of the slope. Between these two layers is a curved slip interface.

According to (2011, Rucker) and (2020, Jiang), the settings for sharp boundary weighting suggest that the program should directly set the smoothness weighting of that boundary to 0, creating a very apparent resistivity contrast.

However, I'm curious if there's a way to manually adjust the structural weighting of the inversion smoothness constraint matrix $\bold{C}$ to a specific quantity, which could lead to a resistivity model with less "sharp" structural contrast. When we are not entirely confident in our structural interface from seismic or GPR.

Your environment

--------------------------------------------------------------------------------
  Date: Wed Feb 28 11:22:14 2024 CST

                OS : Darwin
            CPU(s) : 8
           Machine : x86_64
      Architecture : 64bit
               RAM : 8.0 GiB
       Environment : Jupyter
       File system : apfs

  Python 3.11.6 | packaged by conda-forge | (main, Oct  3 2023, 10:40:37)
  [Clang 15.0.7 ]

           pygimli : 1.4.5
            pgcore : 1.4.0
             numpy : 1.24.4
        matplotlib : 3.8.2
             scipy : 1.11.4
           IPython : 8.18.1
           pyvista : 0.43.0
--------------------------------------------------------------------------------

Steps to reproduce (my code)

# Build a three-layer model for electrical resistivity tomography synthetic test using pygimli package
import matplotlib.pyplot as plt
import numpy as np
from os.path import join

import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import ert

# Model setup
c1 = mt.createCircle(pos=(0, 310),radius=200, start=1.5*np.pi, end=1.7*np.pi,isClosed=True, marker = 2, area=1)
slope = mt.createPolygon([[0.0, 80], [0.0, 110], 
                          [c1.node(12).pos()[0], c1.node(12).pos()[1]], 
                          [c1.node(12).pos()[0], 80]],
                          isClosed=True)
geom = slope + c1
mesh = mt.createMesh(geom,area=1, quality=33)
pg.show(mesh, markers=True, showMesh=True)

# Synthetic data generation
electrode_x = np.linspace(start=0, stop=c1.node(12).pos()[0], num=25)
electrode_y = np.linspace(start=110, stop=c1.node(12).pos()[1], num=25)
# Plot the eletrode position
ax, _ = pg.show(slope, markers=False, showMesh=False)
ax.plot(electrode_x, electrode_y,'kv',label='Electrode')

scheme = ert.createData(elecs=np.column_stack((electrode_x, electrode_y)),
                           schemeName='dd')

# Create a map to set resistivity values in the appropriate regions
# [[regionNumber, resistivity], [regionNumber, resist３ivity], [...]
rhomap = [[1, 150.],
          [2, 50.]]

# Forward modelling
data = ert.simulate(mesh, scheme=scheme, res=rhomap, noiseLevel=1,
                    noiseAbs=1e-6, 
                    seed=1337) #seed : numpy.random seed for repeatable noise in synthetic experiments 

# Inversion using structural constrain mesh
c2 = mt.createCircle(pos=(0, 310),radius=200, start=1.53*np.pi, end=1.67*np.pi,isClosed=False, marker = 1)

plc = slope + c2
mesh3 = mt.createMesh(plc,
                      area=1,
                      quality=33)    # Quality mesh generation with no angles smaller than X degree
ax,_ = pg.show(mesh3)

# Creat the ERT Manager
mgr3 = ert.ERTManager(data)
inv3 = mgr3.invert(mesh=mesh3, lam=100, verbose=True)
mgr3.showResultAndFit(cMap='jet')

The results

[halbmy]

This is indeed a very good question of both practical and methodological relevance. As for any other data, structural information also bears uncertainty that should be accounted for. A few ideas:

• instead of putting 0 on the boundary one can set a different value, e.g. 0.2
• one could add several (sub-parallel) interfaces
• one could define an interface region where smoothness is decreased (reflectivity)
• a set of independent inversions with varying interfaces could be analysed statistically
• a separate region with geostatistical regularization resembling the interfaces

or any combination of it.

I'd be interested in other peoples experience, knowing about existing papers. Maybe we should make a discussion about this issue?

[p56075607]

Thank you for your helpful advice. Practically, if I want to set a structural weighting value of 0.2 to a specific boundary (in my case, the curve slip interface), how should I adjust the inversion settings?

[halbmy]

One can set such an interface weight of a forward operator by

fop.regionManager().setInterfaceConstraint(marker, weight)

(The marker should maybe have a marker different from 1 to be distinguished from the outer boundary.)

In case of a manager, the forward operator is mgr.fop, and one should set the mesh before calling the inversion:

mgr = ert.ERTManager(data)
mgr.setMesh(mesh)
mgr.fop.setInterfaceConstraint(2, 0.2)
mgr.invert(lam=100, verbose=True)

Please try it and report back, because I've never used it.

[p56075607]

Thank you very much for the useful advice! I found that the forward operator of ERT Manager has no setting for the setInterfaceConstraint so I turn into the ERTModelling class to use its regionManager() method to set Interface Constraint and run the inversion for testing different weighting values.

# Set such an interface weight of a FORWARD OPERATOR
fop = ert.ERTModelling()
fop.setMesh(mesh3)
fop.regionManager().setInterfaceConstraint(2, 0.5)
fop.setData(data)
inv3 = pg.Inversion(fop=fop, verbose=True)
transLog = pg.trans.TransLog()
inv3.modelTrans = transLog
inv3.dataTrans = transLog

# Run the inversion with the preset data. The Inversion mesh will be created
# with default settings.
constrained_model = inv3.run(data['rhoa'], data['err'],lam=100)

[p56075607]

And after the SEG webinar I saw the demonstration of the same inversion setting of the forward operator directly using ertManager. So both these two method can use to set the boundary of the inversion forward operator!

# The same inversion setting can be done with the ERTManager
mgr3 = ert.ERTManager(data)
mgr3.setMesh(mesh3)
mgr3.fop.regionManager().setInterfaceConstraint(2, 0.5)
mgr3.invert(lam=100, verbose=True)

[halbmy]

Exactly! So in this case it is not necessary to create the fop by yourself. In the next version, we will provide

mgr.inv.setInterfaceConstraint(2, 0.5)

just like

mgr.inv.setInterRegionConstraint(2, 3, 0.5)

[p56075607]

Thanks! I'm appreciate and look forward to it!

[halbmy]

Would be interesting to see a systematic study how the strength influences the results.

I just saw that your synthetic modelling includes severe modelling errors showing unrealistically low apparent resistivities for larger dipole separations. Reason are wrong boundary conditions. You should append a triangle boundary with "world boundary conditions" around your model to avoid this. There are several examples, just search for appendTriangleBoundary.

[p56075607]

Thank you for the advice for the appendTriangleBoundary mesh. The synthetic data results become more realistic!

And the change of the smoothness weighting results of this synthetic model are here, I adjust the $W_s$ from 0 to 1, and the main influences appear near the boundary!

[halbmy]

Great! Thank you for sharing this comparison.

[halbmy]

With the last version v1.5.1, both interface and interregion constraints can be set as above. Thank you again for raising and discussing this very interesting issue that helped getting insight and improving the code.