gimli-org / gimli

Geophysical Inversion and Modeling Library :earth_africa:
https://www.pygimli.org
Other
344 stars 131 forks source link

Question about structural constrain inversion technique structural weighting adjustment #658

Open p56075607 opened 2 months ago

p56075607 commented 2 months ago

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, resist3ivity], [...]
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

Untitled

halbmy commented 2 months ago

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:

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 commented 2 months ago

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? image

halbmy commented 2 months ago

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 commented 1 month ago

Thank you very much for the useful advice! I found that the forward operator of ERT Manager has no setting for the setInterfaceConstraintso I turn into the ERTModellingclass 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 commented 1 month ago

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 commented 1 month ago

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 commented 1 month ago

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

halbmy commented 1 month ago

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 commented 2 weeks ago

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

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!

slope_synthetic_Ws1

halbmy commented 1 week ago

Great! Thank you for sharing this comparison.