GTA · Placer Mining Detection Overview Data Z-Score Fusion Co-Reg Temporal PCA Results

Detecting Placer Gold Mining from Sentinel-2 & Sentinel-1

A multi-method, multi-sensor change-detection pipeline for the Uylgaiin Dugang site, Mongolia.

Abstract

Placer gold mining in northern Mongolia leaves distinctive optical and radar signatures along river corridors: vegetation stripping, dredge ponds, exposed tailings, sediment plumes, and persistent loss of InSAR coherence. We combine freely-available Sentinel-2 multispectral and Sentinel-1 SAR imagery into a layered detection pipeline at the Uylgaiin Dugang site (49.78°N, 107.84°E). The pipeline starts from a multi-index Z-score change detector, layers in SAR coherence-loss fusion, applies AROSICS sub-pixel co-registration to remove geometric false positives along riverbanks, and culminates in a temporally dense time-series cube analysed by joint spectral-temporal PCA with a per-pixel change-point detector that identifies the moment a pixel switches from one terrain class to another. Cloud handling is unified across stages: SCL masking, dilated cloud edges, frame drop based on AOI valid fraction, and a temporal MAD outlier filter.

This document is interactive. Each method section has a "Theory of the method →" button that opens a deep-dive popup with the underlying math, algorithm pseudo-code, and a small SVG diagram. Click any of them whenever you want the derivation; ignore them and read the prose if you do not.

§ 1
Site & Problem
Uylgaiin Dugang AOI, what placer mining looks like from space, why off-the-shelf change detectors fail.
§ 2
Data Sources
Sentinel-2 L2A, Sentinel-1 SLC/GRD via HyP3, ALOS PALSAR L-band, optional PlanetScope SuperDove.
§ 3
Multi-Index Z-Score
NDVI, MNDWI, BSI, NDTI, NDWI, NBR — standardized differences and fused anomaly maps.
§ 4
Optical × SAR Fusion
Adding InSAR coherence loss as a confirmatory signal in cloud-prone summer months.
§ 5
Co-Registration
AROSICS Global as the production sub-pixel aligner; benchmarked against six alternatives.
§ 6
Temporal PCA & Change Points
Joint spectral-temporal eigendecomposition; per-pixel sliding-window dominant-PC tracking.
§ 7
Results Gallery
Time-lapse, drilldowns, riverbed mask, SAR & coregistered comparisons.

§ 1 The Site & the Problem

Placer gold mining strips a river corridor of vegetation, dredges the alluvium, dumps the rejected gravel as tailings, and leaves behind a chain of turbid settling ponds. From space, every one of those steps is visible — if you know which spectral indices to combine and how to filter out the natural seasonal cycle. Our test site is Uylgaiin Dugang in Töv Province, Mongolia (49.7767°N, 107.8357°E), an active placer site whose extent grew measurably between 2017 and 2020. The AOI is a tight 5×4.5 km box in EPSG:32648 (UTM zone 48N).

ParameterValue
Site nameUylgaiin Dugang Mining
Centre (lat, lon)49.7767°N, 107.8357°E
AOI bbox (EPSG:4326)[107.80, 49.76, 107.87, 49.80]
Reference gridEPSG:32648 (UTM 48N), 10 m
Mining seasonApr–Dec, peak Jun–Aug
Baseline windowJun–Aug 2017
Peak-activity windowJun–Aug 2020

Why off-the-shelf change detection is not enough

A naive Δ-NDVI between two summer composites flags everything: snow patches that haven't melted, cloud shadows the SCL classifier missed, fields harvested at different dates, and the natural sub-pixel drift of the Copernicus reference grid between processing baselines. Below the noise floor sits the actual mining signal, which has three properties that distinguish it: (a) it sits along a river corridor, (b) it persists rather than recovering, and (c) the optical signature is corroborated by SAR coherence loss because the dredge head physically moves. The pipeline below is built around these three properties.

§ 2 Data Sources

All sensors used are free at the point of access. The pipeline scaffolds a fifth (PlanetScope 3 m via the Planet Education & Research program) for the user's pending application, but every figure in this document is generated from the four free sources alone.

SensorBands / ModeResolutionRevisitRole in pipeline
Sentinel-2 L2A 13 bands incl. SWIR (B11, B12) 10–20 m ~5 d Optical backbone — spectral indices, time-series cube
Sentinel-1 SLC C-band IW, VV+VH, complex ~5×20 m ~12 d InSAR coherence loss via ASF HyP3
Sentinel-1 GRD C-band IW backscatter 20 m ~12 d dB-change for "new bright reflectors" (equipment, structures)
ALOS PALSAR mosaic L-band HH/HV annual 25 m annual Long-term mining-area expansion trend
PlanetScope SuperDove (scaffolded) 8 bands, no SWIR 3 m ~daily Pending Planet E&R access — sub-pit detail
Sentinel-2 L2A 10 m · 6 bands Sentinel-1 SLC via HyP3 InSAR Sentinel-1 GRD backscatter dB ALOS PALSAR L-band annual Indices + Z-score NDVI / MNDWI / BSI ... Coherence loss before − after ΔdB equipment mask new strong scatterers Long-term L-band trend multi-year expansion Optical × SAR fusion weighted Z-score boost Co-registration AROSICS Global Detection polygons GeoJSON · GeoPackage Temporal-PCA cube change-time map · drilldowns
Pipeline overview — sources (left) flow through per-sensor processing into the fusion / co-registration stages (centre) and out to the two main products: detection polygons and the temporal-PCA cube.

§ 3 Multi-Index Z-Score Change Detection

The first detector layer compares two Sentinel-2 median composites (Before, After), computes six spectral indices for each, takes the differences, standardizes per index, and weighted-sums them into a single anomaly map. The map is thresholded, morphologically cleaned, vectorised into polygons, and given a confidence boost if the polygon sits near a pre-existing water feature (placer mining hugs rivers).

IndexFormulaMining signatureSign in fusion
NDVI(NIR − R) / (NIR + R)Vegetation stripping → large negative Δinverted
MNDWI(G − SWIR1) / (G + SWIR1)New ponds, dredge water → large positive Δpositive
BSI((SWIR1+R) − (NIR+B)) / (sum)Bare soil / tailings → positive Δpositive
NDTI(R − G) / (R + G)Sediment plumes / turbidity → positive Δpositive
NDWI(G − NIR) / (G + NIR)Open water increase → positive Δpositive
NBR(NIR − SWIR2) / (NIR + SWIR2)General disturbance / burn → negative Δinverted

§ 4 Optical × SAR Fusion

The optical Z-score detector struggles in two situations: (a) summer cloud cover blanks out critical Jun–Aug acquisitions, and (b) early-season false positives along snowmelt edges. Sentinel-1 InSAR coherence is independent of clouds and cleanly highlights "actively-disturbed ground": stable steppe maintains coherence >0.7 between repeat passes; an active dredge drops below 0.3. We fuse the two by adding a SAR boost on top of the optical Z-score.

SAR + optical 4×3 panel analysis
SAR + optical analysis figure for the Uylgaiin Dugang AOI. Row 1 — Before / After RGB and a coherence-loss overlay (cyan polygons = optical detections). Row 2 — Before / After SAR coherence and the coherence-loss difference map (red = disturbance). Row 3 — Before / After SAR amplitude (dB) and the amplitude change map (red = new strong reflectors, a proxy for equipment or structures). Row 4 — Bright-stretched RGB to bring out tailings and brightness change. Cyan/lime outlines mark the optical Z-score detections; the SAR coherence loss independently confirms several of them.

§ 5 Sub-Pixel Co-Registration

Sentinel-2 L2A products inherit Copernicus's global reference grid, but residual sub-pixel drift of even ~3 m produces visible "edge halos" along narrow rivers in the ΔMNDWI map and is the dominant source of false positives in riverine change detection. We benchmark seven co-registration methods in a separate CoRegistration project and pick AROSICS Global as the production aligner: it estimates a single rigid (Δx, Δy) translation in image-space using FFT phase correlation on a stable matching band (NIR), so it removes geometric mis-registration without absorbing real channel motion into a deformation field.

Uncorrected vs co-registered change detection
Z-score detections before vs after AROSICS Global. Same threshold, same Before/After pair. After alignment, edge-halo false positives along the riverbank shrink (total area drops 13%), and mean confidence rises from 0.86 to 0.93 because the retained pixels represent real change rather than misalignment of a third of a pixel.

§ 6 Temporal PCA & Change-Point Detection

The before/after pipeline above only sees two snapshots per year. To exploit the full Sentinel-2 cadence we build a temporally dense data cube: every available S2 L2A scene over a one-year window, SCL-masked, AROSICS-aligned to a single reference, cloud-cleaned (dilation + frame drop + temporal MAD outlier rejection), and resampled onto a regular 5-day temporal grid. The result is a 4-D array $X \in \mathbb{R}^{T \times B \times H \times W}$ with $T \approx 36$ timesteps, $B = 6$ spectral bands, and ~243k pixels.

We then run joint spectral-temporal PCA: per pixel, stack all bands and times into a length-$TB$ "spectral-temporal fingerprint", standardize, build the $TB \times TB$ covariance, and eigendecompose. The principal components are length-$TB$ patterns that capture the dominant modes of how a pixel evolves over the year. Pixels sharing the same dominant PC share a temporal land-cover class — water, vegetation, snow, bare soil, mining. A per-pixel sliding-window detector then finds the time at which the dominant PC switches, which is exactly the question "when did this pixel stop behaving like X and start behaving like Y?"

Sentinel-2 RGB time-lapse over the AOI
1 / 36
Sentinel-2 RGB time-lapse over the 6-month test cube (test_2020 label, Apr–Sep 2020 after cloud cleanup). Snow recedes, vegetation greens up and senesces; you can spot dredge ponds appearing in the central river corridor. Use the controls to play/pause, step frame-by-frame, or scrub through the sequence.
Top-margin change pixel drilldown
Drilldown on the global top-margin change pixel. Top — raw 6-band intensity through the year, with the red change point and the blue/red dashed lines marking the four cleanest before/after frames the algorithm picked. Middle — windowed dominant joint-PC trajectory, switching from PC1 to PC7 at the change time. Bottom — the four 200×200 pixel zoom panels with a yellow crosshair pinpointing the target pixel without occluding its colour.
Riverbed mask + drilldown locations
Riverbed mask derived from per-pixel temporal MNDWI counts (≥7 of 36 timesteps water-like, after morphological cleanup), overlaid in blue on a mid-year true-colour basemap. Yellow ring = the global top-margin change pixel; cyan ring = the strongest change pixel that lies inside the riverbed mask.
Riverbed drilldown
Drilldown on the strongest riverbed change pixel. Same plot anatomy as above. The 200×200 zoom panels make it clear this pixel sits on the river meander; the temporal-band plot shows the spring snow-melt collapse and the post-melt rebound where the dredge head moved through.

§ 7 Results & Validation Numbers

Headline numbers (test_2020 cube, 6-month window)

Per-detection validation thumbnails
Per-detection validation panel. For every Z-score detection polygon we extract a Before / After / Δ-NDVI / Δ-MNDWI thumbnail. This is the manual-QA view used to label the example polygons in detections.geojson.

Where the cloud filter mattered most

Before adding the cloud-cleanup stage to timeseries_cube.py, the joint PCA was dominated by cirrus and cloud-edge artefacts: PC1 captured only 30% of the variance and the change-point detector flagged 35% of the AOI. The drilldown popups looked like clouds. After dilating the SCL mask by 2 px, dropping any frame with <85% valid AOI pixels, and rejecting per-pixel temporal outliers above 4·MAD, PC1 jumped to 73% of variance, the change-pixel count fell to 1.26% of the AOI, and the surviving drilldowns are real ground transitions instead of cloud edges. The per-frame KEEP/DROP decisions are logged to cube_meta.json for audit.

§ 8 Next Steps

  1. Run on the full year. The pipeline scales linearly: 70–75 unique-date Sentinel-2 scenes over a 1-year window are downloaded, AROSICS-aligned, cloud-cleaned, and fed into the temporal-PCA cube without any code changes — just a bigger label.
  2. PlanetScope SuperDove (3 m, daily) via the Planet Education & Research program. The downloader and PlanetScope-specific change detector are already scaffolded under download_planetscope.py and change_detection_planetscope.py; they wait on API key issuance. With 3 m daily we move from "blob detection" to per-dredge tracking.
  3. Refine the AROSICS local mode in the companion CoRegistration project. The temporal-PCA cube is currently aligned via AROSICS Global (a single rigid translation); swapping in the locally-varying Bayesian-tiled aligner would handle pairs across very different baseline reprocessings.
  4. Label the joint PCs by hand. The change-point detector reports "PC1 → PC7" transitions; an interpretive table mapping each PC to a temporal class ("snow-melt", "stable steppe", "disturbed bare ground" ...) would allow filtering for the specific transitions diagnostic of mining onset.

Theory — Multi-Index Z-Score Change Detection

1. Per-index difference

For each spectral index $I$ (NDVI, MNDWI, BSI, NDTI, NDWI, NBR), compute the per-pixel difference between the After median composite and the Before median composite:

$$\Delta I(x, y) \;=\; I_{\text{after}}(x, y) \;-\; I_{\text{before}}(x, y)$$

2. Z-score standardization

Each Δ-index map is standardized by its scene-wide mean and standard deviation, so that "anomalous" pixels can be compared across indices with different natural scales:

$$Z_I(x, y) \;=\; \frac{\Delta I(x, y) \;-\; \overline{\Delta I}}{\sigma_{\Delta I}}$$

where $\overline{\Delta I}$ and $\sigma_{\Delta I}$ are computed over all valid pixels in the AOI.

3. Weighted fusion with sign correction

Each index has a known sign for the mining signature: vegetation indices (NDVI, NBR) drop, water and bare-soil indices (MNDWI, BSI, NDTI, NDWI) rise. We invert the sign of the vegetation indices and compute a single fused anomaly map as a weighted sum:

$$F(x, y) \;=\; \sum_{I \in \mathcal{I}} w_I \cdot s_I \cdot Z_I(x, y), \qquad s_I = \begin{cases} -1 & I \in \{\text{NDVI, NBR}\} \\ +1 & \text{otherwise} \end{cases}$$

The default weights $w_I$ favour the most discriminating indices: NDVI = 0.25, MNDWI = 0.20, BSI = 0.20, NBR = 0.15, NDTI = 0.10, NDWI = 0.10. They are tunable in change_detection.py.

4. Threshold + connected-component labelling

Apply a fused-Z threshold (default 2.5), morphologically close then open the binary mask with an 8-connectivity structuring element to remove speckle, then label connected components. Polygons below 0.5 hectare are discarded.

5. Riparian proximity boost

Placer mining sits on rivers. We compute the Euclidean distance from every detection centroid to the nearest pixel where the Before-period MNDWI exceeds zero. Detections within 2 km get a +0.15 confidence boost; those beyond 5 km lose 0.15. The mining signature is so strongly riparian that this single geometric prior is enough to reliably reorder the detection list.

References

Theory — Optical × SAR Fusion

1. SAR coherence loss

For two co-registered Sentinel-1 SLC acquisitions, the inteferometric coherence is

$$\gamma(x, y) \;=\; \frac{\left| \langle s_1 s_2^{*} \rangle \right|} {\sqrt{\langle |s_1|^2 \rangle \,\langle |s_2|^2 \rangle}}$$

where the angle brackets denote spatial averaging over a small window. Stable terrain (rocks, buildings, dry steppe) has $\gamma$ close to 1; vegetation has $\gamma \sim 0.4$; actively disturbed ground (a working dredge) drops to $\gamma < 0.3$. Coherence loss between an early-summer pair and a late-summer pair is the SAR fingerprint of mining activity:

$$\Delta\gamma(x, y) \;=\; \gamma_{\text{before}}(x, y) - \gamma_{\text{after}}(x, y)$$

2. Resampling SAR to the optical grid

HyP3 InSAR products are typically multi-looked to ~80 m. We bilinearly resample the coherence-loss map onto the optical 10 m UTM grid using rasterio.warp.reproject so that the fusion is pointwise in pixel space.

3. Additive fusion

Rather than multiply (which would suppress optical-only or SAR-only detections), we add a SAR boost on top of the optical Z-score, capped at +1.5 standard deviations:

$$F_{\text{fused}}(x, y) \;=\; F_{\text{opt}}(x, y) \;+\; w_{\text{sar}} \cdot \min\!\left(\frac{\Delta\gamma(x, y)}{\tau_\gamma},\; 1\right) \cdot 1.5$$

with $w_{\text{sar}} = 0.35$ and $\tau_\gamma = 0.3$. The fused score is then thresholded at 2.0 (slightly looser than the optical-only 2.5) and vectorised the same way.

Why additive, not multiplicative?

Multiplicative fusion would require both signals to fire simultaneously to detect anything — which destroys the very thing SAR is good for: penetrating cloud cover that hides the optical signal. Additive fusion lets either modality detect a polygon independently while giving extra weight to polygons that both confirm.

References

Theory — AROSICS Global Co-Registration

1. Phase cross-correlation in the Fourier domain

Given a reference image $f(x, y)$ and a target image $g(x, y) = f(x - \Delta x, y - \Delta y)$ that differs only by a rigid translation, the Fourier shift theorem says $G(u, v) = F(u, v) \cdot e^{-2\pi i (u\Delta x + v\Delta y)}$. The normalized cross-power spectrum

$$R(u, v) \;=\; \frac{F(u, v) \, G^{*}(u, v)}{\left| F(u, v) \, G^{*}(u, v) \right|}$$

is therefore exactly $e^{2\pi i (u\Delta x + v\Delta y)}$. Its inverse Fourier transform is a Dirac at $(\Delta x, \Delta y)$, which is recovered as the brightest peak in the correlation surface. Sub-pixel accuracy is obtained by upsampling the peak via a matrix-multiply DFT (Guizar-Sicairos et al., 2008).

2. Why a single global shift

A more flexible aligner (dense optical flow, locally-varying affine) can pull real changes — a moving river meander, a new dredge pond — into its deformation field, erasing the very signal we want to detect. A single rigid (Δx, Δy) translation cannot do that. It corrects sub-pixel drift of the Copernicus reference grid without ever absorbing land-cover change.

3. Choosing the matching band

AROSICS picks one band (here NIR / B08) for the FFT match. NIR is the right choice for placer mining sites because:

4. Application to all bands

Once the (Δx, Δy) shift is estimated from the matching band, AROSICS applies the same rigid shift to every band of the multispectral stack and writes the shifted output as a fresh GeoTIFF. We then resample back onto the project reference grid with cubic interpolation to keep all downstream rasters pixel-aligned.

Reference NIR B08, before Target NIR B08, after Phase correlation FFT, sub-pixel peak → (Δx, Δy) Apply shift to all 6 bands
AROSICS Global flow — one (Δx, Δy) is estimated on the NIR band and applied identically to every other band.

References

Theory — Temporal PCA & Change-Point Detection

1. The data cube

Stack every cloud-cleaned, AROSICS-aligned, gap-filled Sentinel-2 scene into a 4-D tensor

$$X \;\in\; \mathbb{R}^{T \times B \times H \times W}$$

where $T$ is the number of regular 5-day timesteps over the year, $B = 6$ is the number of spectral bands (B02, B03, B04, B08, B11, B12), and $H \times W$ is the AOI grid. For each pixel $(y, x)$ we have a length-$TB$ "spectral-temporal fingerprint":

$$\mathbf{x}_{(y,x)} \;=\; \big(X_{0,0,y,x},\; X_{0,1,y,x},\; \ldots,\; X_{T-1,B-1,y,x}\big)^{\!\top}$$

2. Standardisation

Per pixel and per band, subtract the temporal mean and divide by the temporal standard deviation. This makes the PCs reflect patterns of change, not absolute brightness, so a bright tailings pile and a dim river bed can land in the same PC if they evolve similarly through the year.

3. Joint spectral-temporal PCA (T-mode EOF)

Reshape the standardised cube into a $(TB) \times N$ matrix $M$ where $N = HW$ and each column is one pixel's full fingerprint. The $TB \times TB$ joint covariance is

$$C \;=\; \frac{1}{N}\, M\, M^{\top}$$

Eigendecompose $C = U \Lambda U^{\!\top}$. The columns of $U$ are orthonormal length-$TB$ spectral-temporal eigenvectors — the dominant modes of how a pixel evolves over the year. Project every pixel onto the top $K$ PCs:

$$A \;=\; U_{:,1:K}^{\top}\, M \;\;\in\;\; \mathbb{R}^{K \times N}$$

Each column of $A$ is a pixel's loading on the top-$K$ PCs. The pixel's dominant PC is $\operatorname{argmax}_{k}\, |A_{k,i}|$, and pixels sharing a dominant PC have a similar annual signature — a temporal land-cover class.

Cube X T × B × H × W ~36 × 6 × 465 × 522 Reshape M (TB) × N ~216 × 243k C = M Mᵀ / N (TB) × (TB) eigh → U, Λ Top-K PCs U[:, :K] + scores A
Joint spectral-temporal PCA: every pixel becomes a length-$TB$ vector; covariance over $N$ pixels gives $TB$ eigenvectors, each a "fingerprint" of one terrain class.

4. Sliding-window dominant-PC tracking

A pixel that switches land-cover class part-way through the year has a bimodal fingerprint — vegetation-like in $[0, t^{*}]$, disturbed-like in $[t^{*}, T]$. To pinpoint $t^{*}$ we slide a window of $W$ timesteps across each pixel's fingerprint. In the joint case the window covers $W \cdot B$ consecutive flat elements starting at $t \cdot B$. Inside each window we project the pixel slice onto the same-window slice of every PC, normalised to unit length:

$$p_t^{(k)} \;=\; \left| \hat{\mathbf{x}}_t \cdot \hat{\mathbf{u}}_t^{(k)} \right|, \qquad \hat{\mathbf{x}}_t = \frac{\mathbf{x}_t}{\|\mathbf{x}_t\|},\;\; \hat{\mathbf{u}}_t^{(k)} = \frac{\mathbf{u}_t^{(k)}}{\|\mathbf{u}_t^{(k)}\|}$$

The dominant PC at window centre $t$ is $k_t^{*} = \operatorname{argmax}_{k}\, p_t^{(k)}$, and the contrast margin is the gap between the top two projections. A change point is the time where $k_t^{*}$ switches and the combined margin is largest:

$$t^{*} \;=\; \operatorname*{arg\,max}_{t : \, k_t^{*} \neq k_{t+1}^{*}}\; \frac{p_t^{(k_t^{*})} - p_t^{(k_{t+1}^{*})} \;+\; p_{t+1}^{(k_{t+1}^{*})} - p_{t+1}^{(k_t^{*})}}{2}$$

5. Why this is the right detector for placer mining

A naive "before vs after" detector cannot distinguish a snowmelt edge from a dredge pond — both look like new water. Joint spectral-temporal PCA learns the full annual signature of every land-cover class in the AOI: snowmelt has a distinctive March–May peak followed by recovery; vegetation has a sigmoidal green-up; mining-disturbed ground has a sustained drop in NIR with a sustained rise in SWIR. These are all different PCs, and pixels that switch from one to another are exactly the events we are looking for.

6. Cloud cleanup is part of the model

Residual clouds break the temporal-PCA logic spectacularly: a single cirrus pixel that survives SCL masking adds a huge isolated outlier to a pixel's fingerprint, which the PCA then captures as a "principal mode" and the change-point detector flags as a class switch. We added three stages before interpolation:

  1. Dilate the per-frame any-band-NaN mask by 2 px and reapply.
  2. Drop entire scenes with < 85% valid AOI pixels after dilation.
  3. Reject per-pixel temporal outliers above 4 robust z-scores (median + MAD).

On the test_2020 cube this dropped 22 of 39 raw scenes but raised PC1's variance share from 30% to 73% and dropped the change-pixel count from 84,725 (35% of AOI) to 3,057 (1.26%). The surviving change-pixels are no longer cloud edges.

References