"""SimPEG 3-D magnetic susceptibility inversion — L2 vs sparse regularization. Two dipping tabular bodies beneath undulating topography. Forward-model TMI on an airborne grid, add noise consistent with the assigned uncertainty, then invert the SAME data twice: smooth (L2) and compact (sparse, L0 on the model). Demonstrates octree meshing, the integral (geoana/numba) formulation, sensitivity weighting, and the recovered-model difference the regularization choice makes. """ import json, time import numpy as np from discretize import TreeMesh from discretize.utils import active_from_xyz import simpeg from simpeg import (maps, data, data_misfit, regularization, optimization, inverse_problem, directives, inversion) from simpeg.potential_fields import magnetics as mag t0 = time.time() rng = np.random.default_rng(7) log = {"versions": {"simpeg": simpeg.__version__, "numpy": np.__version__}, "formulation": "magnetics.Simulation3DIntegral (geoana engine, numba parallel)"} # ---------------------------------------------------------------- topography def topo_z(x, y): return 40.0*np.exp(-((x-150.)**2+(y+100.)**2)/(2*260.**2)) \ - 25.0*np.exp(-((x+220.)**2+(y-180.)**2)/(2*200.**2)) tx = np.linspace(-900, 900, 101); TX, TY = np.meshgrid(tx, tx) topo = np.c_[TX.ravel(), TY.ravel(), topo_z(TX.ravel(), TY.ravel())] # ---------------------------------------------------------------- survey sx = np.arange(-500, 501, 50.0) # 50 m station / line spacing SX, SY = np.meshgrid(sx, sx) rx_loc = np.c_[SX.ravel(), SY.ravel(), topo_z(SX.ravel(), SY.ravel()) + 40.0] # 40 m drape rx = mag.Point(rx_loc, components="tmi") H0_amp, H0_inc, H0_dec = 50000.0, 65.0, 25.0 # nT, degrees src = mag.UniformBackgroundField(receiver_list=[rx], amplitude=H0_amp, inclination=H0_inc, declination=H0_dec) survey = mag.Survey(src) log["survey"] = {"nD": int(survey.nD), "n_stations": int(rx_loc.shape[0]), "station_spacing_m": 50.0, "drape_height_m": 40.0, "H0_nT": H0_amp, "inclination_deg": H0_inc, "declination_deg": H0_dec} # ---------------------------------------------------------------- octree mesh dh, nb = 25.0, 2**6 # 25 m core cells, 1600 m per axis # nb must give an extent that CONTAINS the survey (+-500 m) with padding: # 2**6 * 25 m = 1600 m -> x,y in [-800, 800]. 2**5 would span only +-400 m # and put every receiver outside the mesh. mesh = TreeMesh([[(dh, nb)], [(dh, nb)], [(dh, nb)]], origin="CCN") mesh.refine_surface(topo, level=-1, padding_cells_by_level=[4, 4, 4], finalize=False) mesh.refine_points(rx_loc, level=-1, padding_cells_by_level=[6, 6, 4], finalize=False) mesh.finalize() active = active_from_xyz(mesh, topo) nC = int(active.sum()) cc = mesh.cell_centers[active, :] log["mesh"] = {"type": "TreeMesh (octree)", "n_cells": int(mesh.n_cells), "n_active": nC, "core_cell_size_m": dh, "extent_x": [float(mesh.nodes_x.min()), float(mesh.nodes_x.max())], "extent_y": [float(mesh.nodes_y.min()), float(mesh.nodes_y.max())], "extent_z": [float(mesh.nodes_z.min()), float(mesh.nodes_z.max())]} # ---------------------------------------------------------------- true model # Two dipping tabular bodies (dykes), different susceptibility and dip sense. chi_true = np.zeros(nC) x, y, z = cc[:,0], cc[:,1], cc[:,2] zs = topo_z(x, y) # local surface elevation d1 = z - (zs - 120.0) - 0.55*(x + 120.0) # dips ~29 deg toward -x b1 = (np.abs(d1) < 55.0) & (np.abs(y + 60.0) < 190.0) & (z < zs - 40.0) & (z > zs - 420.0) chi_true[b1] = 0.06 d2 = z - (zs - 150.0) + 0.42*(x - 260.0) # opposite dip b2 = (np.abs(d2) < 45.0) & (np.abs(y - 190.0) < 150.0) & (z < zs - 60.0) & (z > zs - 400.0) chi_true[b2] = 0.035 log["true_model"] = {"body1_chi_SI": 0.06, "body1_cells": int(b1.sum()), "body2_chi_SI": 0.035, "body2_cells": int(b2.sum()), "background_chi_SI": 0.0} # ---------------------------------------------------------------- forward chi_map = maps.IdentityMap(nP=nC) sim = mag.Simulation3DIntegral(mesh, survey=survey, chiMap=chi_map, active_cells=active, store_sensitivities="ram", engine="geoana", numba_parallel=True) d_true = sim.dpred(chi_true) log["sensitivity_matrix"] = {"shape": [int(s) for s in sim.G.shape], "megabytes": round(float(sim.G.nbytes)/1e6, 1)} log["forward"] = {"all_finite": bool(np.all(np.isfinite(d_true))), "tmi_min_nT": float(d_true.min()), "tmi_max_nT": float(d_true.max()), "tmi_rms_nT": float(np.sqrt(np.mean(d_true**2)))} # ---------------------------------------------------------------- noise noise_floor = 2.0 # nT dobs = d_true + rng.normal(scale=noise_floor, size=d_true.size) std = np.full(dobs.size, noise_floor) mag_data = data.Data(survey, dobs=dobs, standard_deviation=std) target_phi_d = 1.0 * survey.nD # phi_d_star = nD (verified) log["uncertainty"] = {"noise_floor_nT": noise_floor, "relative_error": 0.0, "target_phi_d": float(target_phi_d), "target_convention": "chifact * nD (phi_d_star = nD)"} m0 = np.full(nC, 1e-4) def run(kind): dmis = data_misfit.L2DataMisfit(simulation=sim, data=mag_data) if kind == "L2": reg = regularization.WeightedLeastSquares( mesh, active_cells=active, reference_model=m0, alpha_s=1.0, alpha_x=1.0, alpha_y=1.0, alpha_z=1.0) dlist = [directives.UpdateSensitivityWeights(every_iteration=False), directives.BetaEstimate_ByEig(beta0_ratio=5.0, seed=7), directives.BetaSchedule(coolingFactor=2.0, coolingRate=1), directives.TargetMisfit(chifact=1.0), directives.SaveOutputEveryIteration(save_txt=True, name=f"mag_{kind}"), directives.UpdatePreconditioner()] opt = optimization.ProjectedGNCG(maxIter=25, lower=0.0, upper=0.15, maxIterLS=20, maxIterCG=30, tolCG=1e-3) else: reg = regularization.Sparse(mesh, active_cells=active, reference_model=m0, norms=[0.0, 2.0, 2.0, 2.0], gradient_type="total") dlist = [directives.UpdateSensitivityWeights(every_iteration=False), directives.BetaEstimate_ByEig(beta0_ratio=5.0, seed=7), directives.UpdateIRLS(cooling_factor=2.0, cooling_rate=1, chifact_start=1.0, chifact_target=1.0, max_irls_iterations=25, f_min_change=1e-3), directives.SaveOutputEveryIteration(save_txt=True, name=f"mag_{kind}"), directives.UpdatePreconditioner()] opt = optimization.ProjectedGNCG(maxIter=40, lower=0.0, upper=0.15, maxIterLS=20, maxIterCG=30, tolCG=1e-3) ip = inverse_problem.BaseInvProblem(dmis, reg, opt) inv = inversion.BaseInversion(ip, directiveList=dlist) t = time.time() mrec = inv.run(m0) dpred = sim.dpred(mrec) res = (dobs - dpred) / std rec = {"n_iterations": int(opt.iter), "final_phi_d": float(ip.phi_d), "final_phi_m": float(ip.phi_m), "final_beta": float(ip.beta), "achieved_chi": float(ip.phi_d / survey.nD), "reached_target": bool(ip.phi_d < target_phi_d), "residual_sd": float(res.std()), "chi_max_SI": float(mrec.max()), "chi_in_body1_median": float(np.median(mrec[b1])), "chi_in_body2_median": float(np.median(mrec[b2])), "chi_outside_median": float(np.median(mrec[~(b1 | b2)])), "frac_cells_above_1pct": float((mrec > 0.01).mean()), "wall_seconds": round(time.time() - t, 1)} np.save(f"model_{kind}.npy", mrec) np.save(f"dpred_{kind}.npy", dpred) print(f"[{kind}] chi={rec['achieved_chi']:.3f} reached={rec['reached_target']} " f"iters={rec['n_iterations']} {rec['wall_seconds']}s", flush=True) return rec log["inversion_L2"] = run("L2") log["inversion_sparse"] = run("sparse") np.save("chi_true.npy", chi_true); np.save("dobs.npy", dobs) np.save("d_true.npy", d_true); np.save("std.npy", std) np.save("cell_centers_active.npy", cc); np.save("active.npy", active) np.save("rx_loc.npy", rx_loc) np.save("topo_grid.npy", np.c_[TX.ravel(), TY.ravel(), topo[:,2]]) mesh.write_UBC("mesh.msh") if hasattr(mesh, "write_UBC") else None log["wall_seconds_total"] = round(time.time() - t0, 1) with open("mag3d.json", "w") as fh: json.dump(log, fh, indent=2) print(json.dumps(log, indent=2))