11.27.5. Profile-Level Anomaly Detection#

Machine-Learning Quality Scoring and Quality-Control Confidence Scoring both score data quality at the (station, frequency) grain – one number per cell. That grain is right for deciding whether to mask a single bad frequency, but it cannot flag a station whose entire spectrum is unusual in shape or level even though every individual frequency looks locally reasonable. AnomalyDetector works at the opposite grain: one feature vector per station, built by flattening prepare_z_features()’s (n_components, n_freqs) array, and one anomaly score per station.

The detector is a fully connected autoencoder – the same reconstruction-error idea as Learned Spectral Denoising, but on whole flattened station vectors instead of per-channel sequences, and used as a diagnostic score rather than to produce a cleaned output:

(1)#\[\hat{\mathbf{x}} = g(f(\mathbf{x})), \qquad s_i = \frac{\|\mathbf{x}_i - \hat{\mathbf{x}}_i\|_2^2}{d}\]

where \(f : \mathbb{R}^d \to \mathbb{R}^k\) is the encoder (\(k \ll d\), set by latent_dim) and \(g\) the decoder. A station is flagged when its reconstruction error \(s_i\) exceeds the threshold_percentile-th percentile (default the 95th) of the training scores – so, by construction, roughly 5% of a representative training set is expected to be flagged. When neither PyTorch nor TensorFlow is installed, the same interface falls back to PCA reconstruction via sklearn.decomposition.PCA.

11.27.5.1. Fitting and scoring#

>>> from pycsamt.ai.processing import AnomalyDetector, prepare_z_features
>>> X = prepare_z_features(sites, n_components=4)
>>> n_sites, n_comp, n_freq = X.shape
>>> X_flat = X.reshape(n_sites, n_comp * n_freq)
>>> X_flat.shape
(28, 212)
>>> det = AnomalyDetector(latent_dim=8, channels=(32, 16))
>>> det.fit(X_flat, epochs=150, seed=0, verbose=False)
AnomalyDetector(n_features=212, latent_dim=8, torch, fitted)
>>> det.threshold_
2.4401
>>> flags = det.flag_anomalies(X_flat)
>>> int(flags.sum())
2

With latent_dim=8 and 212 raw features, the bottleneck compresses each station to 4% of its original size – tight enough that a station whose spectrum shape genuinely differs from the rest of the line cannot be reconstructed as well as a typical one.

../../_images/anomaly_summary.png

plot_anomaly_summary() – per-station scores (a) and their pooled distribution (b), with the fitted threshold shared between both panels.#

Two stations cross the threshold, both by a wide margin over the rest of the line: 18-024U and 18-025A – the last two, adjacent stations on the profile. 18-024U is one of the seven stations whose Machine-Learning Quality Scoring score collapsed under the hard Swift-skew rule – a station already flagged as electrically unusual by a completely different method agreeing with the anomaly detector here. Both flagged stations also share a trait neither the rule-based QC score nor a single-channel view would surface directly: their \(\phi_{xy}\) curves are noticeably noisier than the rest of the line (standard deviation \(\approx 28^\circ\) and \(25^\circ\), against a \(15^\circ\) median across all 28 stations) – a station-wide, cross-channel pattern rather than a handful of bad frequencies, exactly the kind of anomaly a whole-spectrum reconstruction-error model is built to catch and a per-frequency QC rule, which only ever looks at one station at a time against a fixed threshold, has no mechanism to compare against its neighbours to find.

../../_images/anomaly_training.png

plot_training_history() on the fitted detector.#

Warning

With only 28 training stations (about 25 after the validation split), this training curve is a textbook small-sample case: training loss keeps falling while validation loss bottoms out early (here, around epoch 7) and drifts back up afterward as the network starts memorizing individual training stations rather than the shared structure across them. fit() already restores the lowest-validation-loss weights automatically, so the deployed model is unaffected – but the curve itself is a reminder that a profile this short gives the Isolation Forest fallback and this network comparably little to generalize from. Prefer more stations, a smaller latent_dim, or the threshold_percentile rule-of-thumb over a literal reading of the loss curve when the training set is this small.

11.27.5.2. Sites in, filtered sites out#

apply() closes the same sites-in / sites-out loop as Learned Spectral Denoising’s apply() and Machine-Learning Quality Scoring’s apply(), but its shape is different: rather than mutating Z.z in place, it drops flagged stations outright via Sites.select(), since the detector scores whole stations, not station-frequency cells. There is therefore no inplace option here – Sites.select() always returns a new container, and so does this method.

>>> clean_sites = det.apply(sites, n_components=4)
>>> len(list(sites)), len(list(clean_sites))
(28, 26)
>>> {s.name for s in sites} - {s.name for s in clean_sites}
{'18-024U', '18-025A'}

n_components (and freq_ref/log_amp when set explicitly) must reproduce the exact feature layout the detector was fitted on – apply() checks the flattened length against the fitted n_features and raises ValueError rather than silently scoring a mismatched vector when they disagree.

Note

The two dropped stations above match the ones read off the summary figure earlier on this page – reran here through a different code path for illustration, not a different result. Because the network’s weight initialization is not seeded (see the warning in Learned Spectral Denoising about fit(seed=...) covering only the data split and synthetic noise, not PyTorch’s own initialization RNG), an independent re-fit of the detector can occasionally flag a different pair of borderline stations even with the same data and threshold_percentile.

11.27.5.3. Parameters#

latent_dim (default 32) is the main capacity knob: too large and the network can reconstruct almost anything, including real anomalies, driving every score toward zero; too small and it cannot represent normal variability either, inflating scores across the board. channels sets the encoder hidden widths (decoder mirrors them). threshold_percentile (default 95.0) trades false positives against false negatives directly – lower it to flag more stations for review, raise it to flag only the most extreme outliers.

Full signatures: AnomalyDetector, plot_anomaly_scores(), plot_anomaly_score_distribution(), and plot_anomaly_summary().

View the complete AnomalyDetector example sourceClick to inspect and copy the complete code
 1def run_anomaly(sites, save_dir: Path) -> None:
 2    print("\n=== AnomalyDetector ===")
 3    X = prepare_z_features(sites, n_components=4)
 4    labels = station_names(sites)
 5    n_sites, n_comp, n_freq = X.shape
 6    X_flat = X.reshape(n_sites, n_comp * n_freq)
 7    X_flat = np.nan_to_num(X_flat, nan=0.0)
 8    print(f"  feature matrix: {X_flat.shape}")
 9
10    det = AnomalyDetector(latent_dim=8, channels=(32, 16))
11    det.fit(X_flat, epochs=150, seed=0, verbose=False)
12    scores = det.transform(X_flat)
13    flags = det.flag_anomalies(X_flat)
14    print(f"  {det!r}")
15    print(f"  threshold_ (95th pct): {det.threshold_:.4f}")
16    print(f"  flagged anomalous: {int(flags.sum())} / {n_sites}")
17    if flags.any():
18        flagged = [labels[i] for i in np.where(flags)[0]]
19        print(f"  stations: {flagged}")
20
21    fig = plot_anomaly_summary(
22        scores, station_labels=labels, threshold=det.threshold_,
23        suptitle="AnomalyDetector -- L18PLT profile-level scores",
24    )
25    _save(fig, save_dir, "anomaly_summary.png")
26
27    ax = plot_training_history(det, title="AnomalyDetector training")
28    _save(ax.figure, save_dir, "anomaly_training.png")