Note
Go to the end to download the full example code.
Straight Ray with Volume Data Misfit Term#
Based on the SEG abstract Heagy, Cockett and Oldenburg, 2014.
Heagy, L. J., Cockett, A. R., & Oldenburg, D. W. (2014, August 5). Parametrized Inversion Framework for Proppant Volume in a Hydraulically Fractured Reservoir. SEG Technical Program Expanded Abstracts 2014. Society of Exploration Geophysicists. doi:10.1190/segam2014-1639.1
This example is a simple joint inversion that consists of a
data misfit for the tomography problem
data misfit for the volume of the inclusions (uses the effective medium theory mapping)
model regularization
True Volume: 0.6240000000000001
/home/vsts/work/1/s/examples/20-published/plot_tomo_joint_with_volume.py:157: FutureWarning:
The defaults for ProjectedGNCG will change in SimPEG 0.26.0. If you want to maintain the previous behavior, explicitly set 'cg_atol=1E-3' and 'cg_rtol=0.0'.
Running inversion with SimPEG v0.24.1.dev29+g256406f92
simpeg.InvProblem will set Regularization.reference_model to m0.
simpeg.InvProblem will set Regularization.reference_model to m0.
simpeg.InvProblem will set Regularization.reference_model to m0.
simpeg.InvProblem will set Regularization.reference_model to m0.
simpeg.InvProblem will set Regularization.reference_model to m0.
simpeg.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using the default solver Pardiso and no solver_opts.***
model has any nan: 0
================================================= Projected GNCG =================================================
# beta phi_d phi_m f |proj(x-g)-x| LS iterCG CG |Ax-b|/|b| CG |Ax-b| Comment
-----------------------------------------------------------------------------------------------------------------
x0 has any nan: 0
0 2.50e-01 8.48e+04 0.00e+00 8.48e+04 3.83e+01 0 0 inf inf
1 2.50e-01 4.40e+03 9.97e-01 4.40e+03 3.67e+01 0 5 4.80e-03 1.03e+01
/home/vsts/work/1/s/simpeg/maps/_property_maps.py:1504: UserWarning:
Maximum number of iterations reached
2 2.50e-01 6.30e+02 8.24e-01 6.31e+02 3.18e+01 0 5 1.73e-03 6.21e-02
3 2.50e-01 1.28e+02 6.79e-01 1.28e+02 2.13e+01 2 5 1.24e+00 3.69e+02
4 2.50e-01 6.83e+01 6.74e-01 6.85e+01 1.19e+01 0 5 4.44e-01 3.79e+01
5 2.50e-01 2.50e+01 6.77e-01 2.52e+01 1.15e+01 0 5 2.93e-01 4.65e+01
6 2.50e-01 1.68e+01 6.72e-01 1.70e+01 1.02e+01 0 5 7.92e-01 4.39e+01 Skip BFGS
7 2.50e-01 1.30e+01 6.68e-01 1.31e+01 8.89e+00 0 5 1.03e+00 3.25e+01 Skip BFGS
8 2.50e-01 1.17e+01 6.68e-01 1.18e+01 9.37e+00 0 5 1.02e+00 2.22e+01
9 2.50e-01 9.29e+00 6.66e-01 9.46e+00 7.88e+00 0 5 2.56e-01 1.75e+01
10 2.50e-01 7.81e+00 6.69e-01 7.97e+00 8.15e+00 0 5 8.62e-01 1.59e+01
11 2.50e-01 5.19e+00 6.69e-01 5.36e+00 6.35e+00 0 5 3.48e-01 1.52e+01
12 2.50e-01 4.57e+00 6.70e-01 4.74e+00 6.10e+00 0 5 7.09e-01 1.38e+01
13 2.50e-01 3.56e+00 6.74e-01 3.73e+00 5.98e+00 0 5 1.20e+00 1.74e+01
14 2.50e-01 3.13e+00 6.73e-01 3.30e+00 5.53e+00 0 5 4.52e-01 9.78e+00
15 2.50e-01 2.75e+00 6.75e-01 2.92e+00 5.35e+00 0 5 1.60e+00 1.53e+01 Skip BFGS
------------------------- STOP! -------------------------
1 : |fc-fOld| = 3.7923e-01 <= tolF*(1+|f0|) = 8.4754e+03
0 : |xc-x_last| = 2.4366e-01 <= tolX*(1+|x0|) = 1.0000e-01
0 : |proj(x-g)-x| = 5.3451e+00 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 5.3451e+00 <= 1e3*eps = 1.0000e-02
1 : maxIter = 15 <= iter = 15
------------------------- DONE! -------------------------
Total recovered volume (no vol misfit term in inversion): 0.018812203265143886
/home/vsts/work/1/s/examples/20-published/plot_tomo_joint_with_volume.py:172: FutureWarning:
The defaults for ProjectedGNCG will change in SimPEG 0.26.0. If you want to maintain the previous behavior, explicitly set 'cg_atol=1E-3' and 'cg_rtol=0.0'.
Running inversion with SimPEG v0.24.1.dev29+g256406f92
simpeg.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using the default solver Pardiso and no solver_opts.***
model has any nan: 0
================================================= Projected GNCG =================================================
# beta phi_d phi_m f |proj(x-g)-x| LS iterCG CG |Ax-b|/|b| CG |Ax-b| Comment
-----------------------------------------------------------------------------------------------------------------
x0 has any nan: 0
0 2.50e-01 1.20e+05 0.00e+00 1.20e+05 1.58e+01 0 0 inf inf
1 2.50e-01 7.50e+03 8.51e-01 7.50e+03 2.72e+01 0 5 6.24e-02 9.16e+01
2 2.50e-01 7.48e+03 8.37e-01 7.48e+03 2.69e+01 5 5 3.92e-01 2.04e+02
3 2.50e-01 7.48e+03 8.19e-01 7.48e+03 2.65e+01 4 5 3.30e-01 1.97e+02 Skip BFGS
4 2.50e-01 7.47e+03 8.13e-01 7.47e+03 2.63e+01 5 5 5.71e-01 2.44e+02 Skip BFGS
5 2.50e-01 7.47e+03 8.11e-01 7.47e+03 2.63e+01 7 5 6.33e-01 2.62e+02 Skip BFGS
6 2.50e-01 7.47e+03 8.08e-01 7.47e+03 2.62e+01 6 5 5.57e-01 2.18e+02 Skip BFGS
------------------------------------------------------------------
0 : ft = 7.4758e+03 <= alp*descent = 7.4700e+03
1 : maxIterLS = 10 <= iterLS = 10
------------------------- End Linesearch -------------------------
The linesearch got broken. Boo.
Total volume (vol misfit term in inversion): 0.027797507446205082
import numpy as np
import scipy.sparse as sp
import matplotlib.pyplot as plt
from simpeg.seismic import straight_ray_tomography as tomo
import discretize
from simpeg import (
maps,
utils,
regularization,
optimization,
inverse_problem,
inversion,
data_misfit,
objective_function,
)
class Volume(objective_function.BaseObjectiveFunction):
r"""
A regularization on the volume integral of the model
.. math::
\phi_v = || \int_V m dV - \text{knownVolume} ||^2
"""
def __init__(self, mesh, knownVolume=0.0, **kwargs):
self.mesh = mesh
self.knownVolume = knownVolume
super().__init__(**kwargs)
@property
def knownVolume(self):
"""known volume"""
return self._knownVolume
@knownVolume.setter
def knownVolume(self, value):
self._knownVolume = utils.validate_float("knownVolume", value, min_val=0.0)
def __call__(self, m):
return (self.estVol(m) - self.knownVolume) ** 2
def estVol(self, m):
return np.inner(self.mesh.cell_volumes, m)
def deriv(self, m):
# return (self.mesh.cell_volumes * np.inner(self.mesh.cell_volumes, m))
return (
2
* self.mesh.cell_volumes
* (self.knownVolume - np.inner(self.mesh.cell_volumes, m))
) # factor of 2 from deriv of ||estVol - knownVol||^2
def deriv2(self, m, v=None):
if v is not None:
return 2 * utils.mkvc(
self.mesh.cell_volumes * np.inner(self.mesh.cell_volumes, v)
)
else:
# TODO: this is inefficent. It is a fully dense matrix
return 2 * sp.csc_matrix(
np.outer(self.mesh.cell_volumes, self.mesh.cell_volumes)
)
def run(plotIt=True):
nC = 40
de = 1.0
h = np.ones(nC) * de / nC
M = discretize.TensorMesh([h, h])
y = np.linspace(M.cell_centers_y[0], M.cell_centers_x[-1], int(np.floor(nC / 4)))
rlocs = np.c_[0 * y + M.cell_centers_x[-1], y]
rx = tomo.Rx(rlocs)
source_list = [
tomo.Src(location=np.r_[M.cell_centers_x[0], yi], receiver_list=[rx])
for yi in y
]
# phi model
phi0 = 0
phi1 = 0.65
phitrue = utils.model_builder.create_block_in_wholespace(
M.gridCC, [0.4, 0.6], [0.6, 0.4], [phi1, phi0]
)
knownVolume = np.sum(phitrue * M.cell_volumes)
print("True Volume: {}".format(knownVolume))
# Set up true conductivity model and plot the model transform
sigma0 = np.exp(1)
sigma1 = 1e4
if plotIt:
fig, ax = plt.subplots(1, 1)
sigmaMapTest = maps.SelfConsistentEffectiveMedium(
nP=1000, sigma0=sigma0, sigma1=sigma1, rel_tol=1e-1, maxIter=150
)
testphis = np.linspace(0.0, 1.0, 1000)
sigetest = sigmaMapTest * testphis
ax.semilogy(testphis, sigetest)
ax.set_title("Model Transform")
ax.set_xlabel(r"$\varphi$")
ax.set_ylabel(r"$\sigma$")
sigmaMap = maps.SelfConsistentEffectiveMedium(M, sigma0=sigma0, sigma1=sigma1)
# scale the slowness so it is on a ~linear scale
slownessMap = maps.LogMap(M) * sigmaMap
# set up the problem and survey
survey = tomo.Survey(source_list)
problem = tomo.Simulation(M, survey=survey, slownessMap=slownessMap)
if plotIt:
_, ax = plt.subplots(1, 1)
cb = plt.colorbar(M.plot_image(phitrue, ax=ax)[0], ax=ax)
survey.plot(ax=ax)
cb.set_label(r"$\varphi$")
# get observed data
data = problem.make_synthetic_data(phitrue, relative_error=0.03, add_noise=True)
dpred = problem.dpred(np.zeros(M.nC))
# objective function pieces
reg = regularization.WeightedLeastSquares(M)
dmis = data_misfit.L2DataMisfit(simulation=problem, data=data)
dmisVol = Volume(mesh=M, knownVolume=knownVolume)
beta = 0.25
maxIter = 15
# without the volume regularization
opt = optimization.ProjectedGNCG(maxIter=maxIter, lower=0.0, upper=1.0)
opt.remember("xc")
invProb = inverse_problem.BaseInvProblem(dmis, reg, opt, beta=beta)
inv = inversion.BaseInversion(invProb)
mopt1 = inv.run(np.zeros(M.nC) + 1e-16)
print(
"\nTotal recovered volume (no vol misfit term in inversion): "
"{}".format(dmisVol(mopt1))
)
# with the volume regularization
vol_multiplier = 9e4
reg2 = reg
dmis2 = dmis + vol_multiplier * dmisVol
opt2 = optimization.ProjectedGNCG(maxIter=maxIter, lower=0.0, upper=1.0)
opt2.remember("xc")
invProb2 = inverse_problem.BaseInvProblem(dmis2, reg2, opt2, beta=beta)
inv2 = inversion.BaseInversion(invProb2)
mopt2 = inv2.run(np.zeros(M.nC) + 1e-16)
print("\nTotal volume (vol misfit term in inversion): {}".format(dmisVol(mopt2)))
# plot results
if plotIt:
fig, ax = plt.subplots(1, 1)
ax.plot(data.dobs)
ax.plot(dpred)
ax.plot(problem.dpred(mopt1), "o")
ax.plot(problem.dpred(mopt2), "s")
ax.legend(["dobs", "dpred0", "dpred w/o Vol", "dpred with Vol"])
fig, ax = plt.subplots(1, 3, figsize=(16, 4))
im0 = M.plot_image(phitrue, ax=ax[0])[0]
im1 = M.plot_image(mopt1, ax=ax[1])[0]
im2 = M.plot_image(mopt2, ax=ax[2])[0]
for im in [im0, im1, im2]:
im.set_clim([0.0, phi1])
plt.colorbar(im0, ax=ax[0])
plt.colorbar(im1, ax=ax[1])
plt.colorbar(im2, ax=ax[2])
ax[0].set_title("true, vol: {:1.3e}".format(knownVolume))
ax[1].set_title(
"recovered(no Volume term), vol: {:1.3e} ".format(dmisVol(mopt1))
)
ax[2].set_title(
"recovered(with Volume term), vol: {:1.3e} ".format(dmisVol(mopt2))
)
plt.tight_layout()
if __name__ == "__main__":
run()
plt.show()
Total running time of the script: (0 minutes 26.532 seconds)
Estimated memory usage: 290 MB