Sensors
- TanDEM-X: Delivers a global DEM at 12 m posting (0.4 m relative vertical accuracy over flat terrain). At site scale, curvature and slope derivatives from this DEM resolve subtle earthwork remnants and ploughzone hollows that correlate with buried feature distributions. No revisit in the conventional sense: the global DEM is a single-epoch product, though new acquisitions for change detection are possible by arrangement.
- Sentinel-2 MSI: 13 spectral bands, 10 m resolution in VNIR, 20 m in SWIR, 5-day revisit at mid-latitudes. NDVI, NDWI and red-edge indices derived from Sentinel-2 time series capture soil-moisture and crop-stress patterns above buried features. Free, open archive from 2015.
- WorldView-3: 0.31 m panchromatic, 1.24 m multispectral VNIR, 3.7 m SWIR (8 SWIR bands). At this resolution, individual feature boundaries visible in surface soil or crop colour can be digitised directly and used as hard constraints in the interpolation model. Tasked commercially; archive depth from 2014.
- ALOS World 3D (AW3D30): Freely available global DEM at 30 m posting, useful as a coarser complement or cross-check to TanDEM-X where the latter requires a commercial licence. Vertical accuracy approximately 5 m RMSE globally, degrading in vegetated terrain.
Why a handful of transects is never enough, and why that is not the surveyor's fault
Ground-penetrating radar is among the most informative non-invasive tools in field archaeology. A 400 MHz antenna resolves features at depths of 0.5 to 2 m with sub-decimetre vertical precision. The problem is areal coverage. A single operator pulling a cart antenna covers perhaps 0.5 to 1 hectare per day at survey-grade transect spacing of 0.5 m. A 10-hectare Roman villa complex therefore demands roughly two weeks of continuous fieldwork before the first interpolated time-slice map can be drawn. Most projects cannot afford that, so surveys are designed as systematic transect grids covering 20 to 40 percent of the site, with the remainder inferred.
That inference is where satellite data become genuinely useful. The argument is not that orbital sensors can see through soil. They cannot, at optical or near-infrared wavelengths, and the SAR subsurface penetration case is handled elsewhere in this library. The argument is that buried features leave surface expressions, however faint, and that micro-topographic curvature, soil-moisture anomalies and crop-stress patterns are spatially continuous covariates that correlate with subsurface structure. If you can model that correlation where GPR has sampled it, you can extrapolate it where it has not.
The covariates: what each satellite layer actually contributes
TanDEM-X at 12 m posting is the workhorse for micro-topography. Second-derivative curvature computed from the DEM highlights subtle positive and negative relief, the kind of 0.2 to 0.5 m mounds and hollows that correspond to wall footings, ditches and floor surfaces that have differentially compacted the overburden. Published work at Roman villa sites in the UK and Germany has shown statistically significant correlation between TanDEM-X curvature anomalies and GPR-detected wall lines at depths of 0.4 to 1.2 m. The 12 m posting is a genuine constraint: features narrower than roughly 24 m are unlikely to be resolved, and the DEM conflates wall and ditch signals at sub-grid spacing.
Sentinel-2 contributes two distinct layers. First, a seasonal NDVI or red-edge chlorophyll index stack captures crop-mark cycles across multiple growing seasons. A single dry summer may not produce legible marks; a five-year median composite from the Sentinel-2 archive (available from 2015) substantially improves signal-to-noise. Second, SWIR-derived soil-moisture indices applied to bare-soil windows reveal differential moisture retention above features, a proxy for buried voids and compacted fills. WorldView-3 adds spatial precision: at 1.24 m multispectral resolution it can resolve individual wall-width anomalies in crop colour, providing hard-constraint training points for the geostatistical model rather than just smooth covariates.
The interpolation workflow: from sampled transects to site-wide probability
The standard published approach treats GPR amplitude or feature-presence at each transect cell as the dependent variable and the satellite-derived raster stack as predictors. Regression kriging is the most commonly applied method: a deterministic regression on the covariates captures the trend, and ordinary kriging of the residuals captures spatial autocorrelation not explained by the satellite layers. The output is a continuous probability surface, typically expressed as the likelihood that a buried feature of a given type exists at each grid cell.
Random forest and gradient-boosted tree classifiers have been applied as non-parametric alternatives, particularly where the relationship between surface expression and subsurface feature is non-linear or varies across a heterogeneous site. These methods require more GPR training samples to avoid overfitting. A minimum of roughly 30 to 50 GPR-confirmed feature instances is a reasonable lower bound for a classifier to generalise reliably, though published studies vary in their reported thresholds.
The probability surface is not a substitute for a GPR map. It is a prioritisation tool: it tells a project manager where to concentrate the next season's ground survey, or where to focus targeted excavation trenches. Used that way, it typically reduces the area requiring full GPR coverage by 30 to 50 percent in published case studies, though those figures depend heavily on site type and covariate quality.
Honest limits: where extrapolation fails and why
Extrapolation uncertainty is not linear. It grows roughly with the square of transect spacing, and beyond two to three times the GPR transect interval the probability estimates carry confidence intervals wide enough to be operationally useless. At a site surveyed on 4 m transect spacing, extrapolation is reasonably constrained within 8 to 12 m of a sampled transect. At 20 m or more from any sampled line, the satellite covariates are doing most of the work and the GPR is providing little more than a distant anchor.
Several site conditions defeat the method entirely. Deep alluvial burial, where features lie beneath 1.5 m of structureless sediment, removes the surface expression that satellite sensors depend on. Intensive modern land use, particularly deep ploughing or construction, destroys the micro-topographic signal in the DEM and homogenises soil colour. Woodland canopy prevents both crop-mark formation and reliable DEM derivation from optical stereo; TanDEM-X penetrates canopy partially but the resulting DEM represents a surface somewhere between ground and crown. In all these cases the covariate layers lose predictive power and the extrapolation should not be trusted.
There is also a subtler problem: the satellite covariates are correlated with surface processes, not subsurface ones. A soil-moisture anomaly may reflect a modern field drain rather than a Roman hypocaust. The model cannot distinguish these without GPR ground-truth, which is precisely what is sparse. Any delivered probability surface should carry an explicit uncertainty layer, not just a point estimate.
Published precedents and what they actually demonstrated
The fusion approach has been applied and published at Bronze Age field systems in the UK, where Sentinel-2 seasonal composites and a lidar-derived DEM (a closer relative of TanDEM-X than an identical substitute) were used to extrapolate magnetometry and GPR transects across a 200-hectare landscape. The resulting probability maps correctly predicted the location of subsequently excavated field boundaries in roughly 70 percent of tested cells, with false-positive rates around 25 percent. Those figures are encouraging for a prospection tool but would be unacceptable as a definitive site map.
Roman villa work in Germany and Italy has used TanDEM-X curvature specifically, finding that wall-line detection rates in the covariate model matched GPR-confirmed walls at rates of 60 to 75 percent when training data covered at least 25 percent of the site area. Detection rates fell to below 50 percent when training coverage dropped below 15 percent. These numbers come from peer-reviewed work published in journals indexed in the Remote Sensing (MDPI) and Archaeological Prospection literature; they are not Satellize figures.
Satellize has applied related covariate-stack methods in the Kingdom of Tonga crop-estimation programme, where geostatistical interpolation of sparse ground-truth points using Sentinel-2 indices is the operational core of the workflow. The archaeological extrapolation problem is structurally similar, though the agronomic version benefits from denser and more regularly distributed ground samples than most heritage budgets allow.
Commissioning a fusion analysis: what to prepare before the satellite data arrive
The satellite data are the easier part to obtain. TanDEM-X DEMs are available commercially through DLR; Sentinel-2 imagery is free from the Copernicus Data Space. The harder preparation is the GPR archive itself. Transect logs need georeferencing to better than 0.5 m horizontal accuracy, which requires differential GPS or total-station control rather than handheld GNSS. Feature picks need to be attributed consistently: a wall reflection and a pipe reflection look similar in a radargram, and a mixed training set degrades the model sharply.
Site-specific covariate selection matters more than generic workflows suggest. A site on chalk downland will respond well to SWIR soil-moisture indices. A site on heavy clay in a humid climate may show no useful bare-soil window across an entire Sentinel-2 archive. Scoping that covariate availability before committing to the fusion approach saves significant analytical effort. A preliminary spectral time-series review of the Sentinel-2 archive over the site, combined with a TanDEM-X curvature preview, is a sensible first commission before any GPR data are incorporated.
Typical figures
| TanDEM-X DEM posting | 12 m (0.4 m relative vertical accuracy, flat terrain) |
| Sentinel-2 spatial resolution | 10 m VNIR, 20 m SWIR, 60 m atmospheric bands |
| Sentinel-2 revisit | 5 days at mid-latitudes (combined Sentinel-2A and 2B) |
| WorldView-3 resolution | 0.31 m pan, 1.24 m VNIR, 3.7 m SWIR |
| Sentinel-2 archive depth | From 2015 (2A), 2017 (2B); freely accessible via Copernicus Data Space |
| TanDEM-X minimum resolvable relief | Approximately 0.2 m relative height difference over flat terrain at 12 m posting |
| Minimum GPR training coverage for reliable extrapolation | Published studies suggest 25 percent site coverage; below 15 percent, model accuracy degrades substantially |
| Extrapolation reliability range | Operationally useful within 2–3x the GPR transect interval; degrades rapidly beyond that |
| Covariate stack delivery format | GeoTIFF raster stack, EPSG-specified projection, 10 or 12 m native grid |
| Probability surface output format | GeoTIFF with accompanying uncertainty layer; compatible with QGIS, ArcGIS, GRASS |
Analytics Satellize can run
| TanDEM-X curvature anomaly map | Second-derivative terrain analysis (profile and plan curvature) from 12 m DEM; thresholded to isolate micro-relief consistent with buried structural features | GeoTIFF raster layer with anomaly polygons as GeoPackage overlay, annotated with curvature magnitude and spatial extent |
| Sentinel-2 multi-season spectral composite stack | Median compositing across 5-year archive; NDVI, red-edge chlorophyll index and SWIR soil-moisture index computed per pixel; bare-soil and crop-mark windows identified separately | Multi-band GeoTIFF stack with per-band seasonal statistics; PDF summary of spectral anomaly windows identified |
| GPR-satellite covariate correlation report | Pearson and Spearman correlation between GPR feature-presence labels and each satellite covariate layer at sampled transect locations; identifies which covariates carry predictive signal at the specific site | PDF report with correlation matrices, scatter plots and ranked covariate importance; informs whether full fusion model is warranted |
| Site-wide subsurface feature probability surface | Regression kriging or random forest trained on GPR transect labels with satellite covariate stack as predictors; cross-validated using leave-one-transect-out procedure | GeoTIFF probability raster (0–1 scale) with accompanying kriging variance or prediction interval layer; GeoPackage of high-probability zones for field prioritisation |
| Survey prioritisation grid | Expected information gain calculated per unsurveyed grid cell using probability surface and uncertainty layer; cells ranked by reduction in uncertainty per unit survey effort | GeoPackage grid with ranked survey-priority scores; exportable as field navigation waypoints |
| WorldView-3 high-resolution feature constraint layer | Manual and semi-automated digitisation of spectral anomalies in 1.24 m VNIR imagery; SWIR band ratios applied to identify differential soil composition at feature boundaries | GeoPackage polygon layer of visually confirmed surface anomalies used as hard constraints in geostatistical model; attribute table includes confidence rating per feature |
Who does the work
We can get this done for you. Satellize runs its own analyst desk and a strong science team. You do not buy a data feed and work out what it means; our people source the imagery, run the analysis described on this page, and hand you the answer with its confidence limits stated. Discuss this requirement.