Error in MarquardtInversion of TDEM

[ahsan435]

Problem description

I am working on the Marquardt Inversion of TDEM code, which is running, but the results are not matching each other.

Your environment

Please provide the output of print(pygimli.Report()) here. If that does not work, please provide some additional information on your:

Date: February 20 23:09:06 2024 China Standard Time

            OS : Windows
        CPU(s) : 32
       Machine : AMD64
  Architecture : 64bit
           RAM : 31.7 GiB
   Environment : Jupyter

Python 3.9.18 | packaged by conda-forge | (main, Dec 23 2023, 16:29:04) [MSC v.1929 64 bit (AMD64)]

       pygimli : 1.4.6
        pgcore : 1.4.0
         numpy : 1.26.4
    matplotlib : 3.8.0
         scipy : 1.11.4
       IPython : 8.15.0
       pyvista : 0.43.3

Intel(R) oneAPI Math Kernel Library Version 2023.2-Product Build 20230613 for Intel(R) 64 architecture applications

Steps to reproduce

import numpy as np
import matplotlib.pyplot as plt
import pygimli as pg
from pygimli.physics.em import VMDTimeDomainModelling
from pygimli.viewer.mpl import drawModel1D
from pygimli.frameworks import MarquardtInversion
from scipy.interpolate import interp1d

# Define the synthetic model parameters
synThk = [300, 100]
synRes = [100, 10, 100]
nlay = len(synRes)  # number of layers
errorabsTEM_data = 1.  # absolute (ppm of primary signal)

fig, (ax1) = plt.subplots(figsize=(2, 5), ncols=1)

drawModel1D(ax1, synThk, synRes, plot='semilogx', color='C0', label="synth")

# MaxMin-10
times = np.logspace(np.log10(1e-3), np.log10(1e-1), 64)

# Initialize VMDTimeDomainModelling
TEM_modeling = VMDTimeDomainModelling(times=times, nLayers=nlay, txArea=625, verbose=True)

# Extract data
TEM_data = TEM_modeling(synThk + synRes)

# Plotting subplots
plt.figure(figsize=(7, 5))
plt.semilogx(times, TEM_data, label='Observed Data', marker='o')
#plt.semilogx(times, invEM.response, label='Inverted Data', marker='x')
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('Apparent Resistivity (Ohm.m)')
plt.show()

# Set up Marquardt inversion
startModel = [100]*nlay

error = 0 * abs(TEM_data)

def run(self, dataVals, errorVals=None, **kwargs):
    if errorVals is None:  # use absoluteError and/or relativeError instead
            absErr = kwargs.pop("absoluteError", 0)
            relErr = kwargs.pop("relativeError",
                                0.01 if np.allclose(absErr, 0) else 0)
            errorVals = pg.abs(absErr / dataVals) + relErr
invTEM = MarquardtInversion(fop=TEM_modeling, verbose=True)
inversion_model = invTEM.run(dataVals=TEM_data, relativeError=error, startModel=startModel)

fig, (ax1) = plt.subplots(figsize=(2, 5), ncols=1)

drawModel1D(ax1, synThk, synRes, plot='semilogx', color='C0', label="synth")
drawModel1D(ax1, model=inversion_model, color='C1', label="EM")

If possible, please add one or more labels to your issue, e.g. if you expect that your issue is rather a question than a problem with the code, please add the label "question".

[halbmy]

MaxMin-10/Promys instrument, 1m above the ground

The MaxMin/Promis instrument is a frequency-domain EM instrument, for which a dedicated operator exists. VMDTimeDomain does not serve for it.

[ahsan435]

Sir, how to modify it for TDEM? Sorry for that I forgot to remove tha 1m from ground, I convert that from joint inversion of DC-EM code so it was written on that but I think so code line is totally base on TDEM.

[halbmy]

1. You can delete the self-created def run function as it is not used.
2. You specify a zero error which does not make sense (and will be overridden), use any low value, e.g. relativeError=0.03
3. For a 3-layer model, the model vector consists of five unknowns (3 resistivities plus 2 thicknesses), your starting vector will be ignored. You should use [100]*(nlay*2-1)

Doing so one can nicely fit the data with chi-square below 1

plt.figure(figsize=(7, 5))
plt.loglog(times, TEM_data, label='Observed Data', marker='o')
plt.semilogx(times, invTEM.response, label='Inverted Data', marker='x')
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('Apparent Resistivity (Ohm.m)')

This means the retrieved model (orange) is equivalent to the synthetic one The resistivity of the upper layer is nicely retrieved as well as the conductance (thickness * conductivity) of the conductor, and of course you cannot resolve the resistor below the conductor signal well, at least not wich such a small loop (25x25m) as the apparent resistivity curve demonstrates.

If you increase the loop to 250x250m (txArea=62500) then you can resolve it.

[ahsan435]

Dear Sir, I am still unable to match results like yours. I made the changes according to your comments but did not get the same results. Can you share your code, please, so I can read it and see how it will work? I am new to coding and learning the PYGIMLI.

import numpy as np
import matplotlib.pyplot as plt
import pygimli as pg
from pygimli.physics.em import VMDTimeDomainModelling
from pygimli.viewer.mpl import drawModel1D
from pygimli.frameworks import MarquardtInversion
from scipy.interpolate import interp1d

# Define the synthetic model parameters
synThk = [300, 100]
synRes = [100, 10, 100]
nlay = len(synRes)  # number of layers
errorabsTEM_data = 1.  # absolute (ppm of primary signal)

fig, (ax1) = plt.subplots(figsize=(7, 5), ncols=1)

drawModel1D(ax1, synThk, synRes, plot='semilogx', color='C0', label="synth")

# MaxMin-10
times = np.logspace(np.log10(1e-3), np.log10(1e-1), 64)

# Initialize VMDTimeDomainModelling
TEM_modeling = VMDTimeDomainModelling(times=times, nLayers=nlay, txArea=62500, verbose=True)

# Extract data
TEM_data = TEM_modeling(synThk + synRes)

Plotting subplots
plt.figure(figsize=(7, 5))
plt.loglog(times, TEM_data, label='Observed Data', marker='o')

# Set up Marquardt inversion
startModel = [100]*(nlay*2-1)  # Start model for 3-layer model with 2 resistivities and 1 thickness

# Specify error
error = np.zeros_like(TEM_data) + 0.03  # Using a low relative error

invTEM = MarquardtInversion(fop=TEM_modeling, verbose=True)
inversion_model = invTEM.run(dataVals=TEM_data, relativeError=error, startModel=startModel)

plt.semilogx(times, invTEM.response, label='Inverted Data', marker='x')
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('Apparent Resistivity (Ohm.m)')
plt.title('TEM Response')
plt.grid(True)
plt.show()

fig, (ax1) = plt.subplots(figsize=(7, 5), ncols=1)

drawModel1D(ax1, synThk, synRes, plot='semilogx', color='C0', label="synth")
drawModel1D(ax1, model=inversion_model, color='C1', label="EM")
plt.title('Synthetic Model vs Inversion Result')
plt.grid(True)
plt.show()

[halbmy]

Just leave out the startmodel in the inversion or increase the first value. 100 Ohmm for the first layer thickness is too small and the inversion stucks. By default, thickness is derived from the time vector. Everything else in your code is correct.

[ahsan435]

Thank you very much, sir. I got the results perfectly now Thanks alot

[ahsan435]

Sir when I am increasing the number of layers and thickness values then both results are not matching in that case what should I do? when I am changing these values # Define the synthetic model parameters synThk = [5, 15, 15] synRes = [1000, 100, 500, 20] nlay = len(synRes) # number of layers errorabsTEM_data = 1. # absolute (ppm of primary signal)

Initialize VMDTimeDomainModelling

times = np.logspace(np.log10(1e-3), np.log10(1e-1), 64)

TEM_modeling = VMDTimeDomainModelling(times=times, nLayers=nlay, txArea=625000, verbose=True)

Generate synthetic data

TEM_data = TEM_modeling(synThk + synRes)

Set up Marquardt inversion

startModel = [1000](nlay2-1) # Start model for 3-layer model with 2 resistivities and 1 thickness

error = np.zeros_like(TEM_data) + 0.03 # Using a low relative error invTEM = MarquardtInversion(fop=TEM_modeling, verbose=True) inversion_model = invTEM.run(dataVals=TEM_data, relativeError=error)

then results are like this only issue in this results

[ahsan435]

Sir Thank you very much for your support and guidance. My issue is resolved. Please close this issue. Thanks

[ahsan435]

Dear Mr. Subenyu, I hope you are doing well, In the VMDTimeDomainModeling, this class directly works on the apparent resistivity When you import this function, it will calculate apparent resistivity from TDEM Pygimli, which is not working on the magnetic field or electric field So if you need apparent resistivity, just import that class directly and if you want to use it without function, then use this in the code

# Rho-App calculation following
Arec = 1
mu0 = 4e-7 * np.pi
rhoa = radius**(4/3) * Arec**(2/3) * mu0**(5/3) / (20**(2/3) * np.pi**(1/3))
rhoa /= times**(5/3) * (em * 2 * np.pi * radius)**(2/3)

Dear adjust according to your code.

Thanks

On Mon, 29 Apr 2024 at 07:55, subenyu @.***> wrote:

Hi ahsan @ahsan435 https://github.com/ahsan435

I have one question that you can get the apparent resisitivity by TEM_data = TEM_modeling(synThk + synRes) and the apparent resisitivity is using the formula ρa=((ATx∗I∗μ0)/(30B))2/3∗4e−7 So do you do you know how can i calculate apparent resistivity by the formula ρa=(ARx∗ATx∗μ0/20/(U/I))2/3∗t−5/3∗4e−7 Thank you very muh.

— Reply to this email directly, view it on GitHub https://github.com/gimli-org/gimli/issues/654#issuecomment-2081704603, or unsubscribe https://github.com/notifications/unsubscribe-auth/A64356MYPAKCCEUWHZX6CDDY7WD6DAVCNFSM6AAAAABDRK5LDGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDAOBRG4YDINRQGM . You are receiving this because you were mentioned.Message ID: @.***>