

In [4]:

    
from SimPEG import *
import simpegDCIP as DC
from simpegem1d import Utils1D
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib



In [5]:

    
import matplotlib
import matplotlib.pyplot as plt
matplotlib.rcParams.update({'font.size': 20, 'text.usetex': True})
# matplotlib.rcParams.update({'font.size': 16, 'font.family': 'arial'})

1D DC inversion of Schlumberger array

This is an example for 1D DC Sounding inversion. This 1D inversion usually use analytic foward modeling, which is efficient. However, we choose different approach to show flexibility in geophysical inversion through mapping. Here mapping ($M$)indicates transformation of our model to a different space:

$$ \mathbf{m} = M(\mathbf{\sigma}) $$

Now we consider a transformation, which maps 3D conductivity model to 1D layer model. That is, 3D distribution of conducitivity can be parameterized as 1D model. Once we can compute derivative of this transformation, we can change our model space, based on the transformation.

Following example will show you how user can implement this set up with 1D DC inversion example. Note that we have 3D forward modeling mesh.

Step1: Generate mesh



In [6]:

    
cs = 25.
npad = 11
hx = [(cs,npad, -1.3),(cs,41),(cs,npad, 1.3)]
hy = [(cs,npad, -1.3),(cs,17),(cs,npad, 1.3)]
hz = [(cs,npad, -1.3),(cs,20)]



In [7]:

    
mesh = Mesh.TensorMesh([hx, hy, hz], 'CCN')



In [8]:

    
mesh.plotGrid()

Step2: Generating model and mapping (1D to 3D)



In [9]:

    
mapping = Maps.ExpMap(mesh)*Maps.Vertical1DMap(mesh)
siglay1 = 1./(100.)
siglay2 = 1./(500.)
sighalf = 1./(100.)
sigma = np.ones(mesh.nCz)*siglay1
sigma[mesh.vectorCCz<=-100.] = siglay2
sigma[mesh.vectorCCz<-150.] = sighalf
mtrue = np.log(sigma)



In [10]:

    
fig, ax = plt.subplots(1,2, figsize = (18*0.8, 7*0.8))
Utils1D.plotLayer(np.log(sigma), mesh.vectorCCz, 'linear', showlayers=True, ax = ax[0])
ax[0].invert_xaxis()
ax[0].set_ylim(-500, 0)
ax[0].set_xlim(-7, -4)
ax[0].set_xlabel('$log(\sigma)$', fontsize = 25)
ax[0].set_ylabel('Depth (m)', fontsize = 25)
ax[0].text(-7., 10., '(a)', fontsize = 30)
dat = mesh.plotSlice((mapping*mtrue), normal='Y', ind = 9, ax = ax[1])
cb = plt.colorbar(dat[0], ax =ax[1])
ax[1].set_title("Vertical section", fontsize = 25)
cb.set_label("Conductivity (S/m)", fontsize = 25)
ax[1].set_xlabel('Easting (m)', fontsize = 25)
ax[1].set_ylabel(' ', fontsize = 25)
ax[1].set_xlim(-1000., 1000.)
ax[1].set_ylim(-500., 0.)
ax[1].text(-1000., 13., '(b)', fontsize = 30)
fig.savefig('mappingDC.png', dpi=200)



In [11]:

    
# fig, ax = plt.subplots(1,1, figsize = (7, 5))
# dat = mesh.plotSlice((mapping*mtrue), normal='Y', ind = 9, ax = ax)
# cb = plt.colorbar(dat[0], ax =ax)
# ax.set_title("Vertical section", fontsize = 16)
# cb.set_label("Conductivity (S/m)", fontsize = 16)
# ax.set_xlabel('Easting (m)', fontsize = 16)
# ax.set_ylabel('Depth (m)', fontsize = 16)
# ax.set_xlim(-1000., 1000.)
# ax.set_ylim(-500., 0.)
# ax.text(-1000., 20., '(b)', fontsize = 20)
# fig.savefig('cond3d.png', dpi=200)

Step3: Design survey: Schulumberger array

$$ \rho_a = \frac{V}{I}\pi\frac{b(b+a)}{a}$$

Let $b=na$, then we rewrite above equation as:

$$ \rho_a = \frac{V}{I}\pi na(n+1)$$

Since AB/2 can be a good measure for depth of investigation, we express

$$AB/2 = \frac{(2n+1)a}{2}$$



In [12]:

    
ntx = 16



In [13]:

    
xtemp_txP = np.arange(ntx)*(25.)-500.
xtemp_txN = -xtemp_txP
ytemp_tx = np.zeros(ntx)
xtemp_rxP = -50.
xtemp_rxN = 50.
ytemp_rx = 0.
abhalf = abs(xtemp_txP-xtemp_txN)*0.5
a = xtemp_rxN-xtemp_rxP
b = ((xtemp_txN-xtemp_txP)-a)*0.5



In [14]:

    
print a
print b









    



100.0
[ 450.  425.  400.  375.  350.  325.  300.  275.  250.  225.  200.  175.
  150.  125.  100.   75.]



In [15]:

    
fig, ax = plt.subplots(1,1, figsize = (12,3))
for i in range(ntx):
    ax.plot(np.r_[xtemp_txP[i], xtemp_txP[i]], np.r_[0., 0.4-0.01*(i-1)], 'k-', lw = 1)
    ax.plot(np.r_[xtemp_txN[i], xtemp_txN[i]], np.r_[0., 0.4-0.01*(i-1)], 'k-', lw = 1)
    ax.plot(xtemp_txP[i], ytemp_tx[i], 'bo')
    ax.plot(xtemp_txN[i], ytemp_tx[i], 'ro')
    ax.plot(np.r_[xtemp_txP[i], xtemp_txN[i]], np.r_[0.4-0.01*(i-1), 0.4-0.01*(i-1)], 'k-', lw = 1)    

ax.plot(np.r_[xtemp_rxP, xtemp_rxP], np.r_[0., 0.2], 'k-', lw = 1)
ax.plot(np.r_[xtemp_rxN, xtemp_rxN], np.r_[0., 0.2], 'k-', lw = 1)
ax.plot(xtemp_rxP, ytemp_rx, 'ko')
ax.plot(xtemp_rxN, ytemp_rx, 'go')
ax.plot(np.r_[xtemp_rxP, xtemp_rxN], np.r_[0.2, 0.2], 'k-', lw = 1)    

ax.grid(True)    
ax.set_ylim(-0.2,0.6)
# ax.set_xlim(-600,600)









    Out[15]:





(-0.2, 0.6)



In [16]:

    
fig, ax = plt.subplots(1,1, figsize = (6,4))
ax.plot(xtemp_txP, ytemp_tx, 'bo')
ax.plot(xtemp_txN, ytemp_tx, 'ro')
ax.plot(xtemp_rxP, ytemp_rx, 'ko')
ax.plot(xtemp_rxN, ytemp_rx, 'go')
ax.legend(('A (C+)', 'B (C-)', 'M (P+)', 'N (C-)'), fontsize = 14)
mesh.plotSlice(np.log10(mapping*mtrue), grid=True, ax = ax, pcolorOpts={'cmap':'binary'})
ax.set_xlim(-600, 600)
ax.set_ylim(-200, 200)
ax.set_title('Survey geometry (Plan view)')
ax.set_xlabel('Easting (m)')
ax.set_ylabel('Northing (m)')
ax.text(-600, 210, '(a)', fontsize = 16)
fig.savefig('DCsurvey.png', dpi = 200)

We generate tx and rx lists:



In [17]:

    
txlist = []
rx = DC.RxDipole(np.r_[xtemp_rxP, ytemp_rx, -12.5], np.r_[xtemp_rxN, ytemp_rx, -12.5])
for i in range(ntx):    
    tx = DC.SrcDipole([rx], [xtemp_txP[i], ytemp_tx[i], -12.5],[xtemp_txN[i], ytemp_tx[i], -12.5])
    txlist.append(tx)
survey = DC.SurveyDC(txlist)

Step4: Set up problem and pair with survey



In [19]:

    
problem = DC.ProblemDC_CC(mesh, mapping=mapping)
problem.pair(survey)
try:
    from pymatsolver import MumpsSolver
    problem.Solver = MumpsSolver
except Exception, e:
    problem.Solver = SolverLU

Step5: Run survey.dpred to comnpute syntetic data



In [20]:

    
%%time
data = survey.dpred(mtrue)









    



CPU times: user 6.05 s, sys: 708 ms, total: 6.76 s
Wall time: 6.01 s

$$ \rho_a = \frac{V}{I}\pi\frac{b(b+a)}{a}$$

To make synthetic example you can use survey.makeSyntheticData, which generates related setups



In [21]:

    
survey.makeSyntheticData(mtrue,std=0.01,force=True)









    Out[21]:





array([ 0.0186522 ,  0.0211157 ,  0.02374998,  0.02677094,  0.02987307,
        0.03446037,  0.03920482,  0.04541832,  0.05344991,  0.062497  ,
        0.07516488,  0.09013051,  0.11364832,  0.1407387 ,  0.18918521,
        0.27485483])



In [22]:

    
appres = data*np.pi*b*(b+a)/a
appres_obs = survey.dobs*np.pi*b*(b+a)/a
fig, ax = plt.subplots(1,2, figsize = (16, 4))
ax[1].semilogx(abhalf, appres, 'k.-')
ax[1].set_xscale('log')
ax[1].set_ylim(100., 180.)
ax[1].set_xlabel('AB/2')
ax[1].set_ylabel('Apparent resistivity ($\Omega m$)')
ax[1].grid(True)
# ax[1].text(100, 183, '(c)', fontsize = 16)
# ax[1].legend(('Observed', 'Predicted'), loc = 1)

dat = mesh.plotSlice((mapping*mtrue), normal='Y', ind = 9, ax = ax[0])
cb = plt.colorbar(dat[0], ax =ax[0])
ax[0].set_title("Vertical section")
cb.set_label("Conductivity (S/m)")
ax[0].set_xlabel('Easting (m)')
ax[0].set_ylabel('Depth (m)')
ax[0].set_xlim(-1000., 1000.)
ax[0].set_ylim(-500., 0.)
# ax[0].text(-1000, 20, '(b)')
# fig.savefig('DCfwd.png', dpi=200)









    Out[22]:





(-500.0, 0.0)

Step6: Run inversion



In [23]:

    
dmis = DataMisfit.l2_DataMisfit(survey)
reg = Regularization.Tikhonov(mesh,mapping=mapping)
opt = Optimization.InexactGaussNewton(maxIter=7,tolX=1e-15)
opt.remember('xc')
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
beta = Directives.BetaEstimate_ByEig(beta0_ratio=1e1)
betaSched = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
inv = Inversion.BaseInversion(invProb, directiveList=[beta,betaSched])



In [24]:

    
m0 = np.log(np.ones(problem.mapping.nP)*sighalf)
mopt = inv.run(m0)









    



SimPEG.InvProblem will set Regularization.mref to m0.
SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
                    ***Done using same solver as the problem***
SimPEG.l2_DataMisfit is creating default weightings for Wd.
============================ Inexact Gauss Newton ============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment   
-----------------------------------------------------------------------------
   0  2.61e-01  6.27e+03  4.48e+00  6.27e+03    8.80e+03      0              
   1  2.61e-01  5.22e+02  3.96e+00  5.23e+02    2.50e+03      0              
   2  5.23e-02  1.67e+02  3.91e+00  1.67e+02    8.21e+02      0              
   3  5.23e-02  7.93e+01  4.08e+00  7.95e+01    2.62e+02      0   Skip BFGS  
   4  1.05e-02  6.35e+01  4.25e+00  6.36e+01    1.71e+02      0   Skip BFGS  
   5  1.05e-02  3.01e+01  5.08e+00  3.02e+01    2.36e+02      0   Skip BFGS  
   6  2.09e-03  2.33e+01  5.32e+00  2.33e+01    1.57e+02      0   Skip BFGS  
   7  2.09e-03  1.01e+01  6.69e+00  1.01e+01    1.34e+02      0   Skip BFGS  
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.3195e+01 <= tolF*(1+|f0|) = 6.2702e+02
0 : |xc-x_last| = 7.9764e-01 <= tolX*(1+|x0|) = 2.6641e-14
0 : |proj(x-g)-x|    = 1.3417e+02 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 1.3417e+02 <= 1e3*eps       = 1.0000e-02
1 : maxIter   =       7    <= iter          =      7
------------------------- DONE! -------------------------



In [25]:

    
# matplotlib.rcParams.update({'font.size': 14, 'text.usetex': True, 'font.family': 'arial'})



In [26]:

    
appres = data*np.pi*b*(b+a)/a
appres_obs = survey.dobs*np.pi*b*(b+a)/a
appres_pred = invProb.dpred*np.pi*b*(b+a)/a
fig, ax = plt.subplots(1,2, figsize = (17, 6))
ax[0].plot(abhalf, appres_obs, 'k.-')
ax[0].plot(abhalf, appres_pred, 'r.-')
ax[0].set_xscale('log')
ax[0].set_ylim(100., 180.)
ax[0].set_xlabel('AB/2', fontsize=25)
ax[0].set_ylabel('Apparent resistivity ($\Omega m$)', fontsize=25)
ax[0].grid(True)
ax[0].legend(('Observed', 'Predicted'), loc = 1, fontsize=20)
ax[0].text(100., 181, '(a)', fontsize = 28)
ax[1].plot(1., 1., 'k', lw = 2)
ax[1].plot(1., 1., 'r', lw = 2)
ax[1].legend(('True', 'Predicted'), loc = 3, fontsize = 20)
Utils1D.plotLayer((np.exp(mopt)), mesh.vectorCCz, 'log', ax = ax[1], **{'lw':2, 'color':'r'})
Utils1D.plotLayer((np.exp(mtrue)), mesh.vectorCCz, 'log', showlayers=True, ax = ax[1], **{'lw':2})
ax[1].set_ylim(-500, 0)
ax[1].set_xlabel('Conductivity (S/m)', fontsize = 25)
ax[1].set_ylabel('Depth (m)', fontsize = 25)
ax[1].text(1e-3, 10., '(b)', fontsize = 28)
# fig.savefig('obspredDC.png', dpi=200)









    Out[26]:





<matplotlib.text.Text at 0x11049df90>



In [27]:

    
# fig, ax = plt.subplots(1,1, figsize = (7.5, 7))
# ax.plot(1., 1., 'k', lw = 2)
# ax.plot(1., 1., 'r', lw = 2)
# ax.legend(('True', 'Predicted'), loc = 3, fontsize = 16)
# Utils1D.plotLayer((np.exp(mopt)), mesh.vectorCCz, 'log', ax = ax, **{'lw':2, 'color':'r'})
# Utils1D.plotLayer((np.exp(mtrue)), mesh.vectorCCz, 'log', showlayers=True, ax = ax, **{'lw':2})
# ax.set_ylim(-500, 0)
# ax.set_xlabel('Conductivity (S/m)', fontsize = 20)
# ax.set_ylabel('Depth (m)', fontsize = 20)
# ax.text(1e-3, 10., '(b)', fontsize = 30)
# fig.savefig('obspred_dc1d_mod.png', dpi=200)



In [ ]: