Physics-Guided Gaussian Process Mapping of Strong-Gradient Radiation Fields from Mobile Robot Surveys: The Role of Sampling Geometry

1. Introduction

Characterizing dose-rate distributions inside nuclear facilities underpins routine operation, decommissioning, and emergency response [1,2]. Mobile robots are increasingly deployed to acquire such measurements in place of human surveyors, following the as-low-as-reasonably-achievable (ALARA) principle [3,4,5]. The surveys conducted after the Fukushima Daiichi accident remain the canonical example [1]. A robotic survey, however, does not produce a dose map directly: it produces sparse, noisy point measurements distributed irregularly along the vehicle trajectory, from which a continuous field must be reconstructed [6,7].

The reconstruction problem is hardest precisely where it matters most. Radiation fields near collimated openings, shielding gaps, and maze entrances exhibit abrupt, strongly anisotropic transitions in which dose rate changes by orders of magnitude over tens of centimeters. Misplacing such a boundary in either direction carries asymmetric costs: an outward error wastes access time, while an inward error exposes workers. A reconstruction that draws a confident but wrong boundary is therefore more dangerous than one that draws an honest, uncertainty-qualified boundary. As mapping outputs are increasingly consumed by downstream planners rather than human experts [8], the ability of a method to report where it should not be trusted becomes as important as its point accuracy.

Two methodological traditions dominate robotic radiation mapping. The geostatistical tradition, comprising kriging and Gaussian process regression (GPR), handles irregular, noisy, clustered data naturally and supplies predictive variance [9,10,11]. West et al. demonstrated GPR mapping of an operating reactor with a mobile robot [9]. Jung et al. adapted the inverse-square law into a GPR kernel for three-dimensional mapping [10]. Zhang et al. introduced multi-kernel weighting to accommodate obstructions [12]. Poisson kriging addresses counting statistics explicitly [13]. These methods, however, typically employ isotropic, stationary priors whose smoothness assumptions blur exactly the abrupt boundaries that strong-gradient scenarios present. The machine-learning tradition, including backpropagation networks [14], U-Net variants [15], and coordinate-input multilayer perceptrons (MLPs), imposes no smoothness prior and can, in principle, represent sharp transitions, but provides no uncertainty output, and its dependence on the quantity and spatial distribution of training samples is rarely examined. Reported results are seldom accompanied by initialization-seed statistics.

A further gap concerns evaluation. Robotic survey data are spatially autocorrelated along trajectories, so the ubiquitous random train/test split largely measures interpolation between immediately adjacent points and can substantially overstate performance in unvisited regions. Few studies of robotic radiation mapping have systematically compared reconstruction methods under both a random hold-out protocol and a spatially blocked protocol, and the calibration of GP uncertainty under the latter has rarely been assessed.

This paper addresses these gaps with three contributions. First, using seven independent robotic surveys of a collimated Cs-137 field and Poisson-noise simulations, we show that the ranking of reconstruction methods is determined by sampling geometry: under uniform sampling the MLP is robustly superior at recovering abrupt boundaries, whereas under trajectory sampling its advantage disappears at visited locations and reverses severely in unvisited regions. Second, we propose a physics-guided GP for strong-gradient fields, in which a parametric collimated-beam template is fitted as the GP mean function and a Matérn 3/2 process models the residual, and show it achieves the highest extrapolation $R^2$ among eight baselines (median $0.61$ versus at most $0.31$), while independently recovering source position, beam direction, and cone half-angle consistent with the actual apparatus across all seven runs. A runtime fit-quality gate additionally makes the correctness of this prior an observable quantity, so the method can flag when its own assumptions break down. Third, we report the calibration of the GP predictive uncertainty under spatial block cross-validation, demonstrating that the resulting maps convey usable confidence information for survey planning and ALARA decision support. In relation to the GPR line of work [9,10,12], our approach extends the physics-informed direction from isotropic inverse-square priors to anisotropic, collimation-aware priors; in relation to neural approaches [14,15], our results delimit the sampling conditions under which they can be recommended.

2. Materials and Methods

2.1. Robotic Survey Platform and Data Acquisition

The survey platform (Figure 1) is a tracked mobile robot composed of an onboard computing unit, a tracked chassis, a light detection and ranging (LiDAR) sensor, a $\gamma$ dose-rate detector, and a lithium battery power system. The Robot Operating System (ROS) [16] serves as the software framework, providing unified management and communication scheduling for all sensors and actuators. The platform is teleoperated over Wi-Fi and sustains more than 2 h of continuous operation, sufficient for indoor facility inspection tasks.

The computing unit is an NVIDIA Jetson Nano edge module with a 128-core GPU and a quad-core ARM Cortex-A57 CPU, which provides adequate processing capacity at low power for sensor fusion, simultaneous localization and mapping (SLAM), and real-time processing of radiation data. A Slamtec RPLIDAR sensor provides 360° range measurements for two-dimensional environment mapping. The inertial measurement unit (IMU) is a GY-52 module integrating a three-axis gyroscope, a three-axis accelerometer, and a three-axis magnetometer for attitude and motion estimation. The $\gamma$ dose-rate detector, which reports the ambient dose equivalent rate ${\dot{H}}^{\ast }\left(10\right)$ (referred to as “dose rate” hereafter), is mounted vertically on top of the platform. It weighs less than 200 g, covers 0.1 µSv/h–10 Sv/h with an energy response of 60 keV–3 MeV and a response time below 3 s, and streams readings to the computing unit over USB [17]. Two-dimensional maps and poses are produced by the gmapping SLAM algorithm [18] fusing LiDAR, IMU, and wheel-odometry data. Dose-rate readings are logged at approximately 1 Hz and paired with the nearest-timestamp SLAM pose.

Localization accuracy was evaluated against an eight-camera motion-capture system (3D accuracy ±0.15 mm) as ground truth. After SE(3) alignment by the Umeyama method [19], the absolute trajectory error (ATE) over five repeated runs had a root-mean-square error (RMSE) of $2.72\pm 0.49$ cm (Figure 2), adequate for dose mapping at the decimeter scale. The median survey speed across the seven radiation runs was 0.14 m/s (max ≤0.3 m/s). Combined with the detector response time, the spatial smearing of an individual reading is approximately 0.4 m, below the characteristic structure scale of the fields considered here.

Radiation experiments were conducted in a laboratory of a national secondary calibration facility. A collimated Cs-137 source (Figure 3; lead collimator, fan-shaped aperture with nominal half-angle of 12°; activity $7.4\times {10}^8$ Bq, i.e., $0.74$ GBq) was oriented toward a wall more than 3 m away. In each of seven independent runs, the operator teleoperated the robot around the room perimeter to build the map, then swept repeatedly across the beam region, yielding 521–2325 georeferenced dose-rate measurements per run after removing pre-convergence poses. All dose rates d (in µSv/h) were transformed as $y={log}_{10}(d+1)$, where the additive constant of 1 µSv/h keeps the transform finite at zero dose and sets the background floor. Readings below 1 µSv/h are thereby mapped to $y\approx 0$ and treated as background. The transform compresses the several-decade dynamic range of the collimated field into a numerically stable target while preserving the ordering of dose rates.

2.2. Physics-Guided Gaussian Process

We model the log-domain field as $y\left(\mathbf{x}\right)=m\left(\mathbf{x}\right)+f\left(\mathbf{x}\right)+\epsilon$, where $m\left(\mathbf{x}\right)$ is a parametric physics template, $f\left(\mathbf{x}\right)\sim \mathcal{GP}(0,k)$ is a residual Gaussian process with a Matérn 3/2 kernel plus white noise, and $\epsilon$ is observation noise. The template represents a single dominant collimated source:

$$m\left(\mathbf{x}\right)={log}_{10}\left({\displaystyle \frac{Cg\left(\theta \right(\mathbf{x}\left)\right)}{r^2\left(\mathbf{x}\right)}}+b+1\right),$$

(1)

where r is the distance to the source position $(s_x,s_y)$, softened as $r^2\to r^2+r_c^2$ with $r_c=0.3$ m to remove the inverse-square singularity in the near field, where the point-source assumption does not hold and measurements can hardly be taken; $\theta$ is the angular offset from the beam direction $\phi$; and

$$g\left(\theta \right)=\ell +(1-\ell )\sigma \left({\displaystyle \frac{{\theta }_{\mathrm{half}}-\theta }{w}}\right)$$

(2)

is a sigmoid-smoothed sector profile with half-angle ${\theta }_{\mathrm{half}}$, transition width w, and leakage factor ℓ; b is the background level. The sigmoid replaces the ideal step-like sector edge with a differentiable transition, reflecting the finite penumbra of a physical collimator while keeping the template amenable to gradient-based least squares. The eight parameters are estimated from training data by robust nonlinear least squares (soft-L1 loss) with multi-start initialization: the source position is initialized at the maximum reading, and the direction at the intensity-weighted bearing of the hottest 15% of points, with additional starts rotated by $\pm {90}^{\circ }$ and ${180}^{\circ }$ to escape directional local minima. GP hyperparameters (signal variance, length scale, noise variance) are estimated by maximizing the marginal likelihood [20] on a subsample of at most 800 points and then held fixed while conditioning on all training points.

The construction itself is a collimation-aware instantiation of a classical idea—GP regression with a parametric (non-zero) mean, known in geostatistics as universal or external-drift kriging [20]. The contribution here lies not in the construction but in the anisotropic template family, the dual-protocol evidence for when such a prior is necessary, and the misspecification analysis of Section 4.

Two design choices merit emphasis. First, the division of labor is deliberate: the template mean governs behavior far from data, where a zero-mean GP reverts to a constant and a coordinate-input neural network (such as the MLP baseline of Section 2.3) extrapolates without constraint, while the residual GP corrects local deviations near data and supplies predictive variance. Second, the method has a defined scope and a runtime adequacy check: it targets scenarios with a single dominant source of known family. Template adequacy is monitored by a quantitative fit-quality gate whose statistic is the normalized robust fit cost $\overline{c}$ (the soft-L1 objective of the template fit divided by the number of training points). Whenever $\overline{c}$ exceeds a threshold $\tau$ calibrated on correctly-specified reference fits ($\tau =1.3\times$ the largest $\overline{c}$ observed for correct templates on held-out surveys), the template is flagged as inadequate—a runtime alarm prompting supplementary measurement or manual review, with a conservative fallback to the plain GP (setting $m\left(\mathbf{x}\right)$ to a constant) as the default. Section 4 quantifies the gate’s detection and false-alarm rates. Multi-source and unknown-source-type scenarios are future work.

2.3. Baseline Methods

Eight baselines span the main traditions: (i) linear interpolation on the Delaunay triangulation with nearest-neighbor fill; (ii) plain GPR with Matérn 3/2 kernel plus constant (bias) kernel and white noise, hyperparameters by marginal likelihood—equivalent to ordinary kriging, i.e. a GP with unknown constant mean [20] (the configuration closest to prior robotic mapping work [9]); (iii) inverse-distance weighting ($k=8$, $p=2$); (iv) random forest regression (200 trees); (v) a coordinate-input MLP (64–32–16, tanh activation, limited-memory BFGS optimizer [21], L2 regularization $\alpha =0.01$, standardized inputs), reported both as a single network and as a seed ensemble (two baselines); (vi) a multi-kernel weighted GP, a marginal-likelihood-fitted weighted sum of Matérn 3/2, Matérn 5/2, and RBF kernels at distinct length scales (with a bias term and white noise), instantiating the multi-kernel weighting proposed for obstructed fields [12]; and (vii) Poisson kriging, an ordinary-kriging GP whose per-point observation variance is set from Poisson counting statistics (proportional to the inverse local count; as the detector reports dose rate rather than raw counts, $1/(d+1)$ serves as the count proxy, clipped to $[1/3,3]$ times its median to avoid a degenerate near-zero nugget) [13]. The clip range and nugget scale shift Poisson kriging’s absolute error by up to ∼0.1 MAE on the sparsest runs but leave it far behind the physics-guided GP under every setting tried (Protocol B MAE $\ge 0.32$ vs. $\le 0.20$), so the comparison is insensitive to this heuristic. All methods are implemented with scikit-learn [22] and SciPy. All methods receive identical training data and are evaluated on identical held-out points. No method is evaluated on an internal split.

2.4. Evaluation Protocols

Simulation (truth-referenced). A collimated Cs-137 source (activity $7.4\times {10}^8$ Bq, detection efficiency 50%, leakage ${10}^{-4}$, half-angle 12°) defines the expected count rate at any location. Measured counts are drawn from a Poisson distribution, so that relative uncertainty varies naturally with position (≈0.02% at 1 m in-beam; ≈26% at the far out-of-beam regions), consistent with counting statistics [13,23]. The detection efficiency and leakage factor set the count level and the out-of-beam floor, respectively. The sampled counts are then log-transformed with the same ${log}_{10}(·+1)$ mapping used for the experimental dose rates, so simulation and experiment are scored on a common scale. Two hundred uniformly distributed samples are drawn per run, and six methods are compared: linear interpolation, GPR with a radial basis function (RBF) kernel, GPR with a Matérn 3/2 kernel plus bias, a multi-kernel weighted GP, Poisson kriging, and the MLP. Reconstructions are scored against the true field on a dense grid. The physics-guided GP is excluded from the simulation comparison because the simulated field is generated by the same collimated-beam model that its template embeds. Evaluating it here would test the method against its own generative assumptions and overstate its performance. It is therefore assessed only on the experimental data, where the template is at best an approximation to the real field. Robustness is assessed over 3 data seeds × 6 network seeds.

Experiments (two protocols). Protocol A (random hold-out, 80/20) measures interpolation quality near visited locations. Protocol B (spatial block cross-validation) partitions each survey area into a $4\times 4$ grid of blocks and performs group 5-fold cross-validation by block, measuring reconstruction in regions the trajectory never visited. Metrics are the mean absolute error (MAE), the RMSE, and the coefficient of determination ($R^2$) in the log domain, which capture average accuracy, large local errors (boundary failures), and spatial-structure recovery, respectively, plus a boundary F1 score: classification agreement, at held-out points, of in-beam status defined by a fixed threshold one decade above background (10 µSv/h). Method differences over the seven runs are tested by two-sided paired Wilcoxon signed-rank tests. To bound computation, GP kernel hyperparameters are estimated once per run (on the Protocol A training set) and reused across cross-validation folds. Template parameters are refitted within each fold from fold-training data only. Because the Protocol A training sample spatially overlaps the Protocol B test blocks, this reuse could in principle leak information favoring the GP-based methods. A control in which kernel hyperparameters were refitted within each fold from fold-training data only yielded indistinguishable results (physics-guided GP MAE $0.189\pm 0.048$ vs. $0.188\pm 0.047$; median $R^2$ 0.61 in both cases; per-run MAE differences $\le 0.005$, in both directions), so the shortcut has no material effect. Uncertainty calibration under Protocol B is assessed by the empirical coverage of $\pm 1\sigma$ and $\pm 2\sigma$ intervals and by the Spearman correlation between predicted $\sigma$ and absolute error. All results are reported as mean ± standard deviation (SD) over the seven runs. No run or seed is excluded. The experimental comparison additionally includes random forest (RF) and inverse-distance weighting (IDW) baselines, together with the physics-guided GP. RF and IDW are reported only for the experiments (being tangential to the smooth-prior-versus-MLP mechanism the simulation isolates), whereas the multi-kernel weighted GP and Poisson kriging are evaluated under both simulation and experiment.

3. Results

3.1. Simulation: Uniform Sampling Favors the MLP

Table 1 summarizes truth-referenced errors under Poisson noise (five data seeds). The MLP dominates: MAE $0.096\pm 0.014$ versus 0.31–0.46 for all other methods, $R^2$$0.947\pm 0.014$ versus 0.81–0.84. Figure 4 shows why: smooth-prior methods blur the beam edge, while the MLP reproduces the sector shape nearly exactly. This advantage is not an artifact of favorable initialization: across 18 runs spanning 3 data seeds and 6 network seeds, MLP $R^2$ remained within 0.932–0.966, never overlapping the GP range (0.808–0.850). Under uniform, well-covering samples, the absence of a smoothness prior is an asset, and MLP instability is negligible.

3.2. Protocol A: Method Parity at Visited Locations

At held-out points along the trajectory (Table 2), the ranking differs from the simulation: linear interpolation and random forest attain the best point metrics ($R^2$ 0.904 and 0.909), with inverse-distance weighting close behind (MAE 0.100, $R^2$ 0.875), and the MLP holds no advantage over the plain GP (MAE 0.112 vs. 0.108; all paired differences $p\ge 0.078$). The physics-guided GP matches the plain GP ($R^2$ 0.884), as expected, since near dense data the residual process dominates the template. Boundary F1 is likewise similar across all nine methods (0.859–0.905; best: linear interpolation $0.905\pm 0.038$; physics-guided GP $0.868\pm 0.069$). Seed sensitivity of the MLP correlates with data volume: the across-seed standard deviation of $R^2$ was ≤0.012 in the six denser runs but 0.069 (range 0.48–0.69) in the sparsest run (521 points). Random splits therefore measure neighborhood interpolation, for which most competent methods suffice.

3.3. Protocol B: Only the Physics-Guided GP Extrapolates

Spatial block cross-validation reverses the simulation ranking (Table 3). Every purely data-driven method degrades severely in unvisited regions; the MLP degrades worst (MAE $0.609\pm 0.305$; median $R^2=-5.7$), exhibiting the streak artifacts visible in Figure 5, and seed-ensembling mitigates but does not repair it (MAE 0.427). The physics-guided GP attains the highest median extrapolation $R^2$ ($0.61$, IQR $[0.51,0.69]$; mean $0.565\pm 0.198$, identical to variant M0 of the ablation in Table 4)—well above the strongest baseline, the multi-kernel weighted GP ($0.31$)—and achieves an MAE of $0.188\pm 0.047$: 35% below linear interpolation, 35% below the multi-kernel weighted GP, 41% below the plain GP, 32% below random forest, 49% below Poisson kriging, 50% below IDW, and 69% below the single MLP. The advantage is consistent over every baseline on all three metrics: the physics-guided GP had the lower error in all seven runs against each of the eight baselines (7/7), which is the most extreme outcome the two-sided paired Wilcoxon signed-rank test can return at $n=7$ (minimum attainable $p=0.016$). Boundary F1 in unvisited regions is $0.694\pm 0.134$ for the physics-guided GP versus 0.43–0.61 for all baselines (Table 3), with the higher per-run F1 against every baseline in at least six of the seven runs. Figure 6 juxtaposes the two protocols and makes the reversal visible at a glance: moving from random hold-out to spatial blocks inflates the mean MAE of every purely data-driven method by a factor of 2.8–5.4 (linear interpolation 3.4×, random forest 2.9×, multi-kernel weighted GP 2.8×, plain GP 3.0×, Poisson kriging 3.2×, IDW 3.7×, single MLP 5.4×), whereas the physics-guided GP degrades by only 1.7×, the smallest factor among all methods.

3.4. Source Parameters Recovery and Calibrated Uncertainty

The template parameters, fitted independently per run, are mutually consistent and physically plausible (Figure 7). For six of the seven runs the fitted beam direction falls within $-2^{\circ }$ to $+{10}^{\circ }$ of the true source orientation, and the fitted half-angle spans 12°–19°. Figure 8 illustrates a representative fit (run 3) overlaid on the survey map: the recovered source lies at the leading edge of the high-dose measurements and the fitted beam cone encloses the observed dose lobe, confirming that the recovered geometry is consistent with the raw data. That the fitted half-angles cluster above rather than around the certificate value of 12° is physically expected: the finite penumbra of a real collimator and residual scattering off surrounding structures broaden the apparent beam beyond its geometric aperture, so an intensity-based fit should recover an effective half-angle slightly larger than the nominal one. The nominal 12° therefore acts as a lower bound on the fitted values, consistent with what is observed.

Run 1 is the sole exception: its fitted direction is the mirror image of the true beam ($180.5^{\circ }$, plotted as $-179.5^{\circ }$ in Figure 7), while its half-angle ($17.5^{\circ }$) remains in the same range as the other runs. This is a symmetry-induced local minimum rather than a fit failure: when the fitted source position lies near the centroid of the trajectory, a beam and its ${180}^{\circ }$-reversed counterpart explain the trajectory-sampled intensities almost equally well, and the least-squares objective admits both as near-degenerate solutions. The magnitude parameters (source position and half-angle) are unaffected by this sign flip, the run’s reconstruction accuracy is on par with the others (Table 2 and Table 3), and the multi-start initialization—which seeds directions rotated by $\pm {90}^{\circ }$ and ${180}^{\circ }$—resolved the ambiguity in every cross-validation fold. The single-run whole-data fit shown here retains the mirrored solution only because it reports the global least-squares optimum without the fold-level restarts. The mapping thus yields a by-product of operational value: an estimate of source position, beam direction, and aperture without any dedicated inversion step, with the one known degeneracy identified and controlled.

Uncertainty calibration under Protocol B (Figure 9) shows empirical coverage of 0.70–0.85 for $\pm 1\sigma$ (nominal 0.683) and 0.84–0.95 for $\pm 2\sigma$ (nominal 0.954). The intervals are slightly overconfident in the worst run but remain usable for confidence mapping. Point-wise discrimination is weaker: the Spearman correlation between predicted $\sigma$ and absolute error ranged from $-0.58$ to $+0.06$, reflecting the stationary kernel’s distance-driven variance. The predictive-$\sigma$ map in Figure 5 illustrates the practical value nonetheless: uncertainty is low along the trajectory and high in unvisited regions, exactly the information a survey planner needs to direct supplementary measurements.

4. Discussion

4.1. Why the Rankings Differ

The simulation and experimental rankings are not contradictory; they are two regimes of one phenomenon. With uniform coverage, every location is near data, reconstruction is interpolation everywhere, and the MLP’s unconstrained flexibility wins at abrupt boundaries. With trajectory sampling, large areas are far from data and reconstruction is extrapolation: the MLP extrapolates without constraint (producing confident artifacts), a zero-mean GP shrinks toward a constant (blurring structure), linear interpolation clamps to nearest values (crude but bounded), and only a physically sensible mean function supplies the missing structure. Our preliminary tests indicated that an isotropic inverse-square prior offered only marginal improvement. This limited benefit is attributable to the anisotropic nature of the collimated radiation field, which is not adequately represented by an isotropic model. These findings underscore that a physics-based prior is effective only when it faithfully captures the underlying physical phenomenon.

4.2. Fairness of the Comparison

A potential objection is that the physics-guided GP wins the extrapolation comparison only because its mean function embeds the correct field family, whereas the purely data-driven baselines encode no structural prior. If any embedded template guaranteed the advantage, the comparison would be circular. To test this, we repeated the Protocol B evaluation with three deliberately altered mean functions, holding the residual-GP machinery and the spatial-block folds fixed (Table 4): an isotropic $1/r^2$ mean with no angular cutoff (M1); the same collimated template but with its beam direction forced ${90}^{\circ }$ away from the fitted value, i.e. correct in structure but wrong in a key parameter (M2); and a constant mean, which reduces the method to the plain GP (M3). The correct template (M0) reproduces the reported result ($R^2=0.565$, MAE 0.188). Every misspecified mean collapses to negative $R^2$ and roughly doubles the MAE (M1: $-1.43$; M2: $-1.98$), landing near the no-prior lower bound (M3: $-5.48$). The advantage is therefore not conferred by embedding a template per se, but by embedding the correct one; a wrong prior is no better than none. Crucially, prior correctness is not privileged external knowledge but an observable quantity, and the fit-quality gate that exploits it can itself be quantified. Across all 35 spatial-block folds (7 runs × 5 folds), the misspecified means incur a normalized fit cost of at least 1.45× (M1) and 2.20× (M2) that of the correct template on the same fold. Using the leave-one-run-out threshold defined in Section 2.2 ($\tau =1.3\times$ the largest correct-template cost among the remaining six runs), the gate detects 66 of 70 misspecified fits (94%) with 0 of 35 false alarms on correct fits. The four missed detections all occur in folds whose held-out blocks contain the beam core, so that the training data span only 65–80% of the run’s dynamic range and scarcely constrain any beam model. Detecting misspecification from such data is intrinsically hard, and this failure mode should be borne in mind operationally. With the gate in place, misspecification raises a runtime alarm rather than passing silently, prompting supplementary measurement or manual review; a conservative fallback to the plain GP (Section 2.3) is available as a default, though the stress test of Section 4 (Table 5) shows the freely-fitted template is not worse than that fallback even under misspecification, so the gate’s value lies in flagging rather than in switching. Consistent with this, the physics-guided GP holds no advantage under Protocol A (Table 2), where dense data make the prior largely irrelevant. Its benefit is confined to the data-starved extrapolation regime where a correct prior is the only available source of structure.

4.3. Gate Robustness to Unseen Misspecification Types

The gate calibration above (94% detection, zero false alarms) is measured on the single-source experimental data with deliberately altered mean functions (M1, M2). Because no multi-source, scatter-dominated, or occluded fields were surveyed, we assess whether the gate generalizes to such qualitatively different misspecifications through a controlled simulation study. Using the same collimated-source physics as Section 2.3, we generated 22 correctly-specified single-source fields and, as three unseen misspecification families, 20 realizations each of (i) multi-source superposition (two or three collimated sources), (ii) an isotropic scatter halo added to a single beam, and (iii) a single beam truncated by a wedge-shaped shadow (occlusion); Figure 10 shows a representative ground-truth field of each type. Each field was sampled along a trajectory-like path (six horizontal sweeps with unvisited bands between them) with Poisson counts, the single-source template was fitted, and its normalized fit cost $\overline{c}$ recorded. Because the simulated fields differ from the experimental data in dynamic range and sampling density, their absolute $\overline{c}$ values are not comparable to the experimental threshold. We therefore re-calibrate the gate internally on the simulated single-source fields, using the same $1.3\times$ leave-out rule, and report detection on the misspecified families (Table 5).

The gate flags 80% of multi-source fields—whose fit cost rises to roughly 80 times the correctly-specified level—at a 9% false-alarm rate on held-out single-source fields, confirming that the mechanism transfers to a structurally distinct misspecification it was never calibrated against. Occlusion is detected less often (40%): a freely-fitted single template still captures the visible lobe, so a localized shadow inflates the global cost only modestly, and this reduced sensitivity should be borne in mind operationally. The isotropic scatter halo is not flagged (0%), which is the correct behavior rather than a miss: the template’s leakage and background terms absorb the diffuse component, so the fit remains excellent (mean absolute error 0.052 versus 0.076 for the plain GP)—the method tolerates scatter rather than being broken by it. Across all four families the physics-guided GP reconstruction never exceeded the plain-GP fallback in mean error, so the gate is best read as an independent misspecification alarm—prompting operator caution or supplementary measurement—while the reconstruction itself degrades gracefully. A real multi-scenario validation remains necessary, but these results indicate the gate is not an artifact of the single dataset on which it was calibrated.

4.4. Relation to Prior Work

Read against the two traditions positioned in Section 1, our results sharpen that positioning rather than restate it. We deliberately retained the MLP as a strong, honestly-reported baseline, and the simulation confirms its genuine superiority under dense uniform sampling. The collimation-aware prior helps only where trajectory sampling leaves large unvisited regions, and among the geostatistical methods it outperforms both plain GPR and the multi-kernel weighted GP in that regime (Table 3). The resulting recommendation is conditional, not partisan: choose by sampling geometry. Where coverage is guaranteed (e.g., gridded post-accident drone sweeps), neural interpolators excel; where surveys are trajectory-constrained, as in the common robotic case, the physics-guided GP is preferable, and its variance output can drive informative path planning [8] to convert extrapolation into interpolation where it matters. We do not include a U-Net radiation-reconstruction baseline [15]: such models require a large corpus of paired sparse-input/dense-truth fields for supervised training, whereas each survey here yields a single field (seven in total), a regime in which a U-Net cannot be trained or compared fairly. We therefore treat it as complementary prior work rather than a baseline under our data budget.

4.5. Practical Implications

For ALARA-driven operations the asymmetry of boundary errors argues for methods that expose their ignorance. Two aspects of the uncertainty output must be distinguished. Its interval coverage is approximately calibrated ($\pm 2\sigma$ coverage 0.84–0.95), so aggregate confidence intervals can be used for regional risk flagging. Its point-wise discrimination is weak, however (Spearman $\sigma$–error correlation as low as $-0.58$ in the worst run), so predicted $\sigma$ should not be read as a reliable per-point error ranking. In practice the spatial $\sigma$ map remains useful because its dominant signal is coverage-driven—low along visited paths, high in unvisited regions—which is precisely what identifies where supplementary measurements would be most informative. The recovered source parameters offer a sanity check against facility records and a seed for formal source inversion [24]. They should be used as such rather than as a standalone source confirmation: as the run-1 mirror-direction degeneracy shows, the beam direction can flip when the source lies near the trajectory centroid, so parameter recovery must be accompanied by the multi-start procedure and a fit-quality gate, and the direction in particular should be cross-checked before operational use.

Radiation counts follow Poisson statistics, so a Poisson-likelihood GP might be expected to fit better than the log-Gaussian model and to sharpen the weak point-wise uncertainty ranking noted above. We tested this directly. Because the detector reports an integrated ambient dose-rate rather than raw counts—so a strict Poisson likelihood is not applicable to the experimental data—the comparison is run in simulation, where the counts are exactly Poisson. Over twelve collimated-source realizations under trajectory sampling, a Laplace-approximation Poisson GP with a log link (constant latent mean, Matérn 3/2 covariance) was compared against the log-Gaussian GP with the same kernel, both scored in the ${log}_{10}(\mathrm{rate}+1)$ domain against the true field (Table 6). The log-Gaussian model had the lower reconstruction error in all twelve realizations (MAE $0.066\pm 0.012$ vs. $0.123\pm 0.060$; paired Wilcoxon $p=0.0005$) and a higher $R^2$, while the Poisson GP improved the $\sigma$–error rank correlation only marginally and not significantly ($0.60\pm 0.11$ vs. $0.50\pm 0.09$; $p=0.13$), at comparable $\pm 2\sigma$ coverage. Two effects explain this. First, the log transform serves to compress the several-decade dynamic range, not to stabilize Poisson variance—for pure Poisson counts the variance-stabilizing transform is a square root, and the residual heteroscedasticity of ${log}_{10}(\mathrm{count}+1)$, largest at low counts, is absorbed by the white-noise kernel term. Second, the Laplace Poisson GP is less numerically stable under trajectory sampling and its mode reverts to the mean rate in unvisited regions, offering no extrapolation advantage over the log-Gaussian mean. A full Poisson (or heteroscedastic) treatment therefore remains a reasonable option for uncertainty ranking, but yields no accuracy or calibration gain here and is not adopted.

4.6. Limitations

First, all experiments involve a single dominant collimated source in one indoor environment. Multi-source superposition, scatter-dominated fields, and non-line-of-sight geometries remain untested on real data, and the template family must be chosen to match the scenario (with fallback to the plain GP when the fit is poor). A simulation stress test (Table 5) indicates that the fit-quality gate detects multi-source misspecification reliably (80%) but is less sensitive to localized occlusion (40%), while correctly tolerating isotropic scatter. A real multi-scenario validation nonetheless remains necessary. Second, $n=7$ runs limits statistical power: the paired Wilcoxon floor is $p=0.016$, and under a Holm–Bonferroni correction across the eight per-baseline comparisons none reaches $\alpha =0.05$ (the $n=7$ floor $p=0.016$ exceeds the corrected threshold $0.05/8=0.006$), so the evidential weight rests on the directional consistency (lower error in all seven runs against every baseline) rather than on nominal significance. Third, predicted $\sigma$ ranks point-wise errors weakly. Heteroscedastic or non-stationary kernels are natural next steps, whereas a Laplace Poisson-likelihood GP, tested in simulation (Table 6), improved point-wise ranking only marginally and non-significantly while degrading reconstruction accuracy, so it is an optional rather than a clearly beneficial extension [13]. Fourth, SLAM accuracy was evaluated in a motion-capture room rather than the radiation laboratory, and detector response combined with platform motion smears readings by ≈0.4 m, bounding the recoverable feature scale. Finally, the direction ambiguity observed in run 1 indicates that template fitting should in practice be accompanied by the multi-start procedure and the fit-quality gate quantified above (94% detection, zero false alarms on this dataset), whose known blind spot—training data that exclude the beam core—should be respected when interpreting gated results.

5. Conclusions

We investigated reconstruction of strong-gradient radiation fields from mobile robot surveys, with dual evaluation protocols separating interpolation from extrapolation, and reached three conclusions. (1) For the single dominant collimated source studied here, sampling geometry, rather than algorithm choice alone, determines method ranking: a coordinate-input MLP is robustly the best boundary reconstructor under uniform sampling (simulation $R^2$ 0.95, stable across 18 seed combinations), yet the worst extrapolator under trajectory sampling. (2) The proposed physics-guided GP, embedding a fitted collimated-beam template as the GP mean, achieved the highest median extrapolation $R^2$ ($0.61$; best baseline $0.31$), reduced unvisited-region MAE by 32–69% relative to all eight baselines (lower error in all 7 runs against every baseline), delineated the beam boundary most accurately (F1 0.69), and recovered source position, direction, and aperture consistent with the actual apparatus across seven runs. (3) Its predictive uncertainty was approximately calibrated ($\pm 2\sigma$ coverage 0.84–0.95) and maps directly onto operational needs, flagging untrustworthy regions and prioritizing supplementary measurements, thereby supporting the precise, decision-ready boundary delineation that field operations and emergency response require. Future work will extend the template family to multiple and unknown source types, exploit the variance output for active survey planning, deepen uncertainty calibration, and generalize to three-dimensional fields.