Issue in the TDEM LCI Inversion

Problem description

I am trying to make the LCI code of TDEM run properly but the results are not good I tried everything but am not getting the best results. Please help me regarding this.

Please describe your issue here.

Your environment

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

Date: Wed Apr 03 22:46:02 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.5.0.post1
        pgcore : 1.5.0
         numpy : 1.26.4
    matplotlib : 3.8.0
         scipy : 1.12.0
       IPython : 8.15.0
        meshio : 5.3.5
        tetgen : 0.6.4
       pyvista : 0.43.4

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

import matplotlib.pyplot as plt
import numpy as np
import pygimli as pg
from pygimli.viewer.mpl import drawModel1D, showStitchedModels
from pygimli.physics.em import VMDTimeDomainModelling

class TDEM2dFOP(pg.core.ModellingBase):
    """TDEM 2d-LCI modelling class based on BlockMatrices."""
    def __init__(self, x, times, txArea, nlay=2, verbose=False):
        super(TDEM2dFOP, self).__init__(verbose)
        self.nlay = nlay
        self.header = {}
        self.pos = None
        self.z = None
        self.topo = None
        self.nx = len(x)
        self.nt = len(times)
        npar = 2 * nlay - 1  # number of parameters per sounding
        self.mesh1d = pg.meshtools.createMesh1D(self.nx, npar)
        self.mesh_ = pg.meshtools.createMesh1D(self.nx, npar)
        self.setMesh(self.mesh_)

        self.fop = pg.physics.em.VMDTimeDomainModelling(times=times, nLayers=nlay, txArea=txArea)

        self.jacobian = pg.matrix.BlockMatrix()
        self.fop_list = []
        for _ in range(self.nx):
            self.fop_list.append(pg.physics.em.VMDTimeDomainModelling(times=times, nLayers=nlay, txArea=txArea))
            n = self.jacobian.addMatrix(self.fop_list[-1].jacobian())
            self.jacobian.addMatrixEntry(n, self.nt * _, npar * _)

        self.jacobian.recalcMatrixSize()
        print(self.jacobian.rows(), self.jacobian.cols())
        self.setJacobian(self.jacobian)

    def response(self, model):
        """Combine forward responses of all soundings."""
        model_array = np.asarray(model).reshape((self.nx, self.nlay * 2 - 1))
        response = pg.Vector(0)
        for model_i in model_array:
            response = pg.cat(response, self.fop.response(model_i))
        return response

    def createJacobian(self, model):
        """Create Jacobian matrix by creating individual Jacobians."""
        model_array = np.asarray(model).reshape((self.nx, self.nlay * 2 - 1))
        for i, model_i in enumerate(model_array):
            self.fop_list[i].createJacobian(model_i)

if __name__ == "__main__":
    # Define the model parameters
    num_positions = 20  # number of sampling positions
    layer_resistivities = [1000, 100, 500, 20]
    num_layers = len(layer_resistivities)  # number of layers
    resistivities = np.ones((num_positions, num_layers)) * layer_resistivities
    x = np.linspace(0, 20, num_positions)[:, None]
    nx = len(x)
    thickness_1 = 5 + 0.2 * np.sin(x * np.pi * 2)
    thickness_2 = 15 + np.sin(x * np.pi * 2)
    thickness_3 = 10 + 0.3 * np.sin(x * np.pi * 2)

    # Time-domain parameters
    num_time_steps = 64
    times = np.logspace(np.log10(1e-4), np.log10(1e-1), num_time_steps)
    transmitter_area = 62500

    true_model = []
    for i in range(num_positions):
        model_i = np.array([thickness_1[i][0], thickness_2[i][0], thickness_3[i][0]] + layer_resistivities)  # True model
        print(model_i)
        true_model.append(model_i)

    try:
        tdem2d_fop = TDEM2dFOP(x, times, transmitter_area, nlay=num_layers)
        data_tem = tdem2d_fop.response(true_model)
        error_tem_abs = 1.  # absolute error (ppm of primary signal)
        error_tem = np.ones_like(data_tem) * error_tem_abs

        # 2D LCI inversion
        # Transformations
        trans_data = pg.trans.Trans()
        trans_thickness = pg.trans.TransLogLU(0.1, 50.)
        trans_resistivity = pg.trans.TransLogLU(10., 2000.)

        for i in range(num_layers - 1):
            tdem2d_fop.region(i).setTransModel(trans_thickness)

        for i in range(num_layers - 1, num_layers * 2 - 1):
            tdem2d_fop.region(i).setTransModel(trans_resistivity)

        # Set constraints
        tdem2d_fop.region(0).setConstraintType(10)
        tdem2d_fop.region(1).setConstraintType(10)

        # Generate starting model by repetition
        starting_model = pg.Vector([10] * (num_layers - 1) + [100] * num_layers)
        model = np.repeat(starting_model, len(x))

        inversion = pg.core.Inversion(data_tem, tdem2d_fop, trans_data)
        inversion.setAbsoluteError(error_tem)
        inversion.setLambda(500)
        inversion.setModel(model)
        inversion.setReferenceModel(model)

        estimated_model = inversion.run()
        estimated_model_array = np.asarray(estimated_model).reshape((nx, num_layers * 2 - 1))
        estimated_model_list = []  # empty array to store the estimated model
        for model_i in estimated_model_array:
            print(model_i)
            estimated_model_list.append(model_i)

        # Plotting the results
        fig, ax = pg.plt.subplots(figsize=(12, 6), ncols=2)

        showStitchedModels(true_model, ax=ax[0], x=None, cMap='Spectral_r', logScale=True, title='True model (Ohm.m)')
        ax[0].set_ylim(-50, 0)
        ax[0].set_xlabel('Distance (m)')
        ax[0].set_ylabel('Depth (m)')
        ax[0].set_title('True model')
        ax[0].legend(loc='lower right', ncol=2, fontsize=10, frameon=False)
        ax[0].grid(True, alpha=0.5)
        ax[0].set_axisbelow(True)

        showStitchedModels(estimated_model_list, ax=ax[1], x=None, cMap='Spectral_r', logScale=True, title='Estimated model (Ohm.m)')
        ax[1].set_ylim(-50, 0)
        ax[1].set_xlabel('Distance (m)')
        ax[1].set_title('Estimated model')
        ax[1].legend(loc='lower right', ncol=2, fontsize=10, frameon=False)
        ax[1].grid(True, alpha=0.5)
        ax[1].set_axisbelow(True)

        plt.show()
    except Exception as e:
        print(f"An error occurred: {e}")

download

gimli-org / gimli

Issue in the TDEM LCI Inversion #693

Problem description

Your environment