Spectral matching workflow: selection, shaping, and suite scaling
In gmsp: Ground Motion Signal Processing

knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      eval = FALSE, purl = FALSE)

This vignette documents the complete spectral-matching workflow: how a suite of recorded ground motions is selected from a processed database, adjusted so that its response spectra match a target uniform-hazard spectrum (UHS), and delivered as traceable file products. Two runners drive the workflow from one JSON file per site:

library(gmsp)

runStage0(file = "gmsp/match.json", root = "/path/to/project")
runMatch(file = "gmsp/match.json", root = "/path/to/project")

Method overview

The workflow has three stages.

Spectral preselection (runStage0()). Every candidate record in the processed database pool is scored by the compatibility of its PSA spectrum (H1 component) against the target band defined by the UHS mean and its lower and upper envelope quantiles. The N most compatible records form the suite. Optional seismological screens (magnitude, distance, provider, one record per event) restrict the pool before ranking.

Modal shaping (per record). Each selected record is decomposed into intrinsic mode functions (see TS2IMF()), the retained modes are re-weighted by positive modal factors (shape.modalFactor) to improve spectral compatibility, and the reconstruction is normalized so the quadratic acceleration norm of the original record is preserved. The in-memory algorithm is exported as fitModalFactor().

Record scaling (suite). One positive coefficient $c_i$ per record (shape.recordFactor) scales the normalized reconstructed history. The in-memory algorithm is exported as fitRecordFactor(). The coefficients solve

$$ \min_{c_i > 0} \; \frac{1}{J}\sum_{j=1}^{J} w_j\, r_j^{2} \;+\; \frac{\lambda_B}{N J}\sum_{i=1}^{N}\sum_{j=1}^{J} v_{ij}^{2}, $$

where $r_j$ is the target-normalized residual of the suite-mean spectrum at period $T_j$, $v_{ij}$ is the excursion of the scaled record spectrum outside the compatibility envelope, $\lambda_B$ (shape.bandWeight) weights envelope containment, and $w_j$ applies shape.deficitWeight when the suite mean falls below the target. The containment term conditions every coefficient individually, so no artificial coefficient bounds are needed; an admissible range, when configured, acts only as a backstop.

Prerequisites

A processed record database with a per-process index, e.g. built with runProcess(); the match config points at it through path.database and processID.
A target table (path.target, CSV or Rds) with columns ID, TR, Vs30, p, Tn, and the spectral value named by target.value. It must provide the mean spectrum and the envelope quantiles (for example p05 and p95) for the configured target.ID and target.envelope.ID.

The site JSON

One JSON file per site drives both runners:

{
  "id": "gac30",
  "path": {
    "database": "~/kashimaDB/gmdb.v2",
    "selection": "gmsp/selection/gac30/selection.csv",
    "target": "oq/data/UHSTable.Rds",
    "out": "gmsp/match/gac30"
  },
  "processID": "F25",
  "stage0": { "N": 35 },
  "target": {
    "ID": "GAC30", "TR": 10000, "Vs30": 1500, "p": "mean",
    "siteID": "gac30", "value": "SaF", "units": "g",
    "envelope": { "ID": "GAC30.max", "lower": "p05", "upper": "p95" }
  },
  "ocid": "H1",
  "xi": 0.05,
  "Tn": { "min": 0.04, "max": 5 },
  "shape": {
    "method": "a3",
    "bandWeight": 2,
    "deficitWeight": 4,
    "modalFactor": { "min": 0.001 }
  },
  "imf": { "includeResidue": false },
  "units": { "source": "mm", "target": "mm", "psa": "mm/s2" },
  "plot": { "audit": true, "standard": true },
  "parallel": true,
  "workers": 14,
  "override": true
}

path.selection is written by runStage0() and read by runMatch(): the selection is a product of the workflow, not a hand-maintained list.

Step 1: spectral preselection

runStage0(file = "gmsp/match.json", root = "/path/to/project")

Outputs, written next to path.selection:

selection.csv — the top-stage0.N records by band compatibility;
stage0.spectral.csv — the full scored pool (band metrics, Mw, Repi, PGA, AI per record) for auditing the selection;
stage0.psa.csv — the pool spectra cache. The first run sweeps the database (roughly an hour for several thousand records, checkpointed in batches); any later run — a changed screen, a different N — reuses the cache and completes in seconds.

Sites of the same project share the target period grid, so the cache can be copied between site selection folders to skip the sweep entirely.

Seismological screening

The candidate block restricts the pool before ranking:

"candidate": {
  "ownerID": ["NGAW", "CESMD"],
  "pool": "gmsp/selection/broad/selection.csv",
  "filter": {
    "EventMagnitude": { "min": 6, "max": 8 },
    "Repi": { "max": 200 }
  },
  "exclude": ["95ef1184f1bea4fa"],
  "oneRecordPerEvent": true
}

ownerID — provider whitelist;
pool — a broad runSelect() output used as the candidate pool, for criteria beyond simple ranges (for example unions of magnitude-distance windows);
filter — numeric min/max ranges over any index column;
exclude — record blacklist;
oneRecordPerEvent — keeps the best-ranked record of each event.

Before fixing N, check that enough genuinely compatible records exist:

library(data.table)
DT <- fread("gmsp/selection/gac30/stage0.spectral.csv")
DT[, .(le05 = sum(bandRMSE <= 0.05, na.rm = TRUE),
       le10 = sum(bandRMSE <= 0.10, na.rm = TRUE),
       le20 = sum(bandRMSE <= 0.20, na.rm = TRUE))]

If the compatible tail is thinner than the desired suite size, reduce N rather than padding the suite with incompatible records.

Step 2: match

runMatch(file = "gmsp/match.json", root = "/path/to/project")

Products under path.out:

data/TSW.csv, PSW.csv, IMW.csv     matched histories, spectra, intensities
data/ITS.csv, IPSA.csv             per-mode traceability products
metadata/                          selection, indexes, the JSON used
match/target.csv                   the target band actually used
match/scale.csv                    c_i per record
match/factors.csv                  modal shaping and energy factors
match/stage.csv                    record-scaling summary diagnostics
match/shaping.csv                  modal shaping QA (see below)
match/metrics.csv, eval.csv        per-record optimization diagnostics
match/mean.csv                     suite-mean spectrum
match/recordCompatibility.csv      per-record band diagnostics
match/compatibility.csv            suite-level band summary

Step 3: acceptance diagnostics

match/stage.csv summarizes the suite solve:

RMSE — target-normalized residual of the suite mean;
pgaRatio — suite-mean PGA over target PGA;
insideFraction, meanViolation, maxViolation — envelope containment of the scaled record spectra (maxViolation should stay below about 1);
scaleMin/Median/Max, nAtMin, nAtMax — the coefficient distribution. With band conditioning the coefficients cluster near their natural scale and nAtMin/nAtMax are zero unless a backstop range was configured and binds — which warrants investigation, not acceptance.

match/shaping.csv guards the modal shaping: per record, the minimum and median ratio of the shaped to the source spectrum across the target grid. Ratios far below one at some periods mean the modal shaping removed real spectral content there (for example by excluding the lowest-frequency mode); the selection compatibility paid for by the preselection must survive the shaping.

The in-memory API

The algorithm core of the match runner is exported for direct use on R objects — no files or JSON involved. Three functions cover the two fitting problems and their acceptance diagnostics:

library(gmsp)
library(data.table)

# One record's IMF rows (from TS2IMF), its source TSW (from AT2TS),
# and the target band: Tn (including 0 = PGA), Sa.low, Sa.mean, Sa.high
# in the PSA units of the signals.
Band <- data.table(Tn = c(0, 0.05, 0.10, 0.20),
                   Sa.low = 0.6 * c(100, 50, 35, 20),
                   Sa.mean = c(100, 50, 35, 20),
                   Sa.high = 1.8 * c(100, 50, 35, 20))

# Modal shaping of one record: fits the shape.modalFactor coefficients.
Fit <- fitModalFactor(.x = IMF, target = Band, source = TSW,
                      key = Selection[1L], factorMin = 0.001)
Fit$coeff     # per-mode factors: b, energyFactor, stageFactor
Fit$metrics   # one-row fit diagnostics (RMSE, band containment, ...)
Fit$PSA       # source / start / final spectra, long format

# Record scaling of the suite: fits one shape.recordFactor per record.
PSA <- rbindlist(lapply(Fits, `[[`, "PSA"))  # all fitted records
Suite <- fitRecordFactor(.x = PSA, target = Band,
                         records = Selection$RecordID,
                         factorMin = 0.5, factorMax = 1.6)
Suite$scale   # RecordID, scaleFactor (c_i)
Suite$stage   # one-row solver diagnostics

# Compatibility diagnostics of the scaled suite against the band.
QA <- matchCompatibility(.x = PSW, target = Band, scale = Suite$scale,
                         factorMin = 0.5, factorMax = 1.6)
QA$record     # per-record PGA position, containment, RMSE
QA$summary    # suite aggregate

runMatch() composes exactly these functions over a project record set, adding file orchestration and the product tree below. Historical note: the two fitting problems appear as "Stage A" and "Stage B" in older logs, and match/stage.csv keeps the frozen product label stage = "B".

Composition with the other runners

runProcess()  builds the processed database the pool comes from
runSelect()   optional broad metadata screen feeding candidate.pool
runStage0()   spectral preselection into path.selection
runMatch()    modal shaping + record scaling into path.out
runPlot()     optional raw QA widgets from the products
runExport()   packages products for delivery

Report-grade figures and tables are rendered outside gmsp by the reporting layer, from the materialized CSV products listed above.