Coupled DC-Magnetic Inversion

[pi1379]

Problem description

Hello, I'm trying to do a coupled inversion with DC and magnetic data based on the inversion script to Ronczka, et al. (2017). "Electric resistivity and seismic refraction tomography: a challenging joint underwater survey at Äspö Hard Rock Laboratory". But it seems that some functions are outdated or based on outdated pybert functions.

My environment

Python 3.11.8 | packaged by conda-forge | (main, Feb 16 2024, 20:40:50) [MSC v.1937 64 bit (AMD64)]

       pygimli : 1.5.1
        pgcore : 1.5.0
         numpy : 1.26.4
    matplotlib : 3.8.4
         scipy : 1.13.0
          tqdm : 4.66.4
       IPython : 8.20.0
       pyvista : 0.43.7

My Code so far

I run this test script with the gallery3d.dat Example data and synthetic magnetic data based on the ert-inversion results from this example. mag_synth.txt

import pygimli as pg
import numpy as np
import pygimli.meshtools as mt
from pygimli.physics import ert
from pygimli.physics.gravimetry import MagneticsModelling

#=LOAD DATA====================================================================
# DC
data_dc = pg.getExampleData("ert/gallery3d.dat")
data_dc["k"] = ert.geometricFactors(data_dc, dim=3)
data_dc["err"] = ert.estimateError(data_dc, relativeError=0.02)
print(data_dc)

# Magnetic
coord_mag = np.loadtxt('mag_synth.txt', delimiter='\t', usecols=(0, 1, 2), skiprows=1, dtype='float', unpack=True, encoding='latin1')
points = coord_mag.T

data_mag = np.loadtxt('mag_synth.txt', delimiter='\t', usecols=3, skiprows=1, dtype='float', unpack=True, encoding='latin1')
abserror_mag = 0.2 # nT

#=3D MESH GENERATION===========================================================
# create 3d mesh
plc = mt.createParaMeshPLC3D(data_dc, paraDepth=12, paraMaxCellSize=3,
                             surfaceMeshQuality=34)
mesh = mt.createMesh(plc, quality=1.3)
print(mesh)
#_ = pg.show(mesh, markers=True, showMesh=True, color='white')

# only use inner mesh for magnetic inversion, otherwise the jacobi matrix has the wrong size
mesh_mag = pg.Mesh()
mesh_mag.createMeshByMarker(mesh,2)
mesh_mag.createNeighborInfos()
#_ = pg.show(mesh_mag, markers=True, showMesh=True, color='white')

#=MAGNETIC IGRF, DEPTH WEIGHTING===============================================
# calculate igrf, set component
F, I, D = 50000, 75, 25  # total field in nT
H = F * np.cos(np.deg2rad(I))
X = H * np.cos(np.deg2rad(D))
Y = H * np.sin(np.deg2rad(D))
Z = F * np.sin(np.deg2rad(I))
igrf = [D, I, H, X, Y, Z, F]
cmp = ["TFA"]

# depth weighting according to Li, Y. & Oldenburg, D. (1996): 3-D inversion of magnetic data. Geophysics 61(2), 394-40
bz = np.array([abs(b.center().z()) for b in mesh_mag.boundaries() if not b.outside()])
z0 = 25
wz = 100 / (bz+z0)**1.5
print(min(wz), max(wz))

#=INVERSION====================================================================
# variables
lambda_dc=30.
lambda_mag=30.
a=0.1
b=0.1
c=1.0
max_iter=10
sep_iter=2
verbose=False

# start inversion class
mgr = ert.ERTManager(data_dc)
fop_mag = MagneticsModelling(mesh=mesh_mag, points=points, cmp=cmp, igrf=igrf)
mag = pg.Inversion(fop=fop_mag, verbose=True)
mag.setConstraintWeights(wz)

if sep_iter > 0:
    print('Running initial separate inversions. ({} iterations.)'.format(sep_iter))
print("Doing initial DC inversion")
mgr.invert(verbose=verbose, lam=lambda_dc, max_iter=sep_iter)

print("Doing initial Magnetic inversion")
mag.run(data_mag, verbose=verbose, lam=lambda_mag, max_iter=sep_iter, absoluteError=abserror_mag)

print('Running structurally coupled inversions. ({} iterations)'.format(max_iter - sep_iter))

statistics = dict(mag=dict(), dc=dict())
statistics['dc']['chi2'] = []
statistics['mag']['chi2'] = []
statistics['dc']['rrms'] = []
statistics['mag']['rrms'] = []
statistics['dc']['rms'] = []
statistics['mag']['rms'] = []
statistics['dc']['iter'] = []
statistics['mag']['iter'] = []

chi2_thresh = 1.2
stop = False  # main stop criterion
while not stop:
    statistics['dc']['chi2'].append(mgr.inv.chi2())
    statistics['mag']['chi2'].append(mag.chi2())
    statistics['dc']['rrms'].append(mgr.inv.relrms())
    statistics['mag']['rrms'].append(mag.relrms())
    statistics['dc']['rms'].append(mgr.inv.absrms())
    statistics['mag']['rms'].append(mag.absrms())
    statistics['dc']['iter'].append(mgr.inv.iter)
    statistics['mag']['iter'].append(mag.iter)

    do_mag_step = mag.chi2() > chi2_thresh
    do_dc_step = mgr.inv.chi2() > chi2_thresh

    mag.fop.region(2).fillConstraintsWeightWithFlatWeight() # not found
    mgr.inv.region(2).fillConstraintsWeightWithFlatWeight() # not found

    roughness_dc = mgr.inv.roughness() # not found
    roughness_mag = mag.roughness() # not found
    cweight_mag = np.minimum((a / (np.absolute(roughness_dc) + a) + b)**c, 1.0)
    cweight_dc = np.minimum((a / (np.absolute(roughness_mag) + a) + b)**c, 1.0)

    if do_dc_step:
        if do_mag_step:
            mgr.inv.setConstraintWeights(cweight_dc)
        print("Doing DC inversion")
        mgr.inv.oneStep() # not found
        #if not mgr.inv.verbose():
        #    mgr.inv.echoStatus()

    if do_mag_step:
        if do_dc_step:
            mag.setConstraintWeights(cweight_mag)
        print("Doing Magnetic inversion")
        mag.oneStep() # not found
        #if not mag.verbose():
        #    mag.echoStatus()

    mag_iter = mag.iter >= max_iter
    dc_iter = mgr.inv.iter >= max_iter

    stop = (not (do_mag_step or do_dc_step)) or mag_iter or dc_iter

resistivity = mgr.inv.model()
susceptibility = mag.model()

Expected behavior

After the initial DC and Magnetic inversion (so far the code worked) it should go in the while loop, calculate cweight with the model roughness and if chi2 of the dc- or magnetic-inversion is above chi2_thres do another inversion step with updated cweights until it reaches one of the abort criterias.

Actual behavior

It seems some of the functions are no longer existing/outdated:

1) .region(2).fillConstraintsWeightWithFlatWeight() I not fully understand what this does, so I haven't found a proper replacement yet.

2) .roughness() it seems to exist in core/src/..., but I can't access the roughness from InversionManager

3) .oneStep() only worked in the LSQRInversion() class but not in the general pg.Inversion or the ERTManager

[halbmy]

Yes, the codes are indeed a bit outdated as they are based on pygimli<=1.0 where the inversion was purely a core (C++) object.

From v1.1 on, the inversion class is purely Python but it still uses a core Inversion under the hood for certain tasks (unless these are replaced) that can be accessed by inv.inv. Whereas 1. starts with a zWeight-defined iteration, 2. and 3. should be available. We should indeed provide an SCCI framework and example. I refer to the newer paper of Skibbe et al. (2021) and the COMET code (https://gitlab.com/Skibbe/comet under snmr/scci) where pyGIMLi 1.1.x was used.

References Skibbe, N., Günther, T. & Müller-Petke, M. (2021): Improved hydrogeophysical imaging by structural coupling of two-dimensional magnetic resonance and electrical resistivity tomography. Geophysics 86 (5), WB135-WB146, https://doi.org/10.1190/geo2020-0593.1

[pi1379]

Thanks Thomas, I applied your corrections, made a few own adjustments to my code and it worked. The Skibbe paper and COMET looks interesting. As soon as I have more time I will try to go deeper in the COMET code.

One more question: In the Ronczka paper you take the roughness vector with Wm for the cWeight calculation and in the Skibbe paper you only take the pureRoughness without Wm. Is there any advantage to apply the latter?

[halbmy]

Great it works! I am looking forward to seeing some results out there and hope we can have an example out there.

In my understanding it is essential to use the actual gradients in the model, i.e. pureRoughness(), for computing the coupling and not the already weighted roughness() where the coupling of the last iteration is already in. But I guess there might be further modifications to improve the matter once we have a flexible framework. Would be good to have @NicoSkibbe involved.

[halbmy]

A few more comments:

1. The flat weight function should not be important and will be added to the pure Python inversion.
2. inv.roughness() exists again from v1.5.2, supporting both weigthed and pure roughness
3. inv.oneStep() should work for all inversions

Independent on this, we are redesigning the Inversion classes to pure Python towards pygimli v1.6. Then is a good point to start over again.

[pi1379]

Hi,

Recently I updated to v1.5.2 and tried your updates to my code (remove the second .inv again), but it doesn't seem to work.

This error also appears for the inv.oneStep() for both the ERT and general Inversion. However the "old" version with the second .inv works as appropriate in v1.5.2, which is fine for me.

My recent script:

import pygimli as pg
import numpy as np
import pygimli.meshtools as mt
from pygimli.physics import ert
from pygimli.physics.gravimetry import MagneticsModelling
#from pygimli.physics.gravimetry import MagManager

#=LOAD DATA====================================================================
# DC
data_dc = pg.getExampleData("ert/gallery3d.dat")
data_dc["k"] = ert.geometricFactors(data_dc, dim=3)
data_dc["err"] = ert.estimateError(data_dc, relativeError=0.02)
print(data_dc)

# Magnetic
coord_mag = np.loadtxt('mag_synth.txt', delimiter='\t', usecols=(0, 1, 2), skiprows=1, dtype='float', unpack=True, encoding='latin1')
points = coord_mag.T

data_mag = np.loadtxt('mag_synth.txt', delimiter='\t', usecols=3, skiprows=1, dtype='float', unpack=True, encoding='latin1')
abserror_mag = 0.2 # nT

#=3D MESH GENERATION===========================================================
# create 3d mesh
plc = mt.createParaMeshPLC3D(data_dc, paraDepth=12, paraMaxCellSize=3,
                             surfaceMeshQuality=34)
mesh = mt.createMesh(plc, quality=1.3)
print(mesh)
#_ = pg.show(mesh, markers=True, showMesh=True, color='white')

# only use inner mesh for magnetic inversion, otherwise the jacobi matrix has the wrong size
#mesh_mag = mesh
mesh_mag = pg.Mesh()
mesh_mag.createMeshByMarker(mesh,2)
#for c in mesh_mag.cells():
#    c.setMarker(1)
mesh_mag.createNeighborInfos()
print(mesh_mag)
#_ = pg.show(mesh_mag, markers=True, showMesh=True, color='white')

#=MAGNETIC IGRF, DEPTH WEIGHTING===============================================
# calculate igrf, set component
F, I, D = 50000, 75, 25  # total field in nT
H = F * np.cos(np.deg2rad(I))
X = H * np.cos(np.deg2rad(D))
Y = H * np.sin(np.deg2rad(D))
Z = F * np.sin(np.deg2rad(I))
igrf = [D, I, H, X, Y, Z, F]
cmp = ["TFA"]

# depth weighting according to Li, Y. & Oldenburg, D. (1996): 3-D inversion of magnetic data. Geophysics 61(2), 394-40
bz = np.array([abs(b.center().z()) for b in mesh_mag.boundaries() if not b.outside()])
z0 = 25
wz = 100 / (bz+z0)**1.5
print(min(wz), max(wz))

#=INVERSION====================================================================
# variables
lambda_dc=30.
lambda_mag=300.
a=0.1
b=0.1
c=1.0
max_iter=15
sep_iter=2
sep_iter_mag=4
max_iter_dc=15
verbose=True

# start inversion class
mgr = ert.ERTManager(data_dc)
fop_mag = MagneticsModelling(mesh=mesh_mag, points=points, cmp=cmp, igrf=igrf)
mag = pg.Inversion(fop=fop_mag, verbose=True)
#mag = MagManager(data='mag_synth.txt')
#mag.fop.setRegionProperties(regionNr=1, background=True)
#mag.fop.setRegionProperties(regionNr=2, startModel=1e-2, limits=[0, 0.8])
#mag.setRegularization(limits=[0, 0.8])
mag.setConstraintWeights(wz)

if sep_iter > 0:
    print('Running initial separate inversions. ({} iterations.)'.format(sep_iter))
print("Doing initial DC inversion")
mgr.invert(mesh=mesh, verbose=verbose, lam=lambda_dc, maxIter=sep_iter)

print("Doing initial Magnetic inversion")
mag.run(data_mag, verbose=verbose, lam=lambda_mag, maxIter=sep_iter_mag, absoluteError=abserror_mag)#, zWeight=wz)

print('Running structurally coupled inversions. ({} iterations)'.format(max_iter - sep_iter))

statistics = dict(mag=dict(), dc=dict())
statistics['dc']['chi2'] = []
statistics['mag']['chi2'] = []
statistics['dc']['rrms'] = []
statistics['mag']['rrms'] = []
statistics['dc']['rms'] = []
statistics['mag']['rms'] = []
statistics['dc']['iter'] = []
statistics['mag']['iter'] = []

chi2_thresh = 1.2
stop = False  # main stop criterion
while not stop:
    statistics['dc']['chi2'].append(mgr.inv.chi2())
    statistics['mag']['chi2'].append(mag.chi2())
    statistics['dc']['rrms'].append(mgr.inv.relrms())
    statistics['mag']['rrms'].append(mag.relrms())
    statistics['dc']['rms'].append(mgr.inv.absrms())
    statistics['mag']['rms'].append(mag.absrms())
    statistics['dc']['iter'].append(mgr.inv.iter)
    statistics['mag']['iter'].append(mag.iter)

    # set chi2 threshold
    do_mag_step = mag.chi2() > chi2_thresh
    do_dc_step = mgr.inv.chi2() > chi2_thresh

    # take roughness for cweights
    #roughness_dc = mgr.inv.inv.roughness()
    roughness_dc = mgr.inv.roughness()  # 1.5.2
    #roughness_mag = mag.inv.roughness()
    roughness_mag = mag.roughness()  # 1.5.2

    # or use pure roughness
    #roughness_dc = mgr.inv.pureRoughness(mgr.inv.model)
    #roughness_mag = mag.pureRoughness(mag.model)

    # calculate cweights
    cweight_mag = np.minimum((a / (np.absolute(roughness_dc) + a) + b)**c, 1.0)
    cweight_dc = np.minimum((a / (np.absolute(roughness_mag) + a) + b)**c, 1.0)

    # or use common cweights
    cweight_new = cweight_mag * cweight_dc

    if do_dc_step:
        if do_mag_step:
            mgr.inv.setConstraintWeights(cweight_new)
        print("Doing DC inversion")
        mgr.inv.oneStep()  # 1.5.2
        #mgr.inv.inv.oneStep()

    if do_mag_step:
        if do_dc_step:
            mag.setConstraintWeights(cweight_new)
        print("Doing Magnetic inversion")
        mag.oneStep()  # 1.5.2
        #mag.inv.oneStep()

    mag_iter = mag.iter >= max_iter
    dc_iter = mgr.inv.iter >= max_iter_dc

    stop = (not (do_mag_step or do_dc_step)) or mag_iter or dc_iter

resistivity = mgr.inv.inv.model()
susceptibility = mag.inv.model()

mesh_mag['res'] = resistivity
mesh_mag['sus'] = susceptibility

mesh_mag.exportVTK('coupled.vtk')

[halbmy]

I don't reall get where you still have to use mgr.inv.inv.

• inv.inv.roughness() is replaced by inv.roughness() (you should call it with weighted=False)
• inv.inv.model() is more easily accessed by inv.model (no brackets)
• in the classic inverison, there is no inv.oneStep(), only in LSQRInversion, so you will still have to call mgr.inv.inv.oneStep(), but this currently is being replaced for version v1.6