Mathematical description of the steps that turn 539 delivered Sentinel-1 frames into the trend, hotspot, and projection maps in the viewer. Equations are written out so each step is auditable. Limitations are listed inline rather than relegated to a footnote.
The pipeline operates on 539 pre-processed Sentinel-1 IW GRD HD acquisitions delivered by Aramco over an AOI in eastern Saudi Arabia (~12.6 × 30.6 km, centred near the Ghawar oilfield). Each acquisition arrives as two 8-bit PNGs (Area 1 north, Area 2 south) that have been clipped to the AOI and amplitude-stretched for visual display.
| Symbol | Meaning |
|---|---|
| $I_t(x,y)$ | Delivered amplitude in DN at pixel $(x,y)$ and acquisition $t$ (range 0–255). |
| $N$ | Number of acquisitions retained after orbit filtering ($N=271$). |
| $t \in [0, T]$ | Time in years since the first retained acquisition. $T \approx 10$. |
| $\sigma^0$ | The calibrated radar backscatter coefficient (dB), not what we have. |
| $L$ | Equivalent number of looks (multi-looking factor in SAR processing). |
The delivered frames are pre-clipped to a fixed lat/lon bounding box, but small geocoding residuals leave the same ground feature at slightly different pixel positions across acquisitions. We align each frame to a reference using normalised cross-correlation (NCC) on a stable template window.
For a template $\mathbf{T}$ taken from the reference frame and a search window $\mathbf{W}_i$ from candidate frame $i$, the NCC at shift $(\Delta x, \Delta y)$ is
\[ \rho_i(\Delta x, \Delta y) \;=\; \frac{\sum_{p}\bigl(T(p) - \bar T\bigr)\bigl(W_i(p + \Delta) - \bar W_i\bigr)} {\sqrt{\sum_p (T(p) - \bar T)^2 \cdot \sum_p (W_i(p+\Delta) - \bar W_i)^2}} \]and the optimal shift is
\[ (\widehat{\Delta x}_i, \widehat{\Delta y}_i) = \mathop{\mathrm{arg\,max}}_{\Delta x, \Delta y \in [-35, 35]} \rho_i(\Delta x, \Delta y). \]We compute shifts independently for Area 1 and Area 2 because they were delivered as separately-clipped tiles and can have slightly different residual offsets. After collecting the per-frame shifts we median-centre them: the per-area median shift across all frames is subtracted from every frame's shift, so the typical frame ends up at $(0,0)$ and only outlier frames pay the warp-induced edge loss.
The frame is then warped with
\[ I_t'(x,y) = I_t\bigl(x - \widehat{\Delta x}_t,\; y - \widehat{\Delta y}_t\bigr), \]using bilinear interpolation and zero-fill at boundaries (OpenCV
warpAffine with BORDER_CONSTANT). Residual jitter
measured by phase correlation on the aligned outputs is
$\le 2$ px in $x$ and $\le 1$ px in $y$ across all 271 retained frames.
The 539 delivered frames span two Sentinel-1 relative orbits with different ground tracks (acquisition times ~14:48 and ~14:57 UTC). Their swath edges fall in opposite corners of the AOI, so mixing them produces a wedge-shaped artefact in any cross-temporal analysis. We retain only the 14:47–14:50 UTC cluster ($N = 271$ frames).
Within a single orbit the swath wedge still drifts slowly with orbital maintenance. Define the per-pixel coverage count over the retained frames as
\[ \mathcal{C}(x,y) \;=\; \sum_{t=1}^{N} \mathbf{1}\!\bigl[I_t'(x,y) > \tau_{\mathrm{OOS}}\bigr], \quad \tau_{\mathrm{OOS}} = 3, \]and the common-coverage mask as
\[ M_{\mathrm{cov}}(x,y) \;=\; \mathbf{1}\!\bigl[\mathcal{C}(x,y) \ge 0.95 \, N\bigr]. \]Every output frame is multiplied by $M_{\mathrm{cov}}$, so all 271 frames display the same active region (the rest is a constant black border). $M_{\mathrm{cov}}$ retains 29.8% of the 1262 × 3066 mosaic.
SAR amplitude is contaminated by multiplicative speckle: if $\sigma^0(x,y)$ is the true backscatter and $n_t$ a unit-mean speckle realisation, then $|S_t(x,y)|^2 \sim \sigma^0(x,y) \cdot n_t$ where for L-look intensity images $L \cdot n_t \sim \chi^2_{2L}$. Speckle has no temporal coherence — across acquisitions $n_t$ realisations are statistically independent.
We exploit that independence by per-pixel temporal multi-look with a rolling median:
\[ \tilde I_t(x,y) \;=\; \operatorname{median}\bigl\{I'_{t-2}(x,y),\;\ldots,\;I'_{t+2}(x,y)\bigr\} \quad (\text{window } w = 5). \]Median is preferred over mean because it is robust to outliers (single noisy frames, residual mis-registration). Temporal multi-looking suppresses multiplicative noise without smearing spatial edges — the failure mode of the spatial Gaussian smoothing that this pipeline used in an earlier iteration.
For each pixel we fit the simple linear model
\[ \tilde I_t(x,y) \;=\; a(x,y) \,+\, b(x,y)\,t \,+\, \varepsilon_t(x,y), \]by ordinary least squares. With $\bar t = \frac{1}{N}\sum_t t$ and $S_{tt} = \sum_t (t - \bar t)^2$, the slope estimator is
\[ \widehat b(x,y) \;=\; \frac{\sum_{t}(t - \bar t)\bigl(\tilde I_t(x,y) - \bar I(x,y)\bigr)}{S_{tt}}, \quad \widehat a(x,y) = \bar I(x,y) - \widehat b(x,y)\,\bar t. \]Negative $\widehat b$ means the surface darkened over the decade — the smoothing signature consistent with sand burial at C-band.
Weekly SAR samples are not independent: surface moisture has memory, speckle realisations are slightly correlated at consecutive sub-pixel positions of the satellite, and the temporal multi-look itself introduces some autocorrelation. The OLS iid standard error of the slope
\[ \mathrm{SE}_{\mathrm{iid}}(\widehat b) \;=\; \sqrt{\dfrac{s_e^2}{S_{tt}}}, \qquad s_e^2 = \tfrac{1}{N-2}\sum_t \bigl(\tilde I_t - \widehat a - \widehat b\,t\bigr)^2 \]therefore understates the true uncertainty. We measure the per-pixel lag-1 autocorrelation of the OLS residuals,
\[ \hat\rho(x,y) \;=\; \frac{\sum_{t=2}^{N} e_t \cdot e_{t-1}}{\sum_{t=1}^{N} e_t^2}, \quad e_t = \tilde I_t - \widehat a - \widehat b\, t, \]and apply the classical effective-sample-size correction (Bayley & Hammersley 1946; see e.g. Kutner et al. 2004):
\[ n_{\mathrm{eff}} \;=\; N \cdot \frac{1 - \hat\rho}{1 + \hat\rho}, \quad \hat\rho \text{ clipped to } [0,\, 0.95]. \]The adjusted slope SE and t-statistic are then
\[ \mathrm{SE}_{\mathrm{adj}}(\widehat b) \;=\; \mathrm{SE}_{\mathrm{iid}}(\widehat b)\,\sqrt{N / n_{\mathrm{eff}}}, \qquad t(x,y) \;=\; \widehat b(x,y) \;/\; \mathrm{SE}_{\mathrm{adj}}(\widehat b). \]The two-sided p-value comes from Student's t-distribution with $n_{\mathrm{eff}} - 2$ degrees of freedom:
\[ p(x,y) \;=\; 2\bigl[1 - F_t\bigl(|t(x,y)|;\; n_{\mathrm{eff}} - 2\bigr)\bigr]. \]The viewer's AOI contains roughly $1.15 \times 10^6$ pixels inside the common-coverage mask. With that many simultaneous tests, uncorrected $p$-values are useless — at $\alpha = 0.05$ alone we would expect tens of thousands of false discoveries. We control the false discovery rate (FDR) using the Benjamini–Hochberg procedure (Benjamini & Hochberg 1995):
At our $\alpha = 0.05$ the empirical cutoff is $p^\star \approx 1.1 \times 10^{-3}$. The significance mask is
\[ M_{\mathrm{sig}}(x,y) \;=\; \mathbf{1}\!\bigl[p(x,y) \le p^\star\bigr] \cdot M_{\mathrm{cov}}(x,y). \]About 2.2% of in-coverage pixels survive this test. Of those, roughly half are darkening ($\widehat b < 0$) and half are brightening — the symmetry is expected under a null where most surfaces are stable.
An FDR-significant pixel represents a $10 \text{m} \times 10 \text{m}$ patch of ground. Rendered at AOI scale that is sub-pixel — invisible. To visualise the spatial concentration of sand-burial signal we apply a 2-D Gaussian kernel density estimate (Silverman 1986). Each significant pixel emits a Gaussian "puff" weighted by $|\widehat b|$.
Let $\mathcal{P}_- = \{(x_i, y_i) : M_{\mathrm{sig}}(x_i, y_i) = 1,\; \widehat b(x_i, y_i) < 0\}$ be the set of significantly-darkening pixels. The burial-density field is
\[ D_-(x,y) \;=\; \sum_{(x_i, y_i) \in \mathcal{P}_-} \bigl|\widehat b(x_i, y_i)\bigr| \cdot K_\sigma(x - x_i,\, y - y_i), \]with the standard isotropic 2-D Gaussian kernel
\[ K_\sigma(u, v) \;=\; \frac{1}{2\pi\sigma^2}\, \exp\!\left(-\frac{u^2 + v^2}{2\sigma^2}\right). \]An analogous $D_+$ is defined over significantly-brightening pixels. We use $\sigma = 30$ px $\approx 300$ m on the ground, chosen to match the natural patch scale of mixed barchan + sand-sheet fields on the Ghawar fringe (Yang et al. 2012; sheet morphologies typically 100–500 m across at this latitude).
In practice $D_\pm$ is computed by convolving the $|\widehat b| \cdot M_{\mathrm{sig}} \cdot \mathbf{1}[\widehat b \lessgtr 0]$ field with the same Gaussian kernel:
\[ D_\pm \;=\; \bigl(\,|\widehat b|\, M_{\mathrm{sig}}\, \mathbf{1}[\,\widehat b \lessgtr 0\,]\,\bigr) \,\ast\, K_\sigma, \]which is implemented via two separable 1-D convolutions
(cv2.GaussianBlur), $\mathcal{O}(HW\sigma)$ per pass.
For rendering we normalise each density field by its own 99th percentile inside the coverage mask, suppress pixels below 5% relative density (the KDE tails), and apply a hysteresis rule $D_- > 1.1\, D_+ \Rightarrow$ red, $D_+ > 1.1\, D_- \Rightarrow$ blue, otherwise transparent.
The KDE in §8 conveys spatial structure but its values are not slopes; they have units of "DN/yr times pixels". For the projection tab we need a slope at every pixel near a burial blob. We use the Nadaraya–Watson kernel regression estimator (Nadaraya 1964; Watson 1964):
\[ \tilde b(x,y) \;=\; \frac{\sum_{(x_i, y_i)\in\mathcal{P}_{\mathrm{sig}}} \widehat b(x_i, y_i)\,K_\sigma(x - x_i,\, y - y_i)} {\sum_{(x_i, y_i)\in\mathcal{P}_{\mathrm{sig}}} K_\sigma(x - x_i,\, y - y_i)}, \]i.e. a kernel-weighted local average of FDR-significant slopes, evaluated everywhere. Equivalently and computationally efficiently:
\[ \tilde b \;=\; \frac{\bigl(\widehat b\, M_{\mathrm{sig}}\bigr) \,\ast\, K_\sigma} {M_{\mathrm{sig}} \,\ast\, K_\sigma + \epsilon}, \]where the denominator $\mathcal{M} = M_{\mathrm{sig}} * K_\sigma$ is the local kernel mass on significant pixels. Unlike a raw Gaussian blur of $\widehat b$, the Nadaraya–Watson form preserves the physical units of $\widehat b$ (DN/yr): in regions densely covered by significant pixels $\tilde b$ approximates the local mean slope, while sparsely-covered neighbourhoods get a slope only if a few significant pixels still contribute meaningful mass.
We zero out $\tilde b$ where $\mathcal{M} < 0.01$ — a 1% local-mass threshold prevents kernel tails from extending burial predictions into regions that contain no FDR-significant evidence.
The viewer's Projection tab answers "by year $N$, how much amplitude darkening has the regional trend accumulated?". Under the assumption that the linear trend continues, the expected accumulated darkening at pixel $(x,y)$ after $N$ years is
\[ \Delta(x, y, N) \;=\; -\,\tilde b(x,y) \cdot N \quad \text{if } \tilde b(x,y) < -\, b_{\min} = -0.3 \text{ DN/yr}, \]and undefined otherwise. The viewer renders $\Delta$ as an alpha-modulated heat ramp: each colored pixel is mixed with the dimmed grayscale baseline by an opacity
\[ \alpha(\Delta) \;=\; 0.08 \,+\, 0.85 \cdot \mathrm{clip}\!\left(\frac{\Delta - \Delta_{\min}}{\Delta_{\max} - \Delta_{\min}},\; 0,\; 1\right), \]with $\Delta_{\min} = 2$ DN (just above speckle noise) and $\Delta_{\max} = 30$ DN (saturation). The hue ramps from yellow at $\Delta = \Delta_{\min}$ through orange to red at $\Delta \ge \Delta_{\max}$.
Independent of the KDE visualisation, we extract a discrete list of burial hotspots by 8-connected component labelling on a thresholded version of the un-smoothed significance map:
\[ H(x,y) \;=\; \mathbf{1}\!\Bigl[ M_{\mathrm{sig}}(x,y) = 1 \;\land\; \widehat b(x,y) < -1.0 \text{ DN/yr} \Bigr]. \]We compute connected components on $H$ using OpenCV's
connectedComponents with 8-neighbour connectivity. Components
with area $\ge 50$ px (= 5,000 m²) are retained and ranked by
the product of cluster area and absolute mean slope. The top 20 (typically fewer, since few clusters survive the FDR-significance pre-filter) are returned to the viewer as click-to-zoom targets.
| Limitation | Where it bites | What would fix it |
|---|---|---|
| 8-bit amplitude-stretched input | Slope units (DN/yr) are arbitrary and possibly time-dependent. Burial threshold of 25 DN is a visual guess. | Re-pull original Level-1 GRDs from Copernicus; calibrate to $\sigma^0$ in dB; redefine threshold by an asymptotic fit to known sand-dune reference pixels. |
| Linear OLS trend model | Sand burial is a step function (regime change) once amplitude approaches the speckle floor; OLS underestimates the rate of sudden burial and over-predicts continued darkening after a pixel has already saturated. | Replace with BFAST / CUSUM change-point detection, or fit a saturating exponential $\sigma^0 = A + (B - A) e^{-kt}$ with $k$ as the rate parameter. |
| Per-pixel independence | Cannot predict burial of currently-stable pixels (kinematic critique: an advancing dune face has slope $= 0$ at its leading edge until it actually arrives). | Extract the binary "sand-like" mask per date, trace dune-front boundaries, compute migration vectors via optical flow, project the boundary forward as a kinematic problem. |
| 50% of frames dropped | Lost statistical power and temporal sampling for event detection. | Re-integrate both relative orbits with per-orbit calibration and incidence-angle covariates after returning to calibrated GRDs. |
| No external validation | "SAR darkening = sand burial" is unverified. Could be vegetation loss, surface paving, moisture change. | Cross-check with Sentinel-2 NDVI on hotspots; pull ALOS-2 L-band coherence stack for buried-infrastructure detection; ground-truth a sample of hotspots in the next field visit. |
| No uncertainty propagation | Projection presents accumulated darkening as a point estimate; the actual confidence interval is wide for low-$|b|$ pixels. | Bootstrap or analytical SE of $\tilde b$ → confidence band on $\Delta(x,y,N)$; render only pixels whose lower-CI burial estimate exceeds the horizon. |