.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "content/examples/08-vrm/plot_inv_vrm_eq.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_content_examples_08-vrm_plot_inv_vrm_eq.py: Method of Equivalent Sources for Removing VRM Responses ======================================================= Here, we use an equivalent source inversion to remove the VRM response from TEM data collected by a small coincident loop system. The data being inverted are the same as in the forward modeling example. To remove the VRM signal we: 1. invert the late time data to recover an equivalent source surface layer of cells. 2. use the recovered model to predict the VRM response at all times 3. subtract the predicted VRM response from the observed data .. GENERATED FROM PYTHON SOURCE LINES 15-18 Import modules -------------- .. GENERATED FROM PYTHON SOURCE LINES 18-37 .. code-block:: Python from simpeg.electromagnetics import viscous_remanent_magnetization as VRM import numpy as np import discretize from simpeg import ( utils, maps, data_misfit, directives, optimization, regularization, inverse_problem, inversion, data, ) import matplotlib.pyplot as plt import matplotlib as mpl .. GENERATED FROM PYTHON SOURCE LINES 38-41 Defining the mesh ----------------- .. GENERATED FROM PYTHON SOURCE LINES 41-48 .. code-block:: Python cs, ncx, ncy, ncz, npad = 2.0, 35, 35, 20, 5 hx = [(cs, npad, -1.3), (cs, ncx), (cs, npad, 1.3)] hy = [(cs, npad, -1.3), (cs, ncy), (cs, npad, 1.3)] hz = [(cs, npad, -1.3), (cs, ncz), (cs, npad, 1.3)] mesh = discretize.TensorMesh([hx, hy, hz], "CCC") .. GENERATED FROM PYTHON SOURCE LINES 49-56 Defining the true model ----------------------- Create xi model (amalgamated magnetic property). Here the model is made by summing a set of 3D Gaussian distributions. And only active cells have a model value. .. GENERATED FROM PYTHON SOURCE LINES 56-81 .. code-block:: Python topoCells = mesh.gridCC[:, 2] < 0.0 # define topography xyzc = mesh.gridCC[topoCells, :] c = 2 * np.pi * 8**2 pc = np.r_[4e-4, 4e-4, 4e-4, 6e-4, 8e-4, 6e-4, 8e-4, 8e-4] x_0 = np.r_[50.0, -50.0, -40.0, -20.0, -15.0, 20.0, -10.0, 25.0] y_0 = np.r_[0.0, 0.0, 40.0, 10.0, -20.0, 15.0, 0.0, 0.0] z_0 = np.r_[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] var_x = c * np.r_[3.0, 3.0, 3.0, 1.0, 3.0, 0.5, 0.1, 0.1] var_y = c * np.r_[20.0, 20.0, 1.0, 1.0, 0.4, 0.5, 0.1, 0.4] var_z = c * np.r_[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] xi_true = np.zeros(np.shape(xyzc[:, 0])) for ii in range(0, 8): xi_true += ( pc[ii] * np.exp(-((xyzc[:, 0] - x_0[ii]) ** 2) / var_x[ii]) * np.exp(-((xyzc[:, 1] - y_0[ii]) ** 2) / var_y[ii]) * np.exp(-((xyzc[:, 2] - z_0[ii]) ** 2) / var_z[ii]) ) xi_true += 1e-5 .. GENERATED FROM PYTHON SOURCE LINES 82-89 Survey ------ Here we must set the transmitter waveform, which defines the off-time decay of the VRM response. Next we define the sources, receivers and time channels for the survey. Our example is similar to an EM-63 survey. .. GENERATED FROM PYTHON SOURCE LINES 89-113 .. code-block:: Python waveform = VRM.waveforms.StepOff() times = np.logspace(-5, -2, 31) # Observation times x, y = np.meshgrid(np.linspace(-30, 30, 21), np.linspace(-30, 30, 21)) z = 0.5 * np.ones(x.shape) loc = np.c_[utils.mkvc(x), utils.mkvc(y), utils.mkvc(z)] # Src and Rx Locations source_listVRM = [] for pp in range(0, loc.shape[0]): loc_pp = np.reshape(loc[pp, :], (1, 3)) receiver_listVRM = [ VRM.Rx.Point(loc_pp, times=times, field_type="dbdt", orientation="z") ] source_listVRM.append( VRM.Src.MagDipole( receiver_listVRM, utils.mkvc(loc[pp, :]), [0.0, 0.0, 0.01], waveform ) ) survey_vrm = VRM.Survey(source_listVRM) .. GENERATED FROM PYTHON SOURCE LINES 114-122 Forward Simulation ------------------ Here we predict data by solving the forward problem. For the VRM problem, we use a sensitivity refinement strategy for cells # that are proximal to transmitters. This is controlled through the *refinement_factor* and *refinement_distance* properties. .. GENERATED FROM PYTHON SOURCE LINES 122-164 .. code-block:: Python # Defining the problem problem_vrm = VRM.Simulation3DLinear( mesh, survey=survey_vrm, active_cells=topoCells, refinement_factor=3, refinement_distance=[1.25, 2.5, 3.75], ) # Predict VRM response fields_vrm = problem_vrm.dpred(xi_true) # Add an artificial TEM response. An analytic solution for the response near # the surface of a conductive half-space (Nabighian, 1979) is scaled at each # location to provide lateral variability in the TEM response. n_times = len(times) n_loc = loc.shape[0] sig = 1e-1 mu0 = 4 * np.pi * 1e-7 fields_tem = -(sig**1.5) * mu0**2.5 * times**-2.5 / (20 * np.pi**1.5) fields_tem = np.kron(np.ones(n_loc), fields_tem) c = ( np.exp(-((loc[:, 0] - 10) ** 2) / (25**2)) * np.exp(-((loc[:, 1] - 20) ** 2) / (35**2)) + np.exp(-((loc[:, 0] + 20) ** 2) / (20**2)) * np.exp(-((loc[:, 1] + 20) ** 2) / (40**2)) + 1.5 * np.exp(-((loc[:, 0] - 25) ** 2) / (10**2)) * np.exp(-((loc[:, 1] + 25) ** 2) / (10**2)) + 0.25 ) c = np.kron(c, np.ones(n_times)) fields_tem = c * fields_tem fields_tot = fields_tem + fields_vrm fields_tot = fields_tot + 0.05 * np.abs(fields_tot) * np.random.normal( size=fields_tot.shape ) .. rst-class:: sphx-glr-script-out .. code-block:: none CREATING T MATRIX CREATING A MATRIX .. GENERATED FROM PYTHON SOURCE LINES 165-172 Inverse Problem --------------- Here, we invert late-time data to recover an equivalent source model. To recover the equivalent source model, only cells at the surface are set as active in the inversion. .. GENERATED FROM PYTHON SOURCE LINES 172-220 .. code-block:: Python # Define problem # survey_inv = VRM.Survey(source_listVRM) actCells = (mesh.gridCC[:, 2] < 0.0) & (mesh.gridCC[:, 2] > -2.0) problem_inv = VRM.Simulation3DLinear( mesh, survey=survey_vrm, active_cells=actCells, refinement_factor=3, refinement_distance=[1.25, 2.5, 3.75], ) survey_vrm.set_active_interval(1e-3, 1e-2) dobs = fields_tot[survey_vrm.t_active] rel_err = 0.05 eps = 1e-11 data_vrm = data.Data( dobs=dobs, survey=survey_vrm, relative_error=rel_err, noise_floor=eps ) # Setup and run inversion dmis = data_misfit.L2DataMisfit(simulation=problem_inv, data=data_vrm) w = utils.mkvc((np.sum(np.array(problem_inv.A) ** 2, axis=0))) ** 0.5 w = w / np.max(w) w = w reg = regularization.Smallness( mesh=mesh, active_cells=actCells, weights={"cell_weights": w} ) opt = optimization.ProjectedGNCG( maxIter=20, lower=0.0, upper=1e-2, maxIterLS=20, tolCG=1e-4 ) invProb = inverse_problem.BaseInvProblem(dmis, reg, opt) directives = [ directives.BetaSchedule(coolingFactor=2, coolingRate=1), directives.TargetMisfit(), ] inv = inversion.BaseInversion(invProb, directiveList=directives) xi_0 = 1e-3 * np.ones(actCells.sum()) xi_rec = inv.run(xi_0) # Predict VRM response at all times for recovered model survey_vrm.set_active_interval(0.0, 1.0) fields_pre = problem_inv.dpred(xi_rec) .. rst-class:: sphx-glr-script-out .. code-block:: none CREATING A MATRIX Running inversion with SimPEG v0.23.0 simpeg.InvProblem will set Regularization.reference_model to m0. /home/vsts/work/1/s/simpeg/simulation.py:197: DefaultSolverWarning: Using the default solver: Pardiso. If you would like to suppress this notification, add warnings.filterwarnings('ignore', simpeg.utils.solver_utils.DefaultSolverWarning) to your script. simpeg.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv. ***Done using same Solver, and solver_opts as the Simulation3DLinear 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 CREATING T MATRIX 0 1.00e+00 3.54e+05 0.00e+00 3.54e+05 9.97e-02 0 ------------------------- STOP! ------------------------- 1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 3.5395e+04 1 : |xc-x_last| = 2.3607e-02 <= tolX*(1+|x0|) = 1.0450e-01 1 : |proj(x-g)-x| = 9.2321e-02 <= tolG = 1.0000e-01 0 : |proj(x-g)-x| = 9.2321e-02 <= 1e3*eps = 1.0000e-02 0 : maxIter = 20 <= iter = 1 ------------------------- DONE! ------------------------- .. GENERATED FROM PYTHON SOURCE LINES 221-224 Plotting -------- .. GENERATED FROM PYTHON SOURCE LINES 224-359 .. code-block:: Python fields_tot = np.reshape(fields_tot, (n_loc, n_times)) fields_vrm = np.reshape(fields_vrm, (n_loc, n_times)) fields_tem = np.reshape(fields_tem, (n_loc, n_times)) fields_pre = np.reshape(fields_pre, (n_loc, n_times)) Fig = plt.figure(figsize=(10, 10)) font_size = 12 # Plot models invMap = maps.InjectActiveCells(mesh, actCells, 0.0) # Maps to mesh topoMap = maps.InjectActiveCells(mesh, topoCells, 0.0) max_val = np.max(np.r_[xi_true, xi_rec]) ax1 = 3 * [None] cplot1 = 2 * [None] xi_mod = [xi_true, xi_rec] map_mod = [topoMap, invMap] titlestr1 = ["True Model (z = 0 m)", "Equivalent Source Model"] for qq in range(0, 2): ax1[qq] = Fig.add_axes([0.15 + 0.35 * qq, 0.7, 0.25, 0.25]) cplot1[qq] = mesh.plot_slice( map_mod[qq] * xi_mod[qq], ind=int((ncz + 2 * npad) / 2 - 1), ax=ax1[qq], grid=True, pcolor_opts={"cmap": "gist_heat_r"}, ) cplot1[qq][0].set_clim((0.0, max_val)) ax1[qq].set_xlabel("X [m]", fontsize=font_size) ax1[qq].set_ylabel("Y [m]", fontsize=font_size, labelpad=-5) ax1[qq].tick_params(labelsize=font_size - 2) ax1[qq].set_title(titlestr1[qq], fontsize=font_size + 2) ax1[2] = Fig.add_axes([0.78, 0.7, 0.01, 0.25]) norm = mpl.colors.Normalize(vmin=0.0, vmax=max_val) cbar14 = mpl.colorbar.ColorbarBase( ax1[2], cmap=mpl.cm.gist_heat_r, norm=norm, orientation="vertical" ) cbar14.set_label( r"$\Delta \chi /$ln$(\lambda_2 / \lambda_1 )$ [SI]", rotation=270, labelpad=15, size=font_size, ) # Plot decays N = x.shape[0] ax2 = 2 * [None] for qq in range(0, 2): ax2[qq] = Fig.add_axes([0.1 + 0.45 * qq, 0.36, 0.35, 0.26]) k = int((N**2 - 1) / 2 - 3 * N * (-1) ** qq) di_tot = utils.mkvc(np.abs(fields_tot[k, :])) di_pre = utils.mkvc(np.abs(fields_vrm[k, :])) di_tem = utils.mkvc(np.abs(fields_tem[k, :])) ax2[qq].loglog(times, di_tot, "k.-") ax2[qq].loglog(times, di_tem, "r.-") ax2[qq].loglog(times, di_pre, "b.-") ax2[qq].loglog(times, np.abs(di_tot - di_pre), "g.-") ax2[qq].set_xlabel("t [s]", fontsize=font_size, labelpad=-10) if qq == 0: ax2[qq].set_ylabel("|dBz/dt| [T/s]", fontsize=font_size) else: ax2[qq].axes.get_yaxis().set_visible(False) ax2[qq].tick_params(labelsize=font_size - 2) ax2[qq].set_xbound(np.min(times), np.max(times)) ax2[qq].set_ybound(1.2 * np.max(di_tot), 1e-5 * np.max(di_tot)) titlestr2 = ( "Decay at X = " + "{:.2f}".format(loc[k, 0]) + " m and Y = " + "{:.2f}".format(loc[k, 1]) + " m" ) ax2[qq].set_title(titlestr2, fontsize=font_size + 2) if qq == 0: ax2[qq].text( 1.2e-5, 54 * np.max(di_tot) / 1e5, "Observed", fontsize=font_size, color="k" ) ax2[qq].text( 1.2e-5, 18 * np.max(di_tot) / 1e5, "True TEM", fontsize=font_size, color="r" ) ax2[qq].text( 1.2e-5, 6 * np.max(di_tot) / 1e5, "Predicted VRM", fontsize=font_size, color="b", ) ax2[qq].text( 1.2e-5, 2 * np.max(di_tot) / 1e5, "Recovered TEM", fontsize=font_size, color="g", ) # Plot anomalies d = [ np.reshape(np.abs(fields_tot[:, 10]), (N, N)), np.reshape(np.abs(fields_tem[:, 10]), (N, N)), np.reshape(np.abs(fields_tot[:, 10] - fields_pre[:, 10]), (N, N)), ] min_val = np.min(np.r_[d[0], d[1], d[2]]) max_val = np.max(np.r_[d[0], d[1], d[2]]) ax3 = 4 * [None] cplot3 = 3 * [None] title_str = ["Observed at t=", "True TEM at t=", "Recov. TEM at t="] for qq in range(0, 3): ax3[qq] = Fig.add_axes([0.07 + 0.28 * qq, 0.05, 0.24, 0.24]) cplot3[qq] = ax3[qq].contourf(x, y, d[qq].T, 40, cmap="magma_r") ax3[qq].set_xticks(np.linspace(-30, 30, 7)) ax3[qq].set_xlabel("X [m]", fontsize=font_size) if qq == 0: ax3[qq].scatter(x, y, color=(0, 0, 0), s=4) ax3[qq].set_ylabel("Y [m]", fontsize=font_size, labelpad=-12) else: ax3[qq].axes.get_yaxis().set_visible(False) ax3[qq].tick_params(labelsize=font_size - 2) ax3[qq].set_xbound(np.min(x), np.max(x)) ax3[qq].set_ybound(np.min(y), np.max(y)) titlestr3 = title_str[qq] + "{:.1e}".format(times[10]) + " s" ax3[qq].set_title(titlestr3, fontsize=font_size + 2) ax3[3] = Fig.add_axes([0.88, 0.05, 0.01, 0.24]) norm = mpl.colors.Normalize(vmin=min_val, vmax=max_val) cbar34 = mpl.colorbar.ColorbarBase( ax3[3], cmap=mpl.cm.magma_r, norm=norm, orientation="vertical", format="%.1e" ) cbar34.set_label("dBz/dt [T/s]", rotation=270, size=font_size, labelpad=15) plt.show() .. image-sg:: /content/examples/08-vrm/images/sphx_glr_plot_inv_vrm_eq_001.png :alt: True Model (z = 0 m), Equivalent Source Model, Decay at X = -9.00 m and Y = 0.00 m, Decay at X = 9.00 m and Y = 0.00 m, Observed at t=1.0e-04 s, True TEM at t=1.0e-04 s, Recov. TEM at t=1.0e-04 s :srcset: /content/examples/08-vrm/images/sphx_glr_plot_inv_vrm_eq_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 17.893 seconds) **Estimated memory usage:** 288 MB .. _sphx_glr_download_content_examples_08-vrm_plot_inv_vrm_eq.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_inv_vrm_eq.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_inv_vrm_eq.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_inv_vrm_eq.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_