What a metal detector hears, 10 Hz to 1 MHz

8. Detailed methodology

Everything the simulations do, written so the method can be reproduced or matched against the literature. Symbols are as in the table at the top of this tab; new ones are defined where they first appear.

8.1 Hardware model

• Grid. Two identical planar grids of n × n square cells of pitch p = L_b/n, where L_b is the board width (m). The RX grid is at height h above the ground and the TX grid at h + 1.6 mm. A pattern is a 0/1 vector a ∈ {0,1}^N, N = n², over cells. The switches route a single-turn current around the boundary of the selected cells. By superposition this equals the sum of the single-cell loops it contains, each carrying the same circulation.
• Configuration. An ordered pair (a, b) of TX and RX patterns. Families used: all aligned squares of size 16/8/4/2 cells (multi-scale, rank 64); all single-cell pairs; the 256 corner-anchored rectangles R(i, j) = {cells (k, l): k < i, l < j}, which have full rank. Reciprocity lets only the pairs a ≤ b be measured.
• Electrical. Loop resistance R = P_e·(R_seg + R_sw)·N_t, where P_e is the number of boundary edges, R_seg = ρ_Cup/(w t) the resistance of one trace segment (trace width w, thickness t), R_sw the on-resistance of one switch and N_t the number of turns. TX current I_a = min(I_max, √(P_max/R_a)).

8.1b Spectral stage (multi-frequency soil fit, explainer tab)

• The detector transmits N_f = 8 tones from 100 Hz to 100 kHz and measures in-phase and quadrature at each. The soil is fitted by least squares as a mix of three spectral shapes: magnetic (χ, flat in-phase), viscous VRM (in-phase ∝ −ln f, constant quadrature) and conductive ground (quadrature ∝ f).
• Target strength. The gold sphere response M(ω) (Wait 1951; conductivity σ_Au = 15 MS/m) is projected out of the span of those shapes. The surviving amplitude per tone, m_eff(a) = ‖residual‖/√N_f, is 0.055–0.09 depending on the radius a (m). The spatial stage uses the dipole strength q = 4πa³ m_eff(a) (m³).
• Clutter. The soil departs from the three shapes by a fraction ε, and varies from patch to patch by v = 10%. The residual susceptibility left after the fit therefore has standard deviation σ_χ = χ₀·v·ε, and a height change survives the fit in the same proportion, h_res = ε.
• The two stages are treated as separable: the spatial simulation runs at one effective frequency (10 kHz) with these residuals.

8.2 Fields and coupling

• Biot–Savart, finite straight segment from p₁ to p₂ carrying 1 A, at point r, with r₁ = r − p₁ and r₂ = r − p₂:
  B(r) = (μ₀/4π) · (|r₁| + |r₂|) · (r₁ × r₂) / ( |r₁||r₂| (|r₁||r₂| + r₁·r₂) )
  Cell field B_k(r) is the sum over its 4 edges; pattern field B_a(r) = N_t Σ_k a_k B_k(r).
• Mutual inductance between straight parallel filaments (Neumann, closed form; perpendicular filaments contribute zero). Segment 1 runs from 0 to ℓ₁ along an axis; segment 2 runs from s to s + ℓ₂ parallel to it, at perpendicular distance d_⊥. With f(x) = x·asinh(x/d_⊥) − √(x² + d_⊥²):
  M = ±(μ₀/4π) [ f(s+ℓ₂) − f(s+ℓ₂−ℓ₁) − f(s) + f(s−ℓ₁) ]
  Cell–cell matrix M_kl between the TX and RX layers; pattern coupling M_ab = N_t² a^T M b.
• Response to a small magnetic volume (first-order Born / reciprocity). A volume dV (m³) with susceptibility χ at r changes the RX voltage for configuration (a, b) by
  δV_ab = (ω I_a / μ₀) · χ dV · B_a(r)·B_b(r)
  The sensitivity kernel B_a·B_b is what the 3D view colours.
• Target. A dipole of strength q (m³) at r_t adds g_ab = (ω I_a/μ₀) q B_a(r_t)·B_b(r_t) to the reading. Each depth d has 7 candidates: the centre r_t = (0, 0, −d) and probes offset by ±0.15 d in x, y and depth.

8.3 Soil, clutter and errors

• Soil discretisation. Voxels in layers whose thickness grows with depth z (m) as (h + z)/κ, so they resolve the sensitivity kernel, which varies on the scale of the distance to the wires. Each layer's lateral size is the largest Δ/k (k a whole number) not above (h + z)/κ, where Δ is the spacing between survey spots; so every spot shift moves the voxels by whole cells, and deep layers are capped at Δ. κ = 2.5–3 for boards up to 30 cm, coarser for larger boards. The domain extends L_b/2 + L_b/3 + 15 cm sideways and 45 cm down.
• Clutter prior. The residual susceptibility is a zero-mean Gaussian field with covariance σ_χ² exp(−|r − r'|/ℓ). The voxel unknowns are the voxel averages, whose covariance Σ is computed exactly by 3 × 3 × 3 sub-point quadrature. Point values at voxel centres are not used: with a rough field and a weighted detector, that understates the clutter (Kaipio & Somersalo's "inverse crime").
• Uniform soil. S0_ab = (ω I_a)(χ₀/2) a^TM^imgb, the exact quasi-static response of a permeable half-space. M^img is the cell–cell mutual inductance to the RX grid's mirror image at distance 2h + 1.6 mm.
• Height error. One unknown per spot, prior standard deviation σ_h = 2 mm and correlation exp(−|Δx|/10 cm) between spots. Its effect on reading (a, b) is c_ab = h_res σ_h ∂S0_ab/∂h.
• Slow drift. Calibration error η S0_ab and coupling drift ζ ω I_a|M_ab| combine to an error e_ab (V). It is one unknown per configuration, the same at every spot. The coupling term uses the un-nulled coupling.
• Fast jitter (V²): (η_f S0_ab)² + (ζ_f ω I_a|M_ab|)², independent per reading and added to the white noise. It does not shrink with time.
• White noise per reading (V²): (4k_BT_KR_b + e_n²)/τ, where τ = 9T/(N_spots N_cfg) is the time per reading (s). In the ADC runs it also includes (ω I_a|M_ab|ν/SNR_ADC)²/(f_s/2)/τ, with ν the analog null factor, SNR_ADC = 10^{(6.02B+1.76)/20} for B effective bits, and f_s the sample rate (Hz).

8.4 Detector, search and figure of merit

• Measurement model: y = α g + G w + C θ + E ξ + n, where:
  • w ~ N(0, Σ) is the soil;
  • θ ~ N(0, H_corr) is the height;
  • ξ ~ N(0, I) is the slow drift, with E holding e_ab in every spot's row for configuration (a, b);
  • n ~ N(0, D) is the white noise.
  The optimal (prewhitened matched-filter) SNR is SNR² = g^TR⁻¹g, with R the covariance of everything but the nugget.
• Exact elimination of the drift. Each configuration's e and D are the same at every spot, so marginalising ξ is the same as centring that configuration's whitened readings across the N spots: r → r − (γ/N)Σ_spots r, where γ = 1 − (1 + N e²/D)^−1/2. The rest is the information form: Λ = [[Σ⁻¹ + G_c^TG_c, G_c^TC_c], [·, H_corr⁻¹ + C_c^TC_c]], where the subscript c marks the centred, whitened rows, and SNR² = g_c^Tg_c − u^TΛ⁻¹u. Systems with at most 6,000 readings are solved directly in measurement space instead.
• Search over position and depth. The filter output's correlation between each depth's centre and its probes gives the roughness Λ_s = diag(λ_x, λ_y, λ_z) (1/m²). With the Euler-characteristic densities ρ₃(t) = (t² − 1)e^−t²/2/(2π)² and ρ₂(t) = t e^−t²/2/(2π)^3/2, the expected false alarms per m² of a search over 5–60 cm are L₃ρ₃(t) + L₂ρ₂(t), where L₃ = ∫√(λ_xλ_yλ_z) dz and L₂ is the top and bottom face term (Worsley et al. 1996). t is set for 0.01 false alarms per m². A brute-force search over 9,261 candidates agreed within 10% (coil) and 1.5× in rate (cell pairs), under 0.1 in t.
• Minimum mass. The nugget is found with 50% probability when its SNR reaches t: q_min = t / SNR_{per unit q}, and m_min solves 4πa³ m_eff(a) = q_min.

8.5 Algorithm (pseudocode)

for each system S (conventional coil | grid family):
    configs (a, b), a ≤ b;  I_a = min(I_max, sqrt(P_max / R_a));  tau = 9 T / (N_spots · N_cfg)
    S0, dS0/dh, coupling M_ab   from exact image-loop / Neumann mutual inductances
    D = white noise + fast jitter;  e = slow drift;  gamma = 1 − (1 + N e²/D)^(−1/2)
    F = whitened rows of every config for the board at the origin, on the voxel grid extended by the spot range
    rows at spot p = F re-indexed (spot shifts are whole voxels):  G_p = F[:, idx_p]
    H = Σ_p S[idx_p, idx_p] − Aᵀ diag(γ(2 − γ)/N) A,   S = FᵀF,  A = Σ_p G_p          # GPU fp32, S in row panels
    Λ = [[Σ⁻¹ + H, h_c], [h_cᵀ, H_corr⁻¹ + C_cᵀC_c]];  factor Λ (GPU fp32, diagonally scaled)
    refine: r = exact centred residuals of all readings (fp64-accurate);  x ← x + Λ⁻¹(Gᵀr − Πx)
            until the objective and the frozen weights' soil variance have both converged;  M = J(x) − correction
    if refinement stalls (tiny noise floor): redo with H accumulated in fp64 on the CPU
    t = search threshold from M;  q_min = t / sqrt(M_centre);  m_min from m_eff(a)
    check: freeze w = R⁻¹g; soil variance of w on a continuum grid of spacing ≤ min(7.5 mm, p/3, h/2)
           K(x) = Σ_δ U_δ ⊛ Φ_δ  (one FFT convolution per TX→RX cell offset δ);  var = w^T(D + EE^T + CHC^T)w + soil(K)

8.6 Numerical checks

• Brute force. On small cases the solver matches a dense R built explicitly in fp64 to 2×10⁻⁴–6×10⁻⁴, the limit set by rounding the rows themselves to fp32.
• fp64 reference. For 16 × 16 cell pairs at unlimited time (the worst-conditioned case), the SNR agrees with a full fp64 accumulation to 4×10⁻⁶ at every depth. The previous fp32 code, without refinement, was off by up to 3×10⁻³ there.
• Soil resolution. Every detector is frozen and its soil variance recomputed on a regular continuum grid; the ratio of model to achieved SNR (the penalty) is reported with each result, along with the mass the frozen detector actually achieves. The grid spacing is at most a third of a cell: on a 15 cm board a 7.5 mm grid missed most of the effect (penalty 1.3 instead of 3.6 for cell pairs), while halving 3.1 mm to 1.6 mm changes corner rectangles by 3%. The check's convolution form matches a direct evaluation to 10⁻⁶.
• Penalty arithmetic. The achieved variance is the sum of the frozen weights' noise, drift, height and fine-grid soil terms, each computed directly. (Subtracting the model's soil term from the model's total, as before, cancels badly when time is unlimited.)
• Unlimited time. T = ∞ removes the electronic noise; the fast jitter keeps the covariance non-singular. Without it, full-rank switch sets improve without limit, because they fit the discretised soil exactly.

8.7 Terms for a literature search

The pieces of this method each have established names. The combination (a switch-reconfigurable coil matrix used as a multi-static electromagnetic induction array, with a joint inversion for near-surface soil clutter, for small-target detection) is what to search for.

• Hardware: reconfigurable / switchable coil matrix; multi-static or MIMO electromagnetic induction (EMI) array; magnetic induction tomography (MIT); coil arrays for UXO classification (e.g. TEMTADS, MetalMapper); induction-balance coil design; programmable bucking/nulling of TX–RX coupling.
• Physics: Born approximation for EMI with magnetic soil; sensitivity kernel / Fréchet derivative B_TX·B_RX; reciprocity theorem for coils; magnetic polarizability of a conducting sphere (Wait 1951); viscous remanent magnetisation (VRM) and frequency-dependent susceptibility of soils; soil compensation for metal detectors.
• Signal processing: generalised least squares; prewhitened matched filter in correlated Gaussian clutter; detection in clutter with a Gaussian-process prior; Woodbury identity; nuisance-parameter marginalisation; Hadamard / multiplex (Fellgett) advantage and its loss under signal-proportional noise; optimal experimental design for linear inverse problems.
• Analogies: MRI receive-coil arrays and parallel imaging (SENSE, g-factor); synthetic-aperture and MIMO radar; electrical impedance tomography.