PF: Gravity: Laguna del Maule Bouguer GravityΒΆ

This notebook illustrates the SimPEG code used to invert Bouguer gravity data collected at Laguna del Maule volcanic field, Chile. Refer to Miller et al 2016 EPSL for full details.

We run the inversion in two steps. Firstly creating a L2 model and then applying an Lp norm to produce a compact model.

Craig Miller

  • ../../../_images/sphx_glr_plot_laguna_del_maule_inversion_001.png
  • ../../../_images/sphx_glr_plot_laguna_del_maule_inversion_002.png
  • ../../../_images/sphx_glr_plot_laguna_del_maule_inversion_003.png

Out:

Downloading https://storage.googleapis.com/simpeg/Chile_GRAV_4_Miller/Chile_GRAV_4_Miller.tar.gz
   saved to: /Users/lindseyjh/git/simpeg/simpeg/examples/04-grav/Chile_GRAV_4_Miller.tar.gz
Download completed!
Begin calculation of distance weighting for R= 3.0
Done 0.0 %
Done 10.0 %
Done 20.0 %
Done 30.0 %
Done 40.0 %
Done 50.0 %
Done 60.0 %
Done 70.0 %
Done 80.0 %
Done 90.0 %
Done 100% ...distance weighting completed!!

SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||
Begin calculation of forward operator: z
Done 0.0 %
Done 10.0 %
Done 20.0 %
Done 30.0 %
Done 40.0 %
Done 50.0 %
Done 60.0 %
Done 70.0 %
Done 80.0 %
Done 90.0 %
Done 100% ...forward operator completed!!


    SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
    ***Done using same Solver and solverOpts as the problem***
model has any nan: 0
=============================== Projected GNCG ===============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment
-----------------------------------------------------------------------------
x0 has any nan: 0
   0  2.14e-03  1.30e+06  1.87e+01  1.30e+06    1.95e+02      0
   1  1.07e-03  3.47e+03  7.18e+06  1.12e+04    1.15e+02      0
   2  5.35e-04  1.71e+03  8.31e+06  6.16e+03    8.50e+01      0   Skip BFGS
   3  2.68e-04  7.96e+02  9.50e+06  3.34e+03    6.92e+01      0   Skip BFGS
   4  1.34e-04  3.39e+02  1.07e+07  1.77e+03    5.66e+01      0   Skip BFGS
   5  6.69e-05  1.32e+02  1.17e+07  9.16e+02    4.12e+01      0   Skip BFGS
Reached starting chifact with l2-norm regularization: Start IRLS steps...
eps_p: 0.5011757218259212 eps_q: 0.5011757218259212
Eps_p: 0.41764643485493436
Eps_q: 0.41764643485493436
delta phim: 1.315e+06
   6  3.34e-05  5.07e+01  2.24e+07  8.01e+02    3.87e+01      0   Skip BFGS
Beta search step
   7  9.00e-05  5.07e+01  2.24e+07  2.07e+03    6.92e+01      0   Skip BFGS
Beta search step
   8  6.26e-05  5.07e+01  2.24e+07  1.45e+03    5.54e+01      0
Eps_p: 0.34803869571244533
Eps_q: 0.34803869571244533
delta phim: 3.051e+00
   9  6.26e-05  9.17e+01  2.55e+07  1.69e+03    4.97e+01      0   Skip BFGS
Eps_p: 0.2900322464270378
Eps_q: 0.2900322464270378
delta phim: 3.035e-01
  10  6.26e-05  1.24e+02  2.99e+07  2.00e+03    1.13e+02      0
Eps_p: 0.24169353868919816
Eps_q: 0.24169353868919816
delta phim: 2.670e-01
  11  6.26e-05  1.38e+02  3.24e+07  2.17e+03    1.63e+02      1
Beta search step
  12  3.92e-05  1.38e+02  3.24e+07  1.41e+03    1.57e+02      1
Eps_p: 0.20141128224099847
Eps_q: 0.20141128224099847
delta phim: 2.664e-01
  13  3.92e-05  7.19e+01  3.67e+07  1.51e+03    1.40e+02      0
Eps_p: 0.16784273520083207
Eps_q: 0.16784273520083207
delta phim: 2.533e-01
  14  3.92e-05  1.22e+02  3.80e+07  1.61e+03    1.66e+02      0
Eps_p: 0.13986894600069338
Eps_q: 0.13986894600069338
delta phim: 2.174e-01
  15  3.92e-05  1.13e+02  4.04e+07  1.70e+03    1.78e+02      1
Beta search step
  16  2.46e-05  1.13e+02  4.04e+07  1.11e+03    1.76e+02      0
Eps_p: 0.11655745500057782
Eps_q: 0.11655745500057782
delta phim: 2.000e-01
  17  2.46e-05  5.34e+01  4.45e+07  1.15e+03    7.29e+01      0
Eps_p: 0.09713121250048153
Eps_q: 0.09713121250048153
delta phim: 1.847e-01
  18  2.46e-05  9.97e+01  4.46e+07  1.20e+03    6.74e+01      0
Eps_p: 0.08094267708373461
Eps_q: 0.08094267708373461
delta phim: 1.488e-01
  19  2.46e-05  7.68e+01  4.66e+07  1.22e+03    9.14e+01      0
Eps_p: 0.06745223090311217
Eps_q: 0.06745223090311217
delta phim: 1.053e-01
  20  2.46e-05  1.20e+02  4.59e+07  1.25e+03    7.66e+01      0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 2.4581e+01 <= tolF*(1+|f0|) = 1.3018e+05
0 : |xc-x_last| = 1.6646e+00 <= tolX*(1+|x0|) = 1.0400e-01
0 : |proj(x-g)-x|    = 7.6594e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 7.6594e+01 <= 1e3*eps       = 1.0000e-02
1 : maxIter   =      20    <= iter          =     20
------------------------- DONE! -------------------------

import os
import shutil
import tarfile
import SimPEG.PF as PF
from SimPEG import Maps, Regularization, Optimization, DataMisfit,\
                   InvProblem, Directives, Inversion
from SimPEG.Utils.io_utils import download
import matplotlib.pyplot as plt
import numpy as np


def run(plotIt=True, cleanAfterRun=True):

    # Start by downloading files from the remote repository
    # directory where the downloaded files are

    url = "https://storage.googleapis.com/simpeg/Chile_GRAV_4_Miller/Chile_GRAV_4_Miller.tar.gz"
    downloads = download(url, overwrite=True)
    basePath = downloads.split(".")[0]

    # unzip the tarfile
    tar = tarfile.open(downloads, "r")
    tar.extractall()
    tar.close()

    input_file = basePath + os.path.sep + 'LdM_input_file.inp'
    # %% User input
    # Plotting parameters, max and min densities in g/cc
    vmin = -0.6
    vmax = 0.6

    # weight exponent for default weighting
    wgtexp = 3.
    # %%
    # Read in the input file which included all parameters at once
    # (mesh, topo, model, survey, inv param, etc.)
    driver = PF.GravityDriver.GravityDriver_Inv(input_file)
    # %%
    # Now we need to create the survey and model information.

    # Access the mesh and survey information
    mesh = driver.mesh
    survey = driver.survey

    # define gravity survey locations
    rxLoc = survey.srcField.rxList[0].locs

    # define gravity data and errors
    d = survey.dobs
    wd = survey.std

    # Get the active cells
    active = driver.activeCells
    nC = len(active)  # Number of active cells

    # Create active map to go from reduce set to full
    activeMap = Maps.InjectActiveCells(mesh, active, -100)

    # Create static map
    static = driver.staticCells
    dynamic = driver.dynamicCells

    staticCells = Maps.InjectActiveCells(
        None, dynamic, driver.m0[static], nC=nC
    )
    mstart = driver.m0[dynamic]

    # Get index of the center
    midx = int(mesh.nCx/2)
    # %%
    # Now that we have a model and a survey we can build the linear system ...
    # Create the forward model operator
    prob = PF.Gravity.GravityIntegral(mesh, rhoMap=staticCells,
                                      actInd=active)
    prob.solverOpts['accuracyTol'] = 1e-4

    # Pair the survey and problem
    survey.pair(prob)

    # Apply depth weighting
    wr = PF.Magnetics.get_dist_wgt(mesh, rxLoc, active, wgtexp,
                                   np.min(mesh.hx)/4.)
    wr = wr**2.

    # %% Create inversion objects
    reg = Regularization.Sparse(mesh, indActive=active,
                                mapping=staticCells, gradientType='total')
    reg.mref = driver.mref[dynamic]
    reg.cell_weights = wr * mesh.vol[active]
    reg.norms = np.c_[0., 1., 1., 1.]
    # reg.norms = driver.lpnorms

    # Specify how the optimization will proceed
    opt = Optimization.ProjectedGNCG(maxIter=20, lower=driver.bounds[0],
                                     upper=driver.bounds[1], maxIterLS=10,
                                     maxIterCG=20, tolCG=1e-3)

    # Define misfit function (obs-calc)
    dmis = DataMisfit.l2_DataMisfit(survey)
    dmis.W = 1./wd

    # create the default L2 inverse problem from the above objects
    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)

    # Specify how the initial beta is found
    betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e-2)

    # IRLS sets up the Lp inversion problem
    # Set the eps parameter parameter in Line 11 of the
    # input file based on the distribution of model (DEFAULT = 95th %ile)
    IRLS = Directives.Update_IRLS(f_min_change=1e-4, maxIRLSiter=40, beta_tol=5e-1)

    # Preconditioning refreshing for each IRLS iteration
    update_Jacobi = Directives.UpdatePreconditioner()

    # Create combined the L2 and Lp problem
    inv = Inversion.BaseInversion(invProb,
                                  directiveList=[IRLS, update_Jacobi, betaest])

    # %%
    # Run L2 and Lp inversion
    mrec = inv.run(mstart)

    if cleanAfterRun:
        os.remove(downloads)
        shutil.rmtree(basePath)

    # %%
    if plotIt:
        # Plot observed data
        PF.Magnetics.plot_obs_2D(rxLoc, d, 'Observed Data')

        # %%
        # Write output model and data files and print misft stats.

        # reconstructing l2 model mesh with air cells and active dynamic cells
        L2out = activeMap * invProb.l2model

        # reconstructing lp model mesh with air cells and active dynamic cells
        Lpout = activeMap*mrec

        # %%
        # Plot out sections and histograms of the smooth l2 model.
        # The ind= parameter is the slice of the model from top down.
        yslice = midx + 1
        L2out[L2out == -100] = np.nan  # set "air" to nan

        plt.figure(figsize=(10, 7))
        plt.suptitle('Smooth Inversion: Depth weight = ' + str(wgtexp))
        ax = plt.subplot(221)
        dat1 = mesh.plotSlice(L2out, ax=ax, normal='Z', ind=-16,
                              clim=(vmin, vmax), pcolorOpts={'cmap': 'bwr'})
        plt.plot(np.array([mesh.vectorCCx[0], mesh.vectorCCx[-1]]),
                 np.array([mesh.vectorCCy[yslice], mesh.vectorCCy[yslice]]),
                 c='gray', linestyle='--')
        plt.scatter(rxLoc[0:, 0], rxLoc[0:, 1], color='k', s=1)
        plt.title('Z: ' + str(mesh.vectorCCz[-16]) + ' m')
        plt.xlabel('Easting (m)')
        plt.ylabel('Northing (m)')
        plt.gca().set_aspect('equal', adjustable='box')
        cb = plt.colorbar(dat1[0], orientation="vertical",
                          ticks=np.linspace(vmin, vmax, 4))
        cb.set_label('Density (g/cc$^3$)')

        ax = plt.subplot(222)
        dat = mesh.plotSlice(L2out, ax=ax, normal='Z', ind=-27,
                             clim=(vmin, vmax), pcolorOpts={'cmap': 'bwr'})
        plt.plot(np.array([mesh.vectorCCx[0], mesh.vectorCCx[-1]]),
                 np.array([mesh.vectorCCy[yslice], mesh.vectorCCy[yslice]]),
                 c='gray', linestyle='--')
        plt.scatter(rxLoc[0:, 0], rxLoc[0:, 1], color='k', s=1)
        plt.title('Z: ' + str(mesh.vectorCCz[-27]) + ' m')
        plt.xlabel('Easting (m)')
        plt.ylabel('Northing (m)')
        plt.gca().set_aspect('equal', adjustable='box')
        cb = plt.colorbar(dat1[0], orientation="vertical",
                          ticks=np.linspace(vmin, vmax, 4))
        cb.set_label('Density (g/cc$^3$)')

        ax = plt.subplot(212)
        mesh.plotSlice(L2out, ax=ax, normal='Y', ind=yslice,
                       clim=(vmin, vmax), pcolorOpts={'cmap': 'bwr'})
        plt.title('Cross Section')
        plt.xlabel('Easting(m)')
        plt.ylabel('Elevation')
        plt.gca().set_aspect('equal', adjustable='box')
        cb = plt.colorbar(dat1[0], orientation="vertical",
                          ticks=np.linspace(vmin, vmax, 4), cmap='bwr')
        cb.set_label('Density (g/cc$^3$)')

        # %%
        # Make plots of Lp model
        yslice = midx + 1
        Lpout[Lpout == -100] = np.nan  # set "air" to nan

        plt.figure(figsize=(10, 7))
        plt.suptitle('Compact Inversion: Depth weight = ' + str(wgtexp) +
                     ': $\epsilon_p$ = ' + str(round(reg.eps_p, 1)) +
                     ': $\epsilon_q$ = ' + str(round(reg.eps_q, 2)))
        ax = plt.subplot(221)
        dat = mesh.plotSlice(Lpout, ax=ax, normal='Z', ind=-16,
                             clim=(vmin, vmax), pcolorOpts={'cmap': 'bwr'})
        plt.plot(np.array([mesh.vectorCCx[0], mesh.vectorCCx[-1]]),
                 np.array([mesh.vectorCCy[yslice], mesh.vectorCCy[yslice]]),
                 c='gray', linestyle='--')
        plt.scatter(rxLoc[0:, 0], rxLoc[0:, 1], color='k', s=1)
        plt.title('Z: ' + str(mesh.vectorCCz[-16]) + ' m')
        plt.xlabel('Easting (m)')
        plt.ylabel('Northing (m)')
        plt.gca().set_aspect('equal', adjustable='box')
        cb = plt.colorbar(dat[0], orientation="vertical",
                          ticks=np.linspace(vmin, vmax, 4))
        cb.set_label('Density (g/cc$^3$)')

        ax = plt.subplot(222)
        dat = mesh.plotSlice(Lpout, ax=ax, normal='Z', ind=-27,
                             clim=(vmin, vmax), pcolorOpts={'cmap': 'bwr'})
        plt.plot(np.array([mesh.vectorCCx[0], mesh.vectorCCx[-1]]),
                 np.array([mesh.vectorCCy[yslice], mesh.vectorCCy[yslice]]),
                 c='gray', linestyle='--')
        plt.scatter(rxLoc[0:, 0], rxLoc[0:, 1], color='k', s=1)
        plt.title('Z: ' + str(mesh.vectorCCz[-27]) + ' m')
        plt.xlabel('Easting (m)')
        plt.ylabel('Northing (m)')
        plt.gca().set_aspect('equal', adjustable='box')
        cb = plt.colorbar(dat[0], orientation="vertical",
                          ticks=np.linspace(vmin, vmax, 4))
        cb.set_label('Density (g/cc$^3$)')

        ax = plt.subplot(212)
        dat = mesh.plotSlice(Lpout, ax=ax, normal='Y', ind=yslice,
                             clim=(vmin, vmax), pcolorOpts={'cmap': 'bwr'})
        plt.title('Cross Section')
        plt.xlabel('Easting (m)')
        plt.ylabel('Elevation (m)')
        plt.gca().set_aspect('equal', adjustable='box')
        cb = plt.colorbar(dat[0], orientation="vertical",
                          ticks=np.linspace(vmin, vmax, 4))
        cb.set_label('Density (g/cc$^3$)')

if __name__ == '__main__':
    run()
    plt.show()

Total running time of the script: ( 4 minutes 58.927 seconds)

Gallery generated by Sphinx-Gallery