mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-30 14:42:59 +08:00
Update example for constrained and compact-constrained MAG inversion
This commit is contained in:
D Fournier
@@ -31,34 +31,16 @@ M_xyz = driver.magnetizationModel
|
||||
|
||||
# Get index of the center
|
||||
midx = int(mesh.nCx/2)
|
||||
midy = int(mesh.nCy/2)
|
||||
midy = int(mesh.nCy/2)+1
|
||||
midz = int(mesh.nCz/2)
|
||||
|
||||
#%% Test script to impose static cells
|
||||
m0 = np.ones(mesh.nC)*1e-3
|
||||
|
||||
# Reshape the model in order to create a static block
|
||||
m0 = np.reshape(m0,(mesh.nCx,mesh.nCy,mesh.nCz), order = 'F')
|
||||
m0[midx,midy,midz] = 0.5
|
||||
m0 = mkvc(m0)
|
||||
|
||||
# Extract cells under topography and create new index for inactive
|
||||
m0 = m0[actv]
|
||||
ind_act = m0!=0.5
|
||||
|
||||
actvCells = Maps.InjectActiveCells(None, ind_act, 0.5, nC=nC)
|
||||
m0 = m0[ind_act]
|
||||
|
||||
#%% Plot obs data
|
||||
PF.Magnetics.plot_obs_2D(rxLoc,d, 'Observed Data')
|
||||
|
||||
#%% Run inversion
|
||||
prob = PF.Magnetics.Problem3D_Integral(mesh, mapping=actvCells, actInd=actv)
|
||||
prob = PF.Magnetics.Problem3D_Integral(mesh, mapping=idenMap, actInd=actv)
|
||||
prob.solverOpts['accuracyTol'] = 1e-4
|
||||
survey.pair(prob)
|
||||
|
||||
# Write out the predicted
|
||||
pred = prob.fields(m0)
|
||||
pred = prob.fields(driver.m0)
|
||||
PF.Magnetics.writeUBCobs('Pred.dat', survey, pred)
|
||||
|
||||
wr = np.sum(prob.G**2.,axis=0)**0.5 / mesh.vol[actv]
|
||||
@@ -72,8 +54,8 @@ plt.title('Distance weighting')
|
||||
plt.xlabel('x');plt.ylabel('z')
|
||||
plt.gca().set_aspect('equal', adjustable='box')
|
||||
|
||||
reg = Regularization.Simple(mesh, indActive=actv, mapping=actvCells)
|
||||
#reg.mref = m0*0
|
||||
reg = Regularization.Simple(mesh, indActive=actv, mapping=idenMap)
|
||||
reg.mref = driver.mref
|
||||
reg.wght = wr
|
||||
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
@@ -90,16 +72,16 @@ update_Jacobi = Directives.Update_lin_PreCond(onlyOnStart=True)
|
||||
inv = Inversion.BaseInversion(invProb, directiveList=[beta,target,beta_init,update_Jacobi])
|
||||
|
||||
# Run inversion
|
||||
mrec = inv.run(m0)
|
||||
mrec = inv.run(driver.m0)
|
||||
|
||||
m_out = actvMap*actvCells*mrec
|
||||
m_out = actvMap*mrec
|
||||
|
||||
# Write result
|
||||
Mesh.TensorMesh.writeModelUBC(mesh,'SimPEG_inv_l2l2.sus',m_out)
|
||||
|
||||
# Plot predicted
|
||||
pred = prob.fields(mrec)
|
||||
PF.Magnetics.plot_obs_2D(rxLoc,pred,'Predicted Data - l2 Inversion')
|
||||
#PF.Magnetics.plot_obs_2D(rxLoc,pred,wd,'Predicted Data')
|
||||
#PF.Magnetics.plot_obs_2D(rxLoc,(d-pred),wd,'Residual Data')
|
||||
|
||||
print "Final misfit:" + str(np.sum( ((d-pred)/wd)**2. ) )
|
||||
@@ -127,6 +109,166 @@ plt.gca().set_aspect('equal', adjustable='box')
|
||||
|
||||
ax = plt.subplot(212)
|
||||
mesh.plotSlice(m_out, ax = ax, normal = 'Y', ind=yslice, clim = (mrec.min(), mrec.max()))
|
||||
plt.title('Cross Section')
|
||||
plt.title('Smooth Unconstrained')
|
||||
plt.xlabel('x');plt.ylabel('z')
|
||||
plt.gca().set_aspect('equal', adjustable='box')
|
||||
plt.gca().set_aspect('equal', adjustable='box')
|
||||
|
||||
#%% Re-run inversion using a starting model with
|
||||
# static cells
|
||||
m0 = np.ones(mesh.nC)*1e-3
|
||||
val = 0.005
|
||||
|
||||
# Reshape the model in order to create a static block
|
||||
m0 = np.reshape(m0,(mesh.nCx,mesh.nCy,mesh.nCz), order = 'F')
|
||||
m0[midx-6,midy,midz+2] = val
|
||||
m0[midx+7,midy,midz+2] = val
|
||||
m0 = mkvc(m0)
|
||||
|
||||
# Extract cells under topography and create new index for inactive
|
||||
m0 = m0[actv]
|
||||
ind_act = m0!=val
|
||||
|
||||
actvCells = Maps.InjectActiveCells(None, ind_act, val, nC=nC)
|
||||
m0 = m0[ind_act]
|
||||
|
||||
# Change the mapping of the problem and run inversion
|
||||
prob.mapping = actvCells
|
||||
|
||||
# Write out the predicted
|
||||
pred = prob.fields(m0)
|
||||
PF.Magnetics.writeUBCobs('Pred.dat', survey, pred)
|
||||
|
||||
reg = Regularization.Simple(mesh, indActive=actv, mapping=actvCells)
|
||||
reg.mref = driver.mref[ind_act]
|
||||
reg.wght = wr
|
||||
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
dmis.Wd = 1/wd
|
||||
opt = Optimization.ProjectedGNCG(maxIter=10,lower=0.,upper=1., maxIterCG= 20, tolCG=1e-3)
|
||||
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
|
||||
|
||||
beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
|
||||
beta_init = Directives.BetaEstimate_ByEig()
|
||||
target = Directives.TargetMisfit()
|
||||
update_Jacobi = Directives.Update_lin_PreCond(onlyOnStart=True)
|
||||
|
||||
inv = Inversion.BaseInversion(invProb, directiveList=[beta,target,beta_init,update_Jacobi])
|
||||
|
||||
# Run inversion
|
||||
mrec = inv.run(m0)
|
||||
|
||||
m_out = actvMap*actvCells*mrec
|
||||
|
||||
# Write result
|
||||
Mesh.TensorMesh.writeModelUBC(mesh,'SimPEG_inv_l2l2_Constrained.sus',m_out)
|
||||
|
||||
# Plot predicted
|
||||
pred = prob.fields(mrec)
|
||||
#PF.Magnetics.plot_obs_2D(rxLoc,pred,'Predicted Data - l2 Inversion')
|
||||
#PF.Magnetics.plot_obs_2D(rxLoc,(d-pred),wd,'Residual Data')
|
||||
|
||||
print "Final misfit:" + str(np.sum( ((d-pred)/wd)**2. ) )
|
||||
|
||||
#%% Plot out a section of the model
|
||||
|
||||
yslice = midx
|
||||
m_out[m_out==-100] = np.nan
|
||||
|
||||
plt.figure()
|
||||
ax = plt.subplot(221)
|
||||
mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-5, clim = (mrec.min(), mrec.max()))
|
||||
plt.plot(np.array([mesh.vectorCCx[0],mesh.vectorCCx[-1]]), np.array([mesh.vectorCCy[yslice],mesh.vectorCCy[yslice]]),c='w',linestyle = '--')
|
||||
plt.title('Z: ' + str(mesh.vectorCCz[-5]) + ' m')
|
||||
plt.xlabel('x');plt.ylabel('z')
|
||||
plt.gca().set_aspect('equal', adjustable='box')
|
||||
|
||||
ax = plt.subplot(222)
|
||||
mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-8, clim = (mrec.min(), mrec.max()))
|
||||
plt.plot(np.array([mesh.vectorCCx[0],mesh.vectorCCx[-1]]), np.array([mesh.vectorCCy[yslice],mesh.vectorCCy[yslice]]),c='w',linestyle = '--')
|
||||
plt.title('Z: ' + str(mesh.vectorCCz[-8]) + ' m')
|
||||
plt.xlabel('x');plt.ylabel('z')
|
||||
plt.gca().set_aspect('equal', adjustable='box')
|
||||
|
||||
ax = plt.subplot(212)
|
||||
mesh.plotSlice(m_out, ax = ax, normal = 'Y', ind=yslice, clim = (mrec.min(), mrec.max()))
|
||||
plt.title('Smooth Constrained')
|
||||
plt.xlabel('x');plt.ylabel('z')
|
||||
plt.gca().set_aspect('equal', adjustable='box')
|
||||
|
||||
#%% Run one more round for sparsity
|
||||
phim = invProb.phi_m_last
|
||||
phid = invProb.phi_d
|
||||
|
||||
# Set parameters for sparsity
|
||||
reg = Regularization.Sparse(mesh, indActive = actv, mapping=actvCells)
|
||||
reg.recModel = mrec
|
||||
reg.mref = driver.mref[ind_act]
|
||||
reg.wght = wr
|
||||
reg.eps_p = eps_p
|
||||
reg.eps_q = eps_q
|
||||
reg.norms = driver.lpnorms
|
||||
|
||||
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
dmis.Wd = wd
|
||||
opt = Optimization.ProjectedGNCG(maxIter=10 , lower=0.,upper=1., maxIterCG= 10, tolCG = 1e-4)
|
||||
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta)
|
||||
|
||||
# Create inversion directives
|
||||
beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
|
||||
update_beta = Directives.Scale_Beta(tol = 0.05)
|
||||
target = Directives.TargetMisfit()
|
||||
IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid )
|
||||
update_Jacobi = Directives.Update_lin_PreCond(onlyOnStart=False)
|
||||
save_log = Directives.SaveOutputEveryIteration()
|
||||
save_log.fileName = 'LogName_blabla'
|
||||
|
||||
inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS,update_beta,update_Jacobi,save_log])
|
||||
|
||||
m0 = mrec
|
||||
|
||||
# Run inversion
|
||||
mrec = inv.run(m0)
|
||||
|
||||
m_out = actvMap*actvCells*mrec
|
||||
|
||||
# Write final model out.
|
||||
Mesh.TensorMesh.writeModelUBC(mesh,'SimPEG_inv_l0l2_Constrained.sus',m_out)
|
||||
|
||||
pred = prob.fields(mrec)
|
||||
|
||||
#%% Plot obs data
|
||||
#PF.Magnetics.plot_obs_2D(rxLoc,pred,'Predicted Data')
|
||||
#PF.Magnetics.plot_obs_2D(rxLoc,d,'Observed Data')
|
||||
print "Final misfit:" + str(np.sum( ((d-pred)/wd)**2. ) )
|
||||
#%% Plot out a section of the model
|
||||
|
||||
yslice = midy
|
||||
|
||||
m_out[m_out==-100] = np.nan
|
||||
|
||||
plt.figure()
|
||||
ax = plt.subplot(221)
|
||||
mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-5, clim = (mrec.min(), mrec.max()))
|
||||
plt.plot(np.array([mesh.vectorCCx[0],mesh.vectorCCx[-1]]), np.array([mesh.vectorCCy[yslice],mesh.vectorCCy[yslice]]),c='w',linestyle = '--')
|
||||
plt.title('Z: ' + str(mesh.vectorCCz[-5]) + ' m')
|
||||
plt.xlabel('x');plt.ylabel('z')
|
||||
plt.gca().set_aspect('equal', adjustable='box')
|
||||
|
||||
ax = plt.subplot(222)
|
||||
mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-8, clim = (mrec.min(), mrec.max()))
|
||||
plt.plot(np.array([mesh.vectorCCx[0],mesh.vectorCCx[-1]]), np.array([mesh.vectorCCy[yslice],mesh.vectorCCy[yslice]]),c='w',linestyle = '--')
|
||||
plt.title('Z: ' + str(mesh.vectorCCz[-8]) + ' m')
|
||||
plt.xlabel('x');plt.ylabel('z')
|
||||
plt.gca().set_aspect('equal', adjustable='box')
|
||||
|
||||
ax = plt.subplot(212)
|
||||
mesh.plotSlice(m_out, ax = ax, normal = 'Y', ind=yslice, clim = (mrec.min(), mrec.max()))
|
||||
plt.title('Compact Constrained')
|
||||
plt.xlabel('x');plt.ylabel('z')
|
||||
plt.gca().set_aspect('equal', adjustable='box')
|
||||
|
||||
|
||||
plt.show()
|
||||
Reference in New Issue
Block a user