Seismic Straight Ray TomographyΒΆ

Example of a straight ray tomography inverse problem.

  • ../../../_images/sphx_glr_plot_seis_straight_ray_tomo_001.png
  • ../../../_images/sphx_glr_plot_seis_straight_ray_tomo_002.png
  • ../../../_images/sphx_glr_plot_seis_straight_ray_tomo_003.png

Out:

('length:', 1.1313708498984762)
('length:', 0.7071067811865476)
Using a seed of:  598
SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||
SimPEG.InvProblem will set Regularization.mref to m0.

    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
============================ Inexact Gauss Newton ============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment
-----------------------------------------------------------------------------
x0 has any nan: 0
   0  5.39e+03  4.55e+02  0.00e+00  4.55e+02    2.42e+02      0
   1  5.39e+03  9.11e+01  3.96e-03  1.12e+02    4.01e+00      0
   2  5.39e+03  9.10e+01  3.96e-03  1.12e+02    2.53e-01      0
   3  6.74e+02  9.10e+01  3.96e-03  9.37e+01    7.39e+01      0   Skip BFGS
   4  6.74e+02  3.83e+01  3.20e-02  5.98e+01    3.71e+01      0
   5  6.74e+02  3.82e+01  3.20e-02  5.98e+01    3.70e+01      0   Skip BFGS
   6  8.43e+01  3.82e+01  3.20e-02  4.09e+01    5.62e+01      0
   7  8.43e+01  3.30e+01  5.75e-02  3.79e+01    5.61e+01      0
   8  8.43e+01  3.09e+01  5.67e-02  3.57e+01    5.60e+01      0
   9  1.05e+01  3.06e+01  5.20e-02  3.12e+01    5.75e+01      0
  10  1.05e+01  2.91e+01  7.06e-02  2.98e+01    5.86e+01      0
  11  1.05e+01  2.89e+01  7.08e-02  2.97e+01    5.89e+01      0
  12  1.32e+00  2.88e+01  7.06e-02  2.89e+01    5.88e+01      0
  13  1.32e+00  2.84e+01  7.80e-02  2.85e+01    5.89e+01      0
  14  1.32e+00  2.73e+01  9.44e-02  2.74e+01    5.91e+01      0   Skip BFGS
  15  1.65e-01  2.72e+01  9.40e-02  2.72e+01    5.92e+01      0   Skip BFGS
  16  1.65e-01  2.72e+01  9.30e-02  2.72e+01    5.92e+01      0
  17  1.65e-01  2.72e+00  1.57e-01  2.75e+00    1.44e+01      0
  18  2.06e-02  1.03e+00  1.72e-01  1.04e+00    7.32e+00      0   Skip BFGS
  19  2.06e-02  7.35e-01  2.11e-01  7.40e-01    5.23e+00      0
  20  2.06e-02  5.95e-01  1.85e-01  5.99e-01    4.95e+00      0   Skip BFGS
  21  2.57e-03  5.76e-01  1.83e-01  5.76e-01    4.81e+00      0
  22  2.57e-03  5.73e-01  1.84e-01  5.73e-01    4.72e+00      0
  23  2.57e-03  3.18e-01  1.68e-01  3.18e-01    4.17e+00      0   Skip BFGS
  24  3.22e-04  2.30e-01  1.69e-01  2.30e-01    3.28e+00      0
  25  3.22e-04  1.97e-01  1.64e-01  1.97e-01    3.12e+00      0
  26  3.22e-04  1.55e-01  1.60e-01  1.55e-01    2.75e+00      0   Skip BFGS
  27  4.02e-05  1.35e-01  1.61e-01  1.35e-01    2.56e+00      0
  28  4.02e-05  1.22e-01  1.59e-01  1.22e-01    2.26e+00      0   Skip BFGS
  29  4.02e-05  1.12e-01  1.60e-01  1.12e-01    2.27e+00      0
  30  5.02e-06  8.45e-02  1.62e-01  8.45e-02    1.89e+00      0   Skip BFGS
  31  5.02e-06  6.56e-02  1.59e-01  6.56e-02    1.79e+00      0
  32  5.02e-06  4.05e-02  1.60e-01  4.05e-02    1.39e+00      0   Skip BFGS
  33  6.28e-07  3.37e-02  1.59e-01  3.37e-02    1.39e+00      0
  34  6.28e-07  3.06e-02  1.58e-01  3.06e-02    1.39e+00      0   Skip BFGS
  35  6.28e-07  2.93e-02  1.58e-01  2.93e-02    1.34e+00      0
  36  7.85e-08  2.87e-02  1.58e-01  2.87e-02    1.33e+00      0   Skip BFGS
  37  7.85e-08  2.82e-02  1.58e-01  2.82e-02    1.32e+00      0
  38  7.85e-08  2.53e-02  1.60e-01  2.53e-02    1.11e+00      0   Skip BFGS
  39  9.81e-09  2.30e-02  1.58e-01  2.30e-02    1.13e+00      0
  40  9.81e-09  1.89e-02  1.55e-01  1.89e-02    1.05e+00      0   Skip BFGS
------------------------- STOP! -------------------------
1 : |fc-fOld| = 4.1382e-03 <= tolF*(1+|f0|) = 4.5647e+01
1 : |xc-x_last| = 4.2267e-02 <= tolX*(1+|x0|) = 3.1000e+00
0 : |proj(x-g)-x|    = 1.0499e+00 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 1.0499e+00 <= 1e3*eps       = 1.0000e-02
1 : maxIter   =      40    <= iter          =     40
------------------------- DONE! -------------------------

from SimPEG import Utils
from SimPEG import Mesh
from SimPEG import Maps
from SimPEG import Survey
from SimPEG import Regularization
from SimPEG import DataMisfit
from SimPEG import Optimization
from SimPEG import InvProblem
from SimPEG import Directives
from SimPEG import Inversion

from SimPEG.SEIS import StraightRay

import numpy as np
import scipy.sparse as sp
import matplotlib.pyplot as plt


def run(plotIt=False):

    O = np.r_[-1.2, -1.]
    D = np.r_[10., 10.]
    x = np.r_[0., 1.]
    y = np.r_[0., 1.]
    print('length:', StraightRay.lengthInCell(O, D, x, y, plotIt=plotIt))
    O = np.r_[0, -1.]
    D = np.r_[1., 1.]*1.5
    print('length:', StraightRay.lengthInCell(O, D, x*2, y*2, plotIt=plotIt))

    nC = 20
    M = Mesh.TensorMesh([nC, nC])
    y = np.linspace(0., 1., nC/2)
    rlocs = np.c_[y*0+M.vectorCCx[-1], y]
    rx = StraightRay.Rx(rlocs, None)

    srcList = [
        StraightRay.Src(loc=np.r_[M.vectorCCx[0], yi], rxList=[rx])
        for yi in y
    ]

    survey = StraightRay.Survey(srcList)
    problem = StraightRay.Problem(M, slownessMap=Maps.IdentityMap(M))
    problem.pair(survey)

    s = Utils.mkvc(Utils.ModelBuilder.randomModel(M.vnC)) + 1.
    survey.dobs = survey.dpred(s)
    survey.std = 0.01

    # Create an optimization program

    reg = Regularization.Tikhonov(M)
    dmis = DataMisfit.l2_DataMisfit(survey)
    opt = Optimization.InexactGaussNewton(maxIter=40)
    opt.remember('xc')
    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
    beta = Directives.BetaSchedule()
    betaest = Directives.BetaEstimate_ByEig()
    inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest])

    # Start the inversion with a model of zeros, and run the inversion
    m0 = np.ones(M.nC)*1.5
    mopt = inv.run(m0)

    if plotIt is True:
        fig, ax = plt.subplots(1, 2, figsize=(8, 4))
        ax[1].plot(survey.dobs)
        ax[1].plot(survey.dpred(m0), 's')
        ax[1].plot(survey.dpred(mopt), 'o')
        ax[1].legend(['dobs', 'starting dpred', 'dpred'])
        M.plotImage(s, ax=ax[0])
        survey.plot(ax=ax[0])
        ax[0].set_title('survey')

        plt.tight_layout()

    if plotIt is True:
        fig, ax = plt.subplots(1, 3, figsize=(12, 4))
        plt.colorbar(M.plotImage(m0, ax=ax[0])[0], ax=ax[0])
        plt.colorbar(M.plotImage(mopt, ax=ax[1])[0], ax=ax[1])
        plt.colorbar(M.plotImage(s, ax=ax[2])[0], ax=ax[2])

        ax[0].set_title('Starting Model')
        ax[1].set_title('Recovered Model')
        ax[2].set_title('True Model')

        plt.tight_layout()

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

Total running time of the script: ( 0 minutes 12.070 seconds)

Gallery generated by Sphinx-Gallery