"""3-D magnetic inversion beneath REAL topography — Mount St Helens, WA. Real inputs ----------- * Topography: NOAA NCEI ETOPO 2022 v1, 15 arc-second surface elevation (tile ETOPO_2022_v1_15s_N60W135_surface.nc), subset and projected to a local tangent plane. 2318 m of genuine relief. * Inducing field: IGRF-14 evaluated at 46.20N, 122.19W, 2 km, 2024-06-01 -> |B| = 52136 nT, I = +67.81 deg, D = +14.80 deg. Synthetic component (deliberate, so a true model exists to validate against) --------------------------------------------------------------------------- A vertical conduit plus a dipping dyke, susceptibilities chosen in the range reported for arc volcanic rocks. The anomaly is forward-modelled from these, then noise is added. There is no measured magnetic survey here, so the recovered model is validated against a KNOWN target rather than asserted against nature. """ 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(11) P = np.load("msh_proj.npz") xg, yg, Zdem = P["x"], P["y"], P["z"] AMP, INC, DEC = float(P["amp"]), float(P["inc"]), float(P["dec"]) log = {"versions": {"simpeg": simpeg.__version__, "numpy": np.__version__}, "real_inputs": { "topography": "NOAA NCEI ETOPO 2022 v1, 15 arc-sec surface elevation, tile N60W135", "dem_cells": [int(Zdem.shape[0]), int(Zdem.shape[1])], "dem_elev_min_m": float(Zdem.min()), "dem_elev_max_m": float(Zdem.max()), "dem_relief_m": float(Zdem.max() - Zdem.min()), "inducing_field": "IGRF-14 @ 46.20N 122.19W, 2 km, 2024-06-01", "amplitude_nT": AMP, "inclination_deg": INC, "declination_deg": DEC, "location": "Mount St Helens, Washington, USA"}, "synthetic_component": "susceptibility bodies (conduit + dyke); the magnetic anomaly is forward-modelled, not measured"} # ---- real topography as an (N,3) point cloud ------------------------------- XG, YG = np.meshgrid(xg, yg) topo = np.c_[XG.ravel(), YG.ravel(), Zdem.ravel()] # ---- survey: terrain-draped airborne line grid, 120 m above ground --------- sx = np.arange(-12000, 12001, 1000.0) sy = np.arange(-14000, 14001, 1000.0) SX, SY = np.meshgrid(sx, sy) # drape: bilinear sample of the real DEM at each station from scipy.interpolate import RegularGridInterpolator dem_i = RegularGridInterpolator((yg, xg), Zdem, bounds_error=False, fill_value=float(Zdem.mean())) z_ground = dem_i(np.c_[SY.ravel(), SX.ravel()]) DRAPE = 120.0 rx_loc = np.c_[SX.ravel(), SY.ravel(), z_ground + DRAPE] rx = mag.Point(rx_loc, components="tmi") src = mag.UniformBackgroundField(receiver_list=[rx], amplitude=AMP, inclination=INC, declination=DEC) survey = mag.Survey(src) log["survey"] = {"nD": int(survey.nD), "line_spacing_m": 1000.0, "station_spacing_m": 1000.0, "drape_above_ground_m": DRAPE, "sensor_elev_min_m": float(rx_loc[:,2].min()), "sensor_elev_max_m": float(rx_loc[:,2].max())} # ---- octree mesh refined on the REAL surface and the survey ---------------- dh, nb = 250.0, 2**7 # 250 m core cells, 32 km per axis mesh = TreeMesh([[(dh, nb)], [(dh, nb)], [(dh, nb)]], origin="CCN") mesh.origin = np.r_[mesh.origin[0], mesh.origin[1], -mesh.h[2].sum() + float(Zdem.max())] mesh.refine_surface(topo, level=-1, padding_cells_by_level=[4, 4, 4], finalize=False) mesh.refine_points(rx_loc, level=-2, padding_cells_by_level=[4, 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), refined on real ETOPO surface", "n_cells": int(mesh.n_cells), "n_active": nC, "core_cell_size_m": dh, "extent_x_m": [float(mesh.nodes_x.min()), float(mesh.nodes_x.max())], "extent_z_m": [float(mesh.nodes_z.min()), float(mesh.nodes_z.max())]} # ---- true model: vertical conduit + dipping dyke --------------------------- chi_true = np.zeros(nC) X, Y, Zc = cc[:,0], cc[:,1], cc[:,2] z_surf_at = dem_i(np.c_[Y, X]) # real ground elevation per cell # conduit: near-vertical pipe under the edifice conduit = (np.hypot(X - 500.0, Y - 300.0) < 1400.0) & (Zc < z_surf_at - 400.0) & (Zc > -4500.0) chi_true[conduit] = 0.075 # dyke: dips ~35 deg to the NE, offset from the conduit dd = Zc - (-1200.0) + 0.70 * (X - 6500.0) dyke = (np.abs(dd) < 700.0) & (np.abs(Y - 4000.0) < 6500.0) \ & (Zc < z_surf_at - 500.0) & (Zc > -5000.0) chi_true[dyke] = 0.040 log["true_model"] = {"conduit_chi_SI": 0.075, "conduit_cells": int(conduit.sum()), "dyke_chi_SI": 0.040, "dyke_cells": int(dyke.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) # smoke-test the Jacobian before spending the run (skill: troubleshooting.md) _ = sim.getJ(np.full(nC, 1e-4)) log["jacobian_ok"] = 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 + uncertainties ------------------------------------------------- floor = 5.0 # nT, typical airborne mag noise dobs = d_true + rng.normal(scale=floor, size=d_true.size) std = np.full(dobs.size, 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": floor, "target_phi_d": float(target_phi_d), "convention": "chifact * 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) dl = [directives.UpdateSensitivityWeights(every_iteration=False), directives.BetaEstimate_ByEig(beta0_ratio=5.0, seed=11), # coolingFactor=2.0 halved beta every iteration and the run jumped # 1108 -> 314 in ONE step, overshooting the target of 725 and fitting # below the noise (residual sd 0.66). A gentler factor lets TargetMisfit # stop the run near chi=1 instead of well past it. directives.BetaSchedule(coolingFactor=1.3, coolingRate=1), directives.TargetMisfit(chifact=1.0), directives.SaveOutputEveryIteration(save_txt=True, name=f"msh_{kind}"), directives.UpdatePreconditioner()] opt = optimization.ProjectedGNCG(maxIter=60, lower=0.0, upper=0.2, 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") dl = [directives.UpdateSensitivityWeights(every_iteration=False), directives.BetaEstimate_ByEig(beta0_ratio=5.0, seed=11), directives.UpdateIRLS(cooling_factor=2.0, cooling_rate=1, chifact_start=1.0, chifact_target=1.0, max_irls_iterations=30, f_min_change=1e-3), directives.SaveOutputEveryIteration(save_txt=True, name=f"msh_{kind}"), directives.UpdatePreconditioner()] opt = optimization.ProjectedGNCG(maxIter=45, lower=0.0, upper=0.2, maxIterLS=20, maxIterCG=30, tolCG=1e-3) ip = inverse_problem.BaseInvProblem(dmis, reg, opt) inv = inversion.BaseInversion(ip, directiveList=dl) t = time.time(); mrec = inv.run(m0) dpred = sim.dpred(mrec); res = (dobs - dpred) / std body = conduit | dyke 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_conduit_median": float(np.median(mrec[conduit])), "chi_dyke_median": float(np.median(mrec[dyke])), "chi_background_median": float(np.median(mrec[~body])), "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"it={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_xyz.npy", topo) np.save("dem_grid.npy", Zdem); np.save("dem_x.npy", xg); np.save("dem_y.npy", yg) log["wall_seconds_total"] = round(time.time() - t0, 1) json.dump(log, open("msh3d.json", "w"), indent=2) print(json.dumps(log, indent=2))