Single‐Sensor Magnetometer Weight And Distance Estimation Script

This document describes a method for estimating the distance, depth, and approximate weight of magnetic targets using single-sensor magnetic data.

The method uses the 2D analytic signal of the magnetic field to locate anomaly peaks and evaluate their shape. From this, the script estimates the sensor-to-source distance, target depth below ground, and an approximate weight.

Why it is used

It is intended for rapid interpretation of marked magnetic anomalies when only one sensor is available. The method is most useful for UXO detection and buried-metal target surveys, but can be used for other geophysical applications where a practical automated estimate is needed for semi-small targets.

What the script outputs

For each detected anomaly, the script produces:

• an estimated distance,
• an estimated weight,
• and corresponding minimum, maximum, and central estimate values.

The following sections describe the input requirements, physical model, algorithm, and outputs in detail.

---

1. Prerequisites and Input Requirements

1.1 Prior Steps

• Anomalies must be marked in the data file before running this script (Mark = 1). More information about setting marks can be found here and here.
• The script exits with a message if no marked rows are found.

1.2 Input File Format

Standard MagNIMBUS CSV file with the following required columns:

Column	Default name	Description
Timestamp	Date + Time (or Timestamp)	Used to split flight lines
Magnetometer field	TMI	Total magnetic intensity (nT)
Altitude AGL	Altitude AGL	Height above ground level (metres)
Latitude / Longitude	Latitude / Longitude	Position (decimal degrees)
Mark	Mark	Anomaly mark (1 = marked row)

All column names are configurable via command-line arguments (see Section 6).

1.3 Background Removal

The script performs spatial lowpass filtering per flight line:

1. The TMI profile is resampled onto a uniform spatial grid (0.5 m step).
2. A 4th-order zero-phase Butterworth filter with cut-off wavelength LOWPASS_WAVELENGTH_M = 15 m (f_cut = 1/λ cycles/m) is applied to obtain the smooth background.
3. A rolling spatial-median pass (window = λ / 0.5 m = 30 samples) is applied after the LP filter for robustness against residual spikes.
4. The background is interpolated back to the original sample positions, stored in the TMI_LPF column, and subtracted from TMI to obtain the residual anomaly used by all subsequent steps.

If a flight line is shorter than λ/2 = 7.5 m or contains fewer than 4 samples, the algorithm falls back to using the line mean as the background estimate.

A legacy "linear" mode (polynomial trend fit to unmarked rows, TREND_DEGREE = 1) remains available by setting BACKGROUND_METHOD = "linear" in the script header.

No pre-processing of the background field is required before running the script.

---

2. Physical Model and Assumptions

2.1 Source Model

The magnetic source is modelled as a point magnetic dipole (compact object, dimensions much smaller than depth). This is a reasonable approximation for UXO and similar metallic objects at depths greater than ~2× the object's largest dimension.

2.2 The Analytical Signal — 2D Gradient Formulation

The script uses the 2D gradient analytic signal (Roest et al. 1992) computed on a regular 2D grid:

$$ A(x, y) = \sqrt{\left(\frac{\partial B}{\partial x}\right)^2 + \left(\frac{\partial B}{\partial y}\right)^2 + \left(\frac{\partial B}{\partial z}\right)^2} $$

where:

• $B(x, y)$ is the detrended magnetic field gridded onto a regular mesh,
• $\partial B/\partial x$, $\partial B/\partial y$ are horizontal gradients computed via a 5-point finite-difference stencil,
• $\partial B/\partial z$ is the vertical derivative computed via 2D FFT (multiplication by the wavenumber $|k| = \sqrt{k_x^2 + k_y^2}$ in the frequency domain).

Physical meaning: $A(x, y)$ is amplitude-independent: its peak lies directly above the source and its width encodes the source depth regardless of magnetisation direction (Roest et al. 1992).

2.3 Two Depth Estimation Methods

Both methods are computed for every anomaly and their results are combined by the selection policy in Section 3, Step 12.

Method A — Half-Width of the Analytic Signal Peak

The full width at half maximum (FWHM) of the local $A$ peak near the mark, sampled along the along-track profile:

$$ W = x_{\mathrm{right}} - x_{\mathrm{left}} $$

where $x_{\mathrm{left}}$, $x_{\mathrm{right}}$ are the nearest crossing points on each side of the AS peak where $A = A_{\max}/2$, found by linear interpolation.

Before computing $W$, an asymmetry cap is applied: if one half-width exceeds the other by more than ASYMMETRY_CAP_RATIO = 1.20, the longer side is trimmed to 1.20 × shorter side. This corrects geological gradients that inflate one side of the AS profile.

$$ d_A = K \cdot W, \quad K = 1.03 $$

$d_A$ is discarded if $d_A >$ MAX_DISTANCE_M = 20 m.

Method B — B/AS Ratio at the AS Peak

Both $B$ and $A$ are evaluated at the AS-peak position (the bin with highest $A$ near the mark centre):

$$ d_B = C_B \cdot \frac{|B(x_{AS_peak})|}{A(x_{AS_peak})}, \quad C_B = 3.0 $$

Physical meaning: for a point dipole, the ratio $|B|/A$ at the anomaly centre is a fixed multiple of the sensor-to-source distance. Evaluating both quantities at the same AS-peak position ensures self-consistency: any systematic AS inflation (e.g. from geology) is largely cancelled in the ratio.

$d_B$ is discarded if $A <$ MIN_AS_THRESHOLD = 1.0 nT/m, $|B| < 1$ nT, or $d_B >$ MAX_DISTANCE_M.

2.4 Theoretical Derivation of Constants

For a dipole of moment $m$ at sensor-to-source distance $D$:

• Peak TMI: $|B_{\max}| = \frac{\mu_0}{4\pi} \frac{2m}{D^3}$
• Peak AS: $A_{\max} = \frac{\mu_0}{4\pi} \frac{6m}{D^4} = \frac{3|B_{\max}|}{D}$

From these:

$$ D = 3 \cdot \frac{|B_{\max}|}{A_{\max}} \quad \Rightarrow \quad C_B = 3.0 $$

The theoretical FWHM of $A(x)$ along a profile through a vertical dipole is $W \approx D$, giving $K_{\mathrm{theoretical}} = 1.0$. The calibrated value $K = 1.03$ is the result of empirical fitting to the Brandenburg 2024 ground-truth dataset (5 targets, depths 0.3–2.5 m), where the 2D grid AS implementation consistently produced slightly wider peaks than the theoretical 1D profile.

2.5 Depth Below Ground Surface

$$ d_{\mathrm{ground}} = d_{\mathrm{sensor}} - \bar{h}_{\mathrm{AGL}} $$

where $\bar{h}{\mathrm{AGL}}$ is the mean altitude AGL over the aggregated bins of the anomaly group. If $d{\mathrm{ground}} < 0$ the value is clipped to 0.

2.6 Weight Estimation Model

Weight estimation runs immediately after depth estimation for each anomaly group (Steps 14–18). It reuses the same data, grid, and anomaly groups already computed during the depth-estimation run.

2.6.1 Peak Field Value for Weight ($b_{\max}$)

The field value used for weight estimation is the residual TMI at the AS-peak position:

$$ |B_{\max}| = |B(x_{AS_peak})| $$

If this value is unavailable (Method B was skipped due to peak_mismatch), the fallback is:

$$ |B_{\max}| = \max_{|x - x_{\mathrm{mark}}| \leq 5,\mathrm{m}} |B_{\mathrm{anom}}(x)| $$

covering all bins within ±AS_PEAK_SEARCH_RADIUS_M = 5 m of the mark centre. This handles cases where the operator pressed the mark button before crossing the anomaly peak.

2.6.2 Dipole Moment from Peak Field and Distance

$$ \boxed{m\ [\mathrm{A·m}^2] = 5 \times 10^{-3} \cdot |B_{\max}|\ [\mathrm{nT}] \cdot D_{\mathrm{weight}}^3\ [\mathrm{m}^3]} $$

where $D_{\mathrm{weight}}$ is the weight distance (Step 14).

Important: The formula uses the sensor-to-source distance $D_\mathrm{weight}$, not the depth below ground. Using ground depth instead would underestimate weight by up to $(\bar{h}_\mathrm{AGL} / D)^3$.

2.6.3 Effective Magnetisation, Volume and Weight

$$ \boxed{M_{\mathrm{source}}\ [\mathrm{kg}] = K_M \cdot |B_{\max}|\ [\mathrm{nT}] \cdot D_{\mathrm{weight}}^3\ [\mathrm{m}^3]} $$

where $K_M = 5 \times 10^{-3} \times \rho / J_{\mathrm{eff}} \approx 0.01847\ \mathrm{kg/(nT·m}^3\mathrm{)}$ with calibrated constants.

Sphere geometry (only supported geometry): the calibrated effective magnetisation is $J_{\mathrm{eff}} = 2{,}112\ \mathrm{A/m}$. — geometric-mean fit over 5 steel objects with known weights (9–26 kg, Brandenburg 2024). Approximately 17.7× the theoretical induced-only value (119 A/m), accounting for remanent magnetisation.

2.6.4 Weight Range (Min / Max)

Min and Max weight estimates are computed by offsetting the central weight distance by $\pm f$, where $f =$ WEIGHT_SPREAD_D_FRAC $= 0.044$:

$$ D_{\min} = D_{\mathrm{weight}} \cdot (1 - f), \quad D_{\max} = D_{\mathrm{weight}} \cdot (1 + f) $$

Because $W \propto D^3$:

$$ \frac{W_{\max}}{W_{\min}} = \left(\frac{1+f}{1-f}\right)^3 \approx 1.30 \quad \Rightarrow \quad \mathrm{weight range} \approx 30,% $$

Using dist_min / dist_max directly for the weight range would amplify a wide distance band by $D^3$, producing a thousands-of-percent weight spread when the two methods disagree. Anchoring on Estimated_Distance_Harmonic with a fixed ±4.4 % offset gives a stable, meaningful weight bracket regardless of method disagreement.

---

3. Algorithm

Step 1 — Load and validate

• Read input CSV.
• Verify required columns (TMI or configured name, Latitude, Longitude, Mark).
• Strip previously computed output columns (_SCRIPT_OUTPUT_COLS) to prevent stale values from a prior run.
• Exit cleanly if no marks are found.

Step 2 — Split into flight lines

Consecutive rows separated by more than LINE_GAP_SECONDS = 30 s are treated as separate flight lines. Each line is processed independently.

Step 3 — Background removal per flight line

Default path (BACKGROUND_METHOD = "lowpass"):

1. Resample TMI profile to a uniform spatial grid (0.5 m step).
2. Apply 4th-order zero-phase Butterworth lowpass filter with LOWPASS_WAVELENGTH_M = 15 m.
3. Apply rolling spatial-median pass (window = 30 samples) for spike robustness.
4. Interpolate background back to original sample positions; store as TMI_LPF; subtract from TMI.

Degenerate cases (line too short or too few samples) subtract the line mean instead.

Legacy "linear" path: fit a degree-TREND_DEGREE polynomial to unmarked rows; subtract.

Step 4 — Coordinate projection

Convert latitude/longitude to local Cartesian (equirectangular, centred on data midpoint):

x_m = R · (lon − lon₀) · cos(lat₀)   [East, metres]
y_m = R · (lat − lat₀)                [North, metres]

Along-track distance is computed via the haversine formula for use in aggregation and halfwidth calculations.

Step 5 — Estimate grid cell size

If DEFAULT_CELL_SIZE_M is None (default):

1. Per-line median GPS step: compute all consecutive point-to-point distances per line (ignoring steps < 0.01 m); take the median across lines.
2. Validity check: if outside [0.5 m, 10 m], fall back to extent / 168.
3. Snap to 0.1 m: max(0.5, round(cs, 1)) — prevents grid-alignment aliasing from sub-0.1 m cell-size variations.

Step 6 — Build 2D grid

1. Grid dimensions: gridWidth = max(2, int((x_max − x_min) / cell_size)), gridHeight = max(2, int((y_max − y_min) / cell_size)).
2. Bin and median: multiple observations in the same cell are reduced to their median.
3. Fill gaps: empty cells within blanking_distance (default 10 × cell_size) are filled via thin-plate-spline interpolation (RBFInterpolator(kernel='thin_plate_spline')) from the ≥ 10 nearest known neighbours. Regularisation smoothing = var(vals) × 10⁻³. Falls back to linear griddata on solver failure or remaining NaNs.
4. Clip: interpolated cells are clamped to ±100 % of the measured data range to prevent TPS overshoot at grid edges.
5. Circular blank mask: cells further than blanking_distance from any data cell are set to NaN.

Step 7 — Compute 2D Analytic Signal

1. Fill remaining NaN with an inverse-distance-weighted priority-queue algorithm (exact reproduction of Java GridInterpolator.interpolate()).
2. Horizontal gradients $\partial B/\partial x$, $\partial B/\partial y$: 5-point finite-difference stencil (falls back to 3-point central, then forward/backward at edges).
3. Vertical derivative $\partial B/\partial z$: 2D FFT multiplication by $|k|$:

   G = fft2(B)           # NaN → 0 before FFT
   G *= sqrt(kx² + ky²)
   dz = real(ifft2(G))
   

4. AS grid: $A = \sqrt{(\partial B/\partial x)^2 + (\partial B/\partial y)^2 + (\partial B/\partial z)^2}$.
5. Cells that were NaN before fill are set to NaN in the AS output.

Step 8 — Sample AS at data positions

Bilinearly interpolate the AS grid at each original data point's $(x_m, y_m)$. NaN values are preserved (not converted to zero) to avoid diluting real AS amplitudes at mark locations.

Step 9 — Spatial aggregation

Bin each flight line into spatial bins of SPATIAL_BIN_M = 0.5 m along the along-track distance. Per bin:

• TMI and AS: mean (NaN excluded).
• AGL: mean of valid values.
• Mark: OR (a bin is marked if any row in it is marked).

The same bin keys are reused to aggregate the AS values, ensuring identical spatial grouping.

Step 10 — Identify anomaly groups

Contiguous runs of marked bins are identified. Each run is one anomaly group.

Step 11 — Extract window and compute methods

For each anomaly group at aggregated positions $[x_s, x_e]$:

• Window: $[x_s - 20,\mathrm{m},\ x_e + 20,\mathrm{m}]$
• Require ≥ MIN_WINDOW_SAMPLES = 20 bins; otherwise skip.
• Edge trim: exclude the outermost 10 % of bins from each end (edge = max(1, int(n × 0.10))) to remove FFT wrap-around artefacts.

Local AS peak search:

1. Restrict to bins within ±ANOMALY_CENTER_RADIUS_M = 15 m of the mark centre. This larger radius handles cases where the operator placed the mark several metres before crossing the actual anomaly centre.
2. Find the bin with the highest $A$ in that sub-range.
3. Fall back to the global window maximum if no bins fall within the search radius.

Method A — Halfwidth:

1. Search outward from the AS peak for the first crossing where $A \leq A_{\max}/2$ (linear interpolation).
2. Apply asymmetry cap: if one half exceeds the other by more than ASYMMETRY_CAP_RATIO = 1.20, trim the longer to 1.20 × shorter.
3. $W = x_{\mathrm{right}} - x_{\mathrm{left}}$, $\quad d_A = K \cdot W$ with $K = 1.03$.

$d_A$ is discarded if $W &lt;$ MIN_FWHM_M = 0.5 m, $W &gt;$ MAX_FWHM_M = 100 m, or $d_A &gt;$ MAX_DISTANCE_M = 20 m.

Method B — B/AS ratio:

Both field value and AS are taken at the AS-peak bin:

$$ d_B = C_B \cdot \frac{|B(x_{AS_peak})|}{A(x_{AS_peak})} $$

$d_B$ is discarded if $A &lt;$ MIN_AS_THRESHOLD = 1.0 nT/m, $|B| &lt; 1$ nT, or $d_B &gt;$ MAX_DISTANCE_M.

Dipole detection: Within ±DIPOLE_CHECK_RADIUS_M = 15 m of the mark centre, flag dipole_anomaly if $B_{\min} &lt; -0.15 \cdot B_{\max}$.

Step 12 — Best-distance selection

Estimated_Distance_Min = $\min(d_A, d_B)$, Estimated_Distance_Max = $\max(d_A, d_B)$ — always the raw method outputs.

Estimated_Distance_Harmonic (central best estimate) is selected by the following rule-based policy, applied in order:

Condition	Rule	Action
$d_B < \bar{h}_{\mathrm{AGL}}$ — AS-ratio result is sub-AGL (physically impossible)	Rule 1	Use $d_A$ (FWHM) alone
$d_A >$ LARGE_FWHM_RATIO $\times, d_B$ — FWHM greatly exceeds AS-ratio, indicating window-inflated halfwidth	Rule 2	Use $d_B$ (AS-ratio) alone
$d_B < d_A$ — AS-ratio underestimates, consistent with 2D-AS inflation by geology or grid noise	Rule 3	Use $d_A$ (FWHM) alone
Both methods agree and are physically valid	Fallback	Harmonic mean: $D_H = 2 d_A d_B / (d_A + d_B)$
Only one method produced a valid result	—	Use that result directly

The harmonic mean pulls the central estimate toward the smaller value compared with the arithmetic mean. When both methods agree ($d_A \approx d_B$), $D_H = D_{\mathrm{AM}}$ (no effect).

Estimated_Depth_Harmonic = $\max(0,\ D_H - \bar{h}_\mathrm{AGL})$; similarly for Min and Max.

Step 13 — Write depth outputs

Fill out_dist_min / max / avg and out_depth_min / max / avg arrays for all marked rows of the group. Proceed to Step 14.

Step 14 — Compute peak field ($b_{\max}$) and weight distance

$b_{\max}$ (peak residual field for weight):

Primary: the residual TMI at the AS-peak position (already found in Step 11):

$$ |B_{\max}| = |B(x_{AS_peak})| $$

Fallback (when the primary value is unavailable due to peak_mismatch): maximum absolute value within ±AS_PEAK_SEARCH_RADIUS_M = 5 m of the mark centre, or the marked bins if no window bins fall within 5 m.

Weight distance ($D_{\mathrm{weight}}$):

When Rule 3 fired (both distances valid, $d_B$ physically valid, $d_B &lt; d_A$, Rule 2 did not fire), use $d_B$ — the B/AS ratio is self-consistent at the AS-peak and cancels the geological inflation that caused the disagreement:

$$ D_{\mathrm{weight}} = d_B \quad \mathrm{(Rule 3 case)} $$

Otherwise:

$$ D_{\mathrm{weight}} = \max!\left(D_H,\ \bar{h}_\mathrm{AGL}\right) $$

The AGL floor ensures a physically valid distance when the depth estimator returns 0.

Step 15 — Compute effective magnetisation

Sphere geometry (only supported geometry): the calibrated effective magnetisation is $J_{\mathrm{eff}} = 2{,}112\ \mathrm{A/m}$ (calibrated constant).

Step 16 — Compute dipole moment

$$ m = K_{\mathrm{DIPOLE}} \cdot |B_{\max}|\ [\mathrm{nT}] \cdot D_{\mathrm{weight}}^3\ [\mathrm{m}^3] $$

Step 17 — Compute volume and weight estimates

$$ W_{\mathrm{harm}} = \frac{\rho \cdot m}{J_{\mathrm{eff}}} $$

$$ W_{\mathrm{min}} = \frac{\rho \cdot K_{\mathrm{DIPOLE}} \cdot |B_{\max}| \cdot \left(D_{\mathrm{weight}} \cdot (1-f)\right)^3}{J_{\mathrm{eff}}}, \qquad W_{\mathrm{max}} = \frac{\rho \cdot K_{\mathrm{DIPOLE}} \cdot |B_{\max}| \cdot \left(D_{\mathrm{weight}} \cdot (1+f)\right)^3}{J_{\mathrm{eff}}} $$

where $f$ is the parameter WEIGHT_SPREAD_D_FRAC, equal to $0.044$.

Step 18 — Write outputs

Append all result columns to marked rows of the group. Write source CSV and targets CSV. See Section 4.

---

4. Outputs

4.1 Additional Columns Added to Source CSV

Column	Description	Unit
TMI_LPF	Lowpass-filtered background subtracted in Step 3	nT
Estimated_Distance_Min	Smaller of $d_A$ and $d_B$	m
Estimated_Distance_Max	Larger of $d_A$ and $d_B$	m
Estimated_Distance_Harmonic	Best-estimate distance — rule-based (Step 12)	m
Estimated_Depth_Min	Depth fromDistance_Min, clipped to 0	m
Estimated_Depth_Max	Depth fromDistance_Max, clipped to 0	m
Estimated_Depth_Harmonic	Depth fromDistance_Harmonic, clipped to 0	m
Estimated_Weight_Min	Lower weight estimate from $D_{\mathrm{weight}} \cdot (1-f)$	kg
Estimated_Weight_Max	Upper weight estimate from $D_{\mathrm{weight}} \cdot (1+f)$	kg
Estimated_Weight_Harmonic	Central weight estimate from $D_{\mathrm{weight}}$	kg

All columns are empty ("") for rows outside any anomaly window.

4.2 Targets CSV ({stem}-targets-as.csv)

One row per marked row. Identical estimation values are written to all marked rows of the same anomaly group.

Contains all original source CSV columns (excluding previously computed output columns) plus all 10 output columns listed in Section 4.1 (TMI_LPF plus the 9 distance/depth/weight estimation columns).

---

5. Internal Constants

Not exposed as CLI arguments; edit the script header to change them.

Constant	Value	Description
K_FACTOR	1.03	FWHM → sensor-to-source distance (Method A); calibrated on Brandenburg 2024 dataset
CB_FACTOR	3.0	B/AS ratio → sensor-to-source distance (Method B); exact theoretical value
WINDOW_PADDING_M	20.0	Window extends this far beyond each mark group (m)
SPATIAL_BIN_M	0.5	Spatial bin size for along-track aggregation (m)
DEFAULT_CELL_SIZE_M	None	Override grid cell size (None = auto-estimated)
MIN_WINDOW_SAMPLES	20	Minimum bins in window for processing
ANOMALY_CENTER_RADIUS_M	15.0	Radius around mark centre for AS peak search (m); large to handle offset marks
AS_PEAK_SEARCH_RADIUS_M	5.0	Radius around mark centre for $b_{\max}$ fallback search (m)
DIPOLE_CHECK_RADIUS_M	15.0	Radius for dipole sign-change detection (m)
METHOD_DISAGREE_THRESHOLD	0.35	Relative distance difference to triggermethod_disagreement flag
LARGE_FWHM_RATIO	2.5	$d_A / d_B$ threshold above which Rule 2 selects $d_B$
ASYMMETRY_CAP_RATIO	1.20	Max allowed ratio between left and right AS half-widths before trimming the longer side
MIN_FWHM_M	0.5	Minimum valid full FWHM (m)
MAX_FWHM_M	100.0	Maximum valid full FWHM (m);wide_anomaly flag triggers when $W > 100\ \mathrm{m}$
MIN_AS_THRESHOLD	1.0	Minimum AS at AS-peak for Method B to be valid (nT/m)
MAX_DISTANCE_M	20.0	Maximum plausible sensor-to-source distance; results above this are discarded (m)
BACKGROUND_METHOD	"lowpass"	Background removal:"lowpass" (Butterworth LP) or "linear" (polynomial fit)
LOWPASS_WAVELENGTH_M	15.0	Spatial cut-off wavelength for lowpass background filter (m)
TREND_DEGREE	1	Polynomial degree for linear detrending (legacy mode)
MIN_TREND_SAMPLES	10	Minimum rows needed for linear background removal (legacy mode)
LINE_GAP_SECONDS	30	Time gap that splits data into separate flight lines (s)
CSV_SEPARATOR	","	CSV field delimiter
MARK_COLUMN	"Mark"	Hard-coded mark column name
K_DIPOLE	5e-3	Dipole moment coefficient A·m²·nT⁻¹·m⁻³ (exact: $\frac{4\pi}{2\mu_0}\times10^{-9}$)
J_EFF\_SPHERE\_AM	2112.0	Calibrated effective magnetisation, sphere geometry (A/m); ≈ 17.7× induced-only value
WEIGHT_SPREAD_D_FRAC	0.044	Distance offset for weight Min/Max: ±4.4 % in $D$ → ≈ 30 % weight range
DENSITY_STEEL	7800.0	Material density (kg/m³)

---

6. Command-Line Parameters

Argument	Default	Description
file_path	(required)	Input CSV file path
--mag-column	TMI	Magnetometer column name
--altitude-agl-column	Altitude AGL	Altitude AGL column (used for depth = distance − AGL, and for weight distance floor)
--output-dir	(same as input)	Directory for targets CSV