Matched-filter retrieval of methane column enhancement

Matched-filter retrieval on a synthetic plume scene¶

The matched filter is the maximum-likelihood detector for a known signature $t$ in additive Gaussian noise with known covariance $\Sigma$ . For a pixel spectrum $x$ , background mean $\mu$ , the abundance estimate is

\hat\varepsilon \;=\; \frac{(x - \mu)^\top \Sigma^{-1} t}{t^\top \Sigma^{-1} t}\,.

(1)

The target $t$ is the radiance (or normalised radiance) perturbation a reference plume enhancement $\Delta\mathrm{VMR}_\text{ref}$ would produce; scaling $\hat\varepsilon$ by $\Delta\mathrm{VMR}_\text{ref}$ recovers an estimate of each pixel’s true enhancement. This notebook runs the complete pipeline end-to-end on a hyperspectral synthetic scene built from the CH $_4$ LUT plus a Gaussian-plume column map (with OU-turbulence meander enabled) — showing:

Scene construction — plug a gauss_puff simulation (with new OU turbulence support) into the nB-LUT → dirty scene → add photon noise.
Target spectrum construction via target_spectrum_normalized_linear (the classical MF assumption).
Background statistics — build_lowrank_covariance_operator returns a structured gaussx.LowRankUpdate covariance plus its background mean in one trimmed pass. gaussx.solve applies the Woodbury identity internally so the matched filter never materialises an $O(B^3)$ inverse.
Retrieval — matched_filter_image_op → $\hat\varepsilon$ → $\Delta X$ map.
Validation — scatter vs. truth, RMSE, per-pixel detection SNR.

The retrieval is hyperspectral here (~500 channels over the $\nu_3$ Q-branch) because the matched filter’s SNR advantage scales as $\sqrt{n_\text{channels}}$ — a 2-band multispectral retrieval would need much stronger plumes to be worth the ceremony.

import warnings
from pathlib import Path


warnings.filterwarnings("ignore")

import jax
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

from plumax.gauss_puff import OUTurbulence, simulate_puff
from plumax.radtran import (
    ObservationGeometry,
    SpectralResponseFunction,
    build_lowrank_covariance_operator,
    build_nb_lut,
    inject_plume,
    matched_filter_image_op,
    matched_filter_snr_op,
    target_bands,
    target_spectrum_normalized_linear,
)


# The gaussx low-rank covariance path needs 64-bit JAX; enable it before
# building any arrays.
jax.config.update("jax_enable_x64", True)

# Built by 01_hapi_lut_ch4.ipynb into a per-user cache, outside the repository.
LUT_PATH = Path.home() / ".cache" / "plumax" / "hapi_lut" / "ch4_absorption_lut.nc"
assert LUT_PATH.exists(), f"Run 01_hapi_lut_ch4.ipynb first — missing {LUT_PATH}"

ch4_lut = xr.open_dataset(LUT_PATH)
nu_grid = ch4_lut["wavenumber"].values
print(f"LUT: {nu_grid.size} channels over {nu_grid.min():.0f}-{nu_grid.max():.0f} cm⁻¹")

LUT: 10000 channels over 4000-4500 cm⁻¹

1. Hyperspectral “instrument” and atmospheric state¶

We build a narrow-band SRF with one band per LUT wavenumber sample — effectively treating the LUT grid as the instrument’s native wavelength axis. This is the hyperspectral limit of the nB-LUT machinery and keeps the demo self-contained (no extra resampling).

geom = ObservationGeometry(sza_deg=30.0, vza_deg=0.0, path_length_cm=8e5)
T_K, p_atm = 260.0, 0.6
VMR_BG = 1.9e-6

# Subsample the LUT axis to a manageable hyperspectral channel count. The
# matched filter's SNR scales as √n_channels, so a few hundred channels is
# usually the sweet spot — more than that adds cost without useful signal.
channel_stride = max(1, nu_grid.size // 200)
nu_inst = nu_grid[::channel_stride]
wl_inst = 1e7 / nu_inst
sort = np.argsort(wl_inst)
srf = SpectralResponseFunction(
    wavelengths_hr_nm=wl_inst[sort],
    band_centers_nm=wl_inst[sort],
    band_widths_nm=np.full(wl_inst.size, 5.0),  # narrow, ~1 cm⁻¹ equivalent
    band_names=tuple(f"c{i:04d}" for i in range(wl_inst.size)),
    srf_type="gaussian",
)
print(
    f"Hyperspectral SRF: {srf.n_bands} channels (every {channel_stride}-th LUT sample)"
)

Hyperspectral SRF: 200 channels (every 50-th LUT sample)

2. Synthetic scene construction¶

Use gauss_puff.simulate_puff with OU turbulence enabled — this is the new feature from this PR. The OU process adds realistic meander on top of the deterministic wind advection, giving a plume with the kind of off-axis wiggles a matched filter has to handle in production.

time_array = np.linspace(0.0, 120.0, 13)
wind_speed = np.full_like(time_array, 3.0)
wind_direction = np.full_like(time_array, 270.0)  # wind from west

turb = OUTurbulence(sigma_fluctuations=0.4, correlation_time=30.0)

puff_ds = simulate_puff(
    emission_rate=5.0,  # kg/s ≈ 18 000 kg/hr — strong for clear retrieval signal
    source_location=(0.0, 0.0, 2.0),
    wind_speed=wind_speed,
    wind_direction=wind_direction,
    stability_class="D",
    domain_x=(-100.0, 600.0, 141),
    domain_y=(-150.0, 150.0, 61),
    domain_z=(0.0, 50.0, 11),
    time_array=time_array,
    release_frequency=1.0,
    turbulence=turb,
    turbulence_seed=0,
)

# Use the last time-step's column field (plume fully formed by then).
col_kg_m2 = puff_ds["column_concentration"].isel(time=-1).values  # (x, y)
M_CH4 = 16.04e-3
delta_X_map = (col_kg_m2 / M_CH4).T  # → (y, x) for image plots

print(
    f"ΔX map: {delta_X_map.shape}, peak {delta_X_map.max():.2f} mol/m², "
    f"nonzero pixels {np.sum(delta_X_map > 0.05)}"
)

ΔX map: (61, 141), peak 10.58 mol/m², nonzero pixels 1275

Build the nB LUT once, inject the plume, and add a touch of photon noise so the retrieval has something to fight against.

nb_lut = build_nb_lut(
    ch4_lut,
    srf,
    T_K=T_K,
    p_atm=p_atm,
    amf=geom.air_mass_factor,
    max_delta_column=20.0,
    n_grid=2001,
)

ny, nx = delta_X_map.shape
clean_scene = np.ones((srf.n_bands, ny, nx), dtype=float)
dirty_scene = inject_plume(clean_scene, delta_X_map, nb_lut)

rng = np.random.default_rng(0)
noise_std = 0.003  # ~0.3% fractional Gaussian noise per pixel per channel
noisy_scene = dirty_scene + rng.normal(scale=noise_std, size=dirty_scene.shape)

fig, axes = plt.subplots(1, 3, figsize=(12, 3.5))
axes[0].imshow(delta_X_map, origin="lower", cmap="viridis")
axes[0].set_title(r"truth $\Delta X$ [mol/m²]")
axes[1].imshow(
    dirty_scene[len(dirty_scene) // 2], origin="lower", cmap="gray", vmin=0.7, vmax=1.0
)
axes[1].set_title("dirty scene (mid channel)")
axes[2].imshow(
    noisy_scene[len(noisy_scene) // 2], origin="lower", cmap="gray", vmin=0.7, vmax=1.0
)
axes[2].set_title(f"+ noise (σ={noise_std:.3f})")
for ax in axes:
    ax.set_xticks([])
    ax.set_yticks([])
plt.tight_layout()
plt.show()

3. Target spectrum¶

Use the linearised target (Maclaurin order 1) at a reference enhancement of 1 ppm. The matched filter will return abundance $\hat\varepsilon$ in units of “multiples of the reference target”, which we map back to $\Delta\mathrm{VMR}$ and then to a column enhancement.

delta_ref = 1e-6  # 1 ppm reference
t_hr = target_spectrum_normalized_linear(
    ch4_lut,
    nu_inst,
    T_K=T_K,
    p_atm=p_atm,
    delta_vmr=delta_ref,
    path_length_cm=geom.path_length_cm,
    amf=geom.air_mass_factor,
)
t_b = target_bands(t_hr, srf, wl_inst)
print(
    f"Target spectrum: {t_b.size} bands, RMS magnitude {np.sqrt((t_b**2).mean()):.3e}"
)

fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(wl_inst[sort], t_b, lw=0.8)
ax.set_xlabel("Wavelength [nm]")
ax.set_ylabel(f"Target $t(\\lambda)$ @ $\\Delta$VMR={delta_ref * 1e6:.1f} ppm")
ax.set_title("Normalised-radiance target spectrum")
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()

Target spectrum: 200 bands, RMS magnitude 6.853e-02

4. Background statistics¶

Trimmed mean + low-rank covariance, both robust to the plume contaminating a small fraction of pixels. Default trim_frac=0.1 drops the top/bottom 10% of pixels per channel by brightness.

cov_op, mu = build_lowrank_covariance_operator(
    noisy_scene,
    rank=8,
    trim_frac=0.05,
    regularization=1e-8,
)
print(
    f"μ: {mu.shape}, "
    f"cov: LowRankUpdate of rank {int(cov_op.rank)} over {cov_op.in_size()} bands"
)

μ: (200,), cov: LowRankUpdate of rank 8 over 200 bands

5. Run the matched filter¶

from plumax.radtran.config import number_density_cm3


# Operator-backed retrieval: gaussx.solve(Σ, t) runs once via the Woodbury
# identity in O(B·rank) instead of an O(B³) dense inverse on Σ. For the
# ~200-channel hyperspectral setup below with rank=8 that's ~25× faster,
# and the savings compound when rank stays small as B grows.
eps_map = matched_filter_image_op(noisy_scene, mu, cov_op, t_b)
snr_map = matched_filter_snr_op(eps_map, cov_op, t_b)

# ε̂ is in units of ΔVMR_ref, so ε̂ · ΔVMR_ref = ΔVMR. To map to a column
# enhancement [mol/m²]:  ΔX = ΔVMR · N_total · L_vert · 1e4 / N_A.
N_total = number_density_cm3(p_atm, T_K)
N_A = 6.02214076e23
delta_X_map_retrieved = eps_map * delta_ref * N_total * geom.path_length_cm * 1e4 / N_A

fig, axes = plt.subplots(1, 3, figsize=(13, 3.5))
axes[0].imshow(delta_X_map, origin="lower", cmap="viridis")
axes[0].set_title(r"truth $\Delta X$")
axes[1].imshow(
    delta_X_map_retrieved,
    origin="lower",
    cmap="viridis",
    vmin=delta_X_map.min(),
    vmax=delta_X_map.max(),
)
axes[1].set_title(r"matched-filter $\hat{\Delta X}$")
im = axes[2].imshow(snr_map, origin="lower", cmap="magma")
axes[2].set_title("detection SNR")
fig.colorbar(im, ax=axes[2], fraction=0.046)
for ax in axes:
    ax.set_xticks([])
    ax.set_yticks([])
plt.tight_layout()
plt.show()

6. Validation¶

Scatter retrieved vs. truth for the plume pixels, compute RMSE, and compare against the no-retrieval baseline (constant $\hat{\Delta X} = 0$ ).

mask = delta_X_map > 0.01  # plume pixels
x = delta_X_map[mask]
y = delta_X_map_retrieved[mask]
rmse = float(np.sqrt(np.mean((y - x) ** 2)))
rmse_zero = float(np.sqrt(np.mean(x**2)))
print(f"MF RMSE     = {rmse:.3f} mol/m²")
print(f"zero RMSE   = {rmse_zero:.3f} mol/m² (no-retrieval baseline)")
print(f"improvement = {rmse_zero / rmse:.1f}×")

fig, ax = plt.subplots(figsize=(5, 5))
ax.scatter(x, y, s=6, alpha=0.4, color="C0")
lim = [0, max(x.max(), y.max()) * 1.05]
ax.plot(lim, lim, "k--", lw=1, label="1:1")
ax.set_xlim(lim)
ax.set_ylim(lim)
ax.set_xlabel(r"truth $\Delta X$ [mol/m²]")
ax.set_ylabel(r"retrieved $\hat{\Delta X}$ [mol/m²]")
ax.set_title(f"RMSE = {rmse:.3f} mol/m²")
ax.legend()
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()

MF RMSE     = 0.903 mol/m²
zero RMSE   = 1.521 mol/m² (no-retrieval baseline)
improvement = 1.7×

What this notebook demonstrates¶

The radtran stack — SRF → forward-model → nB-LUT → matched filter — cleanly composes to an end-to-end synthetic observation-plus-retrieval loop.
The matched filter tracks the plume column map with an RMSE of under 10% of the peak column (≈0.9 vs ≈10.6 mol/m² here, a 1.7× improvement over the no-retrieval baseline) at 0.3% fractional noise and a few hundred channels — even in the presence of OU-turbulence meander.
Target-spectrum construction, background estimation, and retrieval are all driven by a single HAPI-LUT dataset. Swapping the absorber (CO $_2$ , H $_2$ O, N $_2$ O, …) is a one-line change.
The covariance lives as a structured gaussx.LowRankUpdate throughout, and the matched filter calls gaussx.solve once to produce $\Sigma^{-1}\,t$ via Woodbury. No dense inverse is ever materialised; the hot path stays $O(B\cdot k)$ even as the channel count grows.

Next natural extensions: replace the flat “clean scene” with a real Sentinel-2 or EMIT radiance cube, add the PSF / GSD observation operators from jej_vc_snippets/methane_retrieval/lut_obs_op.py, and swap the linearised target for target_spectrum_normalized_nonlinear when analysing strong plumes.