Principles of RRINGG

Overview

RRINGG provides two main functionalities for combining GNSS and InSAR data:

  1. Correction of the InSAR velocity field (-po correction)

  2. Referencing of InSAR time series to a global geodetic frame (-po referencing)

Both rely on fitting a bilinear surface (affine plane) to the differences between GNSS and InSAR observations. The GNSS network supplies the geodetic truth, ensuring that InSAR products are both bias-free and compatible with international reference frames (ITRF2014, ITRF2020).

Part I — Correction of the InSAR velocity field

Objective

Remove systematic long-wavelength biases from an InSAR mean velocity field using co-located GNSS stations as absolute reference. Typical bias sources include orbital ramps, ionospheric gradients, and reference-frame offsets.

Mathematical model

The bias between InSAR and GNSS LOS velocities is modelled as a plane that varies linearly with geographic position:

\[\Delta v(x, y) = a \cdot x + b \cdot y + c\]

where \((x, y)\) are longitude and latitude (decimal degrees) and \(a, b, c\) are the plane parameters estimated by fitting the InSAR–GNSS velocity differences at \(N\) station locations.

The corrected velocity field is then:

\[v_{\text{corr}}(x, y) = v_{\text{InSAR}}(x, y) - \Delta v(x, y)\]

Processing steps

  1. Retrieve GNSS station velocities (East, North, Up) for all stations inside or near the InSAR footprint.

  2. Project each GNSS velocity vector onto the InSAR line of sight using the per-pixel LOS unit vectors.

  3. Sample the InSAR velocity at each station location (windowed average over a configurable neighbourhood).

  4. Compute the difference \(\Delta v_i = v_{\text{InSAR},i} - v_{\text{GNSS,LOS},i}\) at each station.

  5. Fit the plane \((a, b, c)\) to these differences using the configured fitting algorithm (IRLS or RANSAC).

  6. Subtract the fitted plane pixel-wise from the full raster.

Part II — Referencing of InSAR time series

Objective

Convert relative InSAR displacement time series into absolute displacements expressed in a global geodetic frame, using GNSS trajectory models as the absolute reference.

Two methods are available, selected with --referencing_method.

Method 1 — Mean velocity referencing (mean_velocity)

Principle

GNSS motion is assumed to be linear over the observation period. The GNSS LOS velocity at each station is integrated to predict the absolute displacement at each InSAR acquisition epoch \(t_i\):

\[G(x, y, t_i) = v_{\text{GNSS,LOS}}(x, y) \cdot (t_i - t_{\text{ref}})\]

where \(t_{\text{ref}}\) is the reference date (usually the first acquisition). The bias at each epoch is modelled as a plane:

\[B(x, y, t_i) = D_{\text{InSAR,rel}}(x, y, t_i) - G(x, y, t_i) = a_i \cdot x + b_i \cdot y + c_i\]

The absolute displacement is:

\[D_{\text{abs}}(x, y, t_i) = D_{\text{InSAR,rel}}(x, y, t_i) - B(x, y, t_i)\]

Best for: tectonic drift, steady subsidence — environments where deformation is essentially linear.

Limitation: cannot capture non-linear transients (post-seismic relaxation, volcanic inflation, seasonal cycles).

Method 2 — Temporal differential referencing (differential)

Principle

Uses the full GNSS position time series (daily or weekly sampling). A trajectory model (linear trend + annual + semi-annual sinusoids) is fitted to the GNSS time series and evaluated at each InSAR acquisition date to predict the absolute displacement:

\[G^*(x, y, t_i) = G(x, y, t_i) - G(x, y, t_0)\]

where \(t_0\) is the global reference epoch (earliest InSAR acquisition). The bias plane at each epoch is:

\[B(x, y, t_i) = D_{\text{InSAR,rel}}(x, y, t_i) - G^*(x, y, t_i) = a_i \cdot x + b_i \cdot y + c_i\]

Best for: dynamic environments — earthquakes, volcanic unrest, hydrological loading, postseismic deformation.

Limitation: requires dense GNSS time series (daily sampling preferred); sensitive to gaps and offsets in the GNSS record.

Method comparison

Feature

Mean velocity

Differential

GNSS input

Mean velocities only

Full position time series

Deformation model

Linear only

Linear + non-linear

GNSS sampling requirement

Sparse is sufficient

Dense required (daily preferred)

Computational cost

Low

Higher (one plane fit per epoch)

Best for

Tectonic drift, subsidence

Earthquakes, volcanoes, seasonal cycles

Limitation

Misses transients

Sensitive to GNSS noise and gaps

Part III — Station weighting

A priori weights (Fisher-information)

Not all GNSS stations are equally reliable. RRINGG assigns each station an a priori weight based on two independent quality indicators:

Precision — inversely proportional to the squared projected LOS velocity uncertainty \(\sigma_i\):

\[w_{\text{precision},i} \propto \frac{1}{\sigma_i^2}\]

Temporal coverage — proportional to the number of GNSS observation days \(N_i\), capped at twice the minimum required:

\[w_{\text{coverage},i} = \min\!\left(\frac{N_i}{2 \cdot N_{\min}},\; 1\right)\]

The combined a priori weight is the element-wise product:

\[w_{\text{prior},i} = w_{\text{precision},i} \times w_{\text{coverage},i} \;\propto\; \frac{N_i}{\sigma_i^2}\]

This is the Fisher-information formulation: it automatically penalises stations with short or imprecise time series, without requiring manual tuning of per-criterion coefficients.

Weights are normalised to unit mean before entering the plane fit. The --fitting weighted flag (default: on) activates this scheme; setting it to False in config.py applies uniform weights instead.

Robust (IRLS) weights

During iteratively reweighted least squares (see below), each station also receives a robust weight derived from its normalised residual \(r_i = e_i / \hat{\sigma}\), where \(e_i\) is the fitting residual and \(\hat{\sigma}\) is the MAD-based scale estimate.

The final weight used in the normal equations is the product:

\[w_i = w_{\text{prior},i} \times w_{\text{robust},i}\]

Part IV — Plane fitting algorithms

IRLS — Iteratively Reweighted Least Squares (least_squares)

IRLS solves the weighted normal equations

\[\hat{\boldsymbol{\beta}} = (\mathbf{X}^\top \mathbf{W} \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{W} \mathbf{y}\]

at each iteration, updating the diagonal weight matrix \(\mathbf{W}\) from the current residuals using a robust influence function.

Two influence functions are available (CFG.fitting.weight_function):

Huber (default, c = 1.345):

\[\begin{split}w_{\text{robust}}(r) = \begin{cases} 1 & \text{if } |r| \le c \\ c / |r| & \text{if } |r| > c \end{cases}\end{split}\]

Preserves ordinary least-squares behaviour for small residuals; linearly downweights large ones. Breakdown point ≈ 50 %, efficiency ≈ 95 % under normality.

Tukey bisquare (c = 4.685):

\[\begin{split}w_{\text{robust}}(r) = \begin{cases} \left(1 - (r/c)^2\right)^2 & \text{if } |r| \le c \\ 0 & \text{if } |r| > c \end{cases}\end{split}\]

Hard-zeros stations beyond the threshold — more aggressive outlier rejection. Recommended only when gross outliers are expected.

Convergence — iterations stop when \(\max_j |\Delta\hat{\beta}_j| / (|\hat{\beta}_j| + \varepsilon) < \texttt{tolerance}\) or after CFG.fitting.max_iter iterations (default 50).

After convergence, a single k-fold cross-validation (CFG.fitting.cv_folds, default 5) is performed on the final model to compute CV-MSE.

RANSAC (ransac)

Random Sample Consensus draws repeated minimal subsets of stations, fits a plane to each, and counts inliers (stations whose residual falls within a threshold). The iteration with the most inliers wins; a final weighted OLS re-fit is performed on those inliers.

RANSAC is more aggressive than IRLS: stations beyond the inlier threshold receive exactly zero weight (hard rejection vs. continuous downweighting).

When to use RANSAC: dense networks with heterogeneous quality where a significant fraction of stations may be contaminated (e.g., the Italian NGL network in Case 1 of the testing guide).

Limitation: requires enough stations relative to the outlier fraction. If the detected outlier fraction exceeds CFG.fitting.ransac_max_outlier_fraction (default 30 %), the fit is flagged as unreliable. Use least_squares when the station count is low (< 10).

Part V — Spatial distribution and Kriging fallback

Voronoi distribution check

Before fitting, RRINGG assesses whether the GNSS station network provides adequate spatial coverage of the InSAR footprint. It computes a Voronoi diagram of the station positions (projected to metric coordinates), intersects each cell with the InSAR footprint, and checks the coefficient of variation of the cell areas:

\[\text{CV} = \frac{\text{std}(\text{areas})}{\text{mean}(\text{areas})}\]

If \(\text{CV} > \text{CFG.spatial.distribution\_cv\_threshold}\) (default 0.6) and a small fraction of stations dominates a disproportionate share of the footprint area, the distribution is flagged as inhomogeneous.

Ordinary Kriging fallback

When the Voronoi check fails, RRINGG generates a set of synthetic control points on a regular grid (spacing CFG.spatial.grid_step_m, default 120 km) using Ordinary Kriging interpolated from the available GNSS stations. These synthetic points supplement the real stations in the plane fit, providing additional spatial constraints in under-sampled regions.

Kriging weights are derived from the interpolation variance \(\sigma^2_{\text{kriging}}\):

\[w_{\text{kriging},j} \propto \frac{1}{\sigma^2_{\text{kriging},j}}\]

and normalised to the range [1, 10]. A zero variance is clamped to \(10^{-6}\) to avoid division by zero.