Deconvolution Operators#

This page describes the receiver-function deconvolution operators available in MsPASS. Deconvolution estimates an Earth response from an observed component and an estimated source wavelet. In simplified form,

\[d(t) = s(t) * m(t) + n(t),\]

where d is the observed target component, s is the source or reference wavelet, m is the receiver-function estimate, and n is noise. Real source spectra contain notches and real data contain noise, so MsPASS uses regularized inverse filters or sparse iterative methods rather than direct spectral division.

Most users should read this page in four passes:

Skim Essential terms so the method descriptions use familiar words.
Read Wrappers and engines, Common conventions, and Choosing an operator to understand the processing workflow.
Read the method section for the operator you plan to use.
Use Validation and QC workflow to inspect plots and saved metadata.

Essential terms#

target component or target trace: The observed component being deconvolved.
source wavelet or reference wavelet: The estimated source pulse used to build an inverse operator. In many receiver-function workflows this is the vertical component, the P component aligned with the incoming P-wave direction, or an externally prepared stack.
signal_window: The target interval to analyze and usually return.
noise_window: A time interval used to estimate noise for inverse stabilization or iterative stopping.
TimeSeries and PowerSpectrum: MsPASS data objects used by some deconvolution APIs. A TimeSeries stores a sampled time-domain signal. A PowerSpectrum stores a frequency-domain power estimate, commonly used as a reusable noise spectrum.
deconvolution window or receiver-function lag window: The output lag interval configured by parameters such as deconvolution_data_window_start and deconvolution_data_window_end. For the GID wrappers, the signal_window must contain this lag window.
receiver function: The estimated Earth response. Scalar inverse operators return a regularized receiver-function trace directly. GID operators first estimate a sparse impulse response and then shape it into a receiver-function trace.
inverse operator: The regularized operator applied to data or residuals. In GID it is used to form a detection function for choosing candidate spikes. It is not the output shaping wavelet.
generalized iterative deconvolution or GID: A three-component sparse deconvolution family based on Wang and Pavlis (2016). This page uses greedy GID for the classic one-spike-at-a-time iteration and group_sparse for the regularized solver that estimates all lags together.
NS-GID: Noise-stable GID. It is selected with deconvolution_type ns_gid and changes the inverse operator used inside GID.
sparse impulse response: The finite spike series selected by GID. In code this is returned by sparse_output.
shaped receiver function: The sparse impulse response convolved with the output shaping wavelet. This is the finite-bandwidth receiver function returned by getresult for GID.
output shaping wavelet: The wavelet convolved with the sparse impulse response to form the shaped receiver function. In code the preferred method name is output_shaping_wavelet. The older name ideal_output is retained only as a legacy alias.
actual output or resolution kernel: The inverse operator applied to the estimated source wavelet. It describes the resolution of the inverse filter and is returned by actual_output or resolution_kernel.
inverse_wavelet: A diagnostic representation of the inverse filter itself. It can be longer than the returned receiver-function lag window because FFT-based operators work on padded arrays internally.
sparse support: The lag samples retained as nonzero arrivals in the sparse impulse response. A lag group is the three-component coefficient vector at one lag sample.
support threshold: A cutoff used after a sparse solve to decide which lag groups are retained in sparse_output and in the final refit model.
lag-weight penalty: A greedy-GID mechanism that downweights candidate lags after a spike has already been accepted. It changes future candidate selection; it is not an amplitude shrinkage penalty.
group-sparse regularization: A separate GID mode selected with deconvolution_type group_sparse. It estimates the full sparse impulse response with a grouped sparsity penalty instead of using the greedy GID spike picker.
QC metadata: Quality-control metadata stored by wrappers and engines. These scalar fields record processing status, convergence, residual norms, penalty settings, support thresholds, and inverse-stability diagnostics.
CNR: Colored-noise-ratio deconvolution.
SNR: Signal-to-noise ratio.

Wrappers and engines#

MsPASS separates data handling from numerical deconvolution.

RFdeconProcessor and RFdecon: General receiver-function wrappers. They operate on MsPASS data objects, read parameter files, extract wavelet, data, and noise windows, call the selected engine, and save QC metadata in a subdocument. Conventional scalar methods are applied component by component. Generalized iterative deconvolution methods operate on the full three-component seismogram because their spike selection is vector-valued. If a required window is absent, the wrapper returns a killed datum and does not attach the QC subdocument.
CNRRFDecon and CNRArrayDecon: Wrappers for colored-noise-ratio deconvolution. A caller must provide both signal_window and noise_window or provide a precomputed PowerSpectrum noise estimate. If an external wavelet is supplied, it is used for all components; otherwise the configured component of the signal window is used as the source wavelet.
TimeDomainGIDRFDecon and FrequencyDomainGIDRFDecon: Direct wrappers around the GID engines. If signal_window is omitted, the full input time range is used to satisfy loading, while the configured full_data_window remains the returned output interval. If noise_window is omitted, the engine’s parameter-file noise window is used. The analysis interval must contain the configured receiver-function lag window.

GID receiver-function windows#

GID uses four independent time windows. wavelet_window_start and wavelet_window_end isolate the short direct-P source pulse used to build an internal source wavelet. deconvolution_data_window_start and deconvolution_data_window_end define the receiver-function analysis and sparse-search lags; this window must include any late P-to-S conversions of interest, such as P410s and P660s. full_data_window_* defines the returned window and must contain the analysis window. noise_window_* is an independent pre-event interval and can be longer than either signal window.

For example, a mantle-transition-zone RF commonly uses dt=0.05, a source wavelet window of (-5, 20), and analysis/output windows of (-10, 160) with pre-event noise (-200, -10). Configure each selected leaf inverse operator with the same analysis window, target_sample_interval, and shaping_wavelet_dt as the enclosing GID branch. Older parameter files that omit the wavelet-window keys retain their historical behavior: the analysis window is also used for source-wavelet extraction.

Group-sparse support is deliberately narrower than the requested analysis interval when a resolution-kernel column would cross either record edge. The solver requires the complete normalized resolution kernel for every retained lag: partial columns have variable energy and poor refit conditioning, and are therefore rejected rather than interpreted as zero-valued observations. QC records the resulting group_sparse_valid_lag_* interval and the centered_full_fir_support boundary policy. Support is measured on both sides of the resolution kernel’s physical zero-lag sample, so it rejects a tail-only column without discarding a legal near-edge arrival whose full centered kernel is present. Include guard data beyond a late arrival and crop to the desired output window when an arrival near an edge must be resolved. Shortening the analysis interval cannot recover an unsupported edge lag; declare that arrival unsupported when guard data are unavailable.

Group-sparse amplitudes are first obtained from the regularized group-lasso solve and are then optionally de-biased on the selected support. Selected lags can be nearly collinear, so the shared refit measures the Gram condition number with a LAPACK symmetric eigensolve (cubic in support size). Above 1e5 it raises the relative, diagonal-scaled ridge to at least 1e-2; this is a numerical condition guard, not a change to the group-lasso penalty or its default threshold. QC records group_sparse_refit_gram_condition_number, group_sparse_refit_relative_ridge_beta, pre/post residual L2, pre/post maximum amplitudes, and whether the guard or a pre-debias fallback was used. The fallback is used only for a nonfinite or residual-increasing refit.

The lower-level C++ engines do the numerical work. They are useful for tests, diagnostics, and specialized processing tools, but they expect the caller to load the correct data, wavelet, and noise estimates.

Common conventions#

MsPASS treats receiver-function deconvolution as a linear convolution problem. FFT-based operators use zero padding internally, but the returned receiver function is the requested linear lag window, not a wrapped circular-convolution result. For a requested scalar FFT output window of N samples, the working FFT length is at least 2*N - 1 samples, rounded up to the next power of two.

The receiver-function output window is normally configured with deconvolution_data_window_start and deconvolution_data_window_end. Zero lag is the sample whose time is 0.0 relative to the source-wavelet reference time. Diagnostic objects can describe the full padded inverse-operator response: actual_output or resolution_kernel is the inverse operator applied to the source wavelet, while inverse_wavelet is the inverse filter itself. The returned receiver function is the cropped seismic result.

Objective-specific receiver-function presets#

RFdeconProcessor and RFdecon select mtz_plane_wave_lowfreq when called without a parameter file. This is an objective-specific reference preset, not a claim that one parameter set is scientifically optimal for every acquisition or imaging method. It uses a (-10, 120) s analysis/output window, a (-5, 20) s source-wavelet window, and a 0.125 Hz zero-phase Ricker output wavelet. That Ricker bandwidth is a target-specific reference for low-frequency MTZ plane-wave/Kirchhoff imaging, not a generally optimal choice for crustal, high-resolution, or other RF workflows. Its primary standalone reference operators are WaterLevelDecon (water level 0.1) and LeastSquareDecon (damping 1.0).

The installed scalar presets are selected with the preset keyword:

Preset	Intended use	Output shaping
`mtz_plane_wave_lowfreq`	Reference product for MTZ plane-wave/Kirchhoff imaging (selected by default)	Ricker 0.125 Hz
`mtz_plane_wave_highres`	MTZ resolution sensitivity product	Ricker 0.25 Hz
`crustal_plane_wave`	Short-window crustal reference product	Ricker 0.5 Hz
`raw_decon_diagnostic`	Inverse/operator diagnostics only, not direct migration input	none

For example:

from mspasspy.algorithms.RFdeconProcessor import RFdeconProcessor

default = RFdeconProcessor(alg="WaterLevel")
highres = RFdeconProcessor(
    alg="LeastSquares", preset="mtz_plane_wave_highres"
)
historical = RFdeconProcessor(alg="LeastSquares", pf="RFdeconProcessor.pf")

Each preset writes rf_preset into the RF QC metadata. The historical RFdeconProcessor.pf remains selectable for reproducibility, but is a legacy scalar test/example configuration rather than the public scientific default. No-argument GID construction uses the matching long-window/Ricker profile and records gid_workflow=experimental_sensitivity. GID products remain useful sensitivity comparisons; they are not the primary MTZ default.

All image comparisons should hold the output Ricker wavelet fixed. The negative side lobes of a symmetric Ricker wavelet are part of the point-spread function for plane-wave and likely Kirchhoff-style migration; they should not be replaced by an unmatched Butterworth merely because an individual RF looks smoother.

Windows and noise#

Three window concepts appear throughout the deconvolution APIs:

signal_window: The target time interval to analyze and normally the output interval to return.
wavelet or wavelet window: The source estimate used to build the inverse operator. In receiver function workflows this is commonly the vertical or P component, or an external stacked source estimate.

For scalar RFdeconProcessor methods, wavelet_window_start and wavelet_window_end can be distinct from the analysis window. A TimeSeries source wavelet is embedded at its actual time origin on the analysis grid before it is passed to the vector-only scalar engine. Thus a short (-5,20) P pulse paired with an analysis window beginning at -10 s does not acquire a spurious five-second lag shift. A raw numeric vector has no intrinsic time base: pass wavelet_t0 to RFdecon or RFdeconProcessor.loadwavelet to locate its first sample on the physical time axis, or use a TimeSeries. Omitting wavelet_t0 retains the legacy interpretation that vector sample zero is the analysis-window start and emits a warning. The same rule applies to bare external_wavelet vectors passed to the time- and frequency-domain GID wrappers. A TimeSeries wavelet always keeps its own t0 and sample interval.

The GID receiver-function wrappers and engines interpret their configured windows, sample shifts, and output times as P-relative lag coordinates. They therefore reject a UTC analysis datum before loading an external source into a reused engine. Convert event data with ator(P-arrival epoch) before GID RF processing. This is a contract of the GID RF workflow, not a general restriction on TimeSeries or on the datum-independent low-level loadwavelet method.

noise_window: A time interval or spectrum used to stabilize the inverse or to decide when an iterative method should stop.

The word “noise” has several meanings:

source/wavelet noise: Noise used to stabilize the inverse operator denominator. Multitaper, CNR, and NS-GID use this information to avoid excessive inverse gain in weak or noisy source-wavelet bands.
target/data noise: Noise already present in the component being deconvolved. Stable inverse operators reduce amplification of this noise indirectly by regularizing the source-wavelet inverse.
residual-domain noise: Noise used by iterative methods to decide whether another sparse spike is significant and when iteration should stop.

For ns_gid and group_sparse, an external TimeSeries noise estimate can be used for inverse-operator regularization and, when loaded as residual noise, residual-domain stopping or support decisions. An external PowerSpectrum can regularize only the NS-GID inverse operator; a residual-domain noise window or time-domain noise estimate is still needed for sparse iteration stopping or group-sparse support decisions.

Choosing an operator#

No deconvolution operator is best for every receiver-function problem. Start with the failure mode you need to control:

Goal or data condition	Recommended starting point
Simple scalar receiver functions with moderate noise	`LeastSquareDecon` or `WaterLevelDecon` through `RFdecon`
Source spectrum is noisy or has strong notches	`NoiseStableDecon` for scalar tests, or `ns_gid` for GID
Need colored-noise regularization on three components	`CNRDeconEngine` through `CNRRFDecon` or `CNRArrayDecon`
Need sparse three-component arrivals and per-spike diagnostics	`TimeDomainGIDDecon` or `FrequencyDomainGIDDecon`
Noisy GID repeatedly picks the same strong arrival	`deconvolution_type ns_gid` with `adaptive_memory` lag penalty
Want shared sparse support across components	`deconvolution_type group_sparse`

The distributed GID parameter files use deconvolution_type ns_gid with the adaptive_memory lag penalty. Validation sweeps selected that combination because it is more resistant to noise-driven false picks than the legacy least_square inverse. Use group_sparse when a shared-support sparse prior matches the problem, and reserve least_square for explicit legacy comparisons or diagnostics.

Switching GID methods safely#

The shipped GID defaults changed from deconvolution_type least_square to deconvolution_type ns_gid with lag_weight_penalty_function adaptive_memory. New code should switch GID modes through the public Python API instead of editing the installed parameter files. The examples below assume normal parameter-file lookup, meaning PFPATH includes the shipped MsPASS pf directory.

For one-off receiver-function calls, pass the GID options directly to RFdecon:

from mspasspy.algorithms.RFdeconProcessor import RFdecon

# Current shipped default: noise-stable GID plus adaptive-memory lag weights.
rf_default = RFdecon(seis, alg="TimeDomainGID")

# Explicit legacy comparison.  Keep this as a diagnostic, not the default.
rf_legacy = RFdecon(
    seis,
    alg="TimeDomainGID",
    deconvolution_type="least_square",
)

# Group-sparse GID: one regularized solve with shared lag support.
rf_group_sparse = RFdecon(
    seis,
    alg="TimeDomainGID",
    deconvolution_type="group_sparse",
)

# Make the greedy lag-weight penalty stronger or weaker.
rf_weaker_penalty = RFdecon(
    seis,
    alg="TimeDomainGID",
    deconvolution_type="ns_gid",
    lag_weight_penalty_function="adaptive_memory",
    gid_parameters={"lag_weight_penalty_scale_factor": 0.2},
)

For repeated processing, build a configured processor once and reuse it:

from mspasspy.algorithms.RFdeconProcessor import RFdecon, RFdeconProcessor

processor = RFdeconProcessor(
    alg="TimeDomainGID",
    deconvolution_type="group_sparse",
    gid_parameters={"group_sparse_active_threshold": 0.03},
)

rf = RFdecon(seis, alg="TimeDomainGID", engine=processor)

For lower-level workflows that use the GID engines directly, use make_gid_engine or the corresponding parameter helpers:

from mspasspy.algorithms.RFdeconProcessor import (
    make_gid_engine,
    make_gid_pf,
    make_gid_pf_text,
)

engine = make_gid_engine(
    alg="FrequencyDomainGID",
    gid_mode="ns_gid",
    gid_penalty_function="adaptive_memory",
)

pf = make_gid_pf(
    alg="TimeDomainGID",
    deconvolution_type="group_sparse",
    gid_parameters={"group_sparse_active_threshold": 0.03},
)
pf_text = make_gid_pf_text(
    alg="TimeDomainGID",
    deconvolution_type="least_square",
)

make_gid_engine returns the lower-level C++ GID engine object. Pass that object to TimeDomainGIDRFDecon or FrequencyDomainGIDRFDecon; pass an RFdeconProcessor object to the high-level RFdecon(..., engine=...) interface.

These settings are separate layers. The deconvolution_type method layer and its alias gid_mode choose the GID inverse or solver mode, such as ns_gid, least_square, or group_sparse. lag_weight_penalty_function and its alias gid_penalty_function change only the greedy GID lag-selection penalty; they do not enable grouped sparsity. The group-sparse support threshold keys, such as group_sparse_active_threshold, group_sparse_active_threshold_scale, and group_sparse_active_threshold_quantile, are applied after the group-sparse solve to decide which lag groups remain in the sparse output.

The modes above are the recommended user-facing choices. Legacy and advanced inverse modes are still accepted where supported, including water_level, multi_taper, cnr/cnr3c, the noise_stable alias for ns_gid, and the group_lasso/sparse_group_lasso aliases for group_sparse.

Scalar inverse operators#

Scalar operators load one source wavelet and one target trace. They return a regularized inverse-filter receiver function, not a sparse impulse response. The objective-specific scalar presets apply the same Ricker output shaping wavelet to all scalar branches. Older parameter files can still select shaping_wavelet_type none for diagnostic/raw output, but raw and shaped results should not be compared as equivalent migration inputs.

LeastSquareDecon

Frequency-domain damped least-squares deconvolution. The inverse operator has the form

\[S_g^{-1}(\omega) = \frac{\overline{S(\omega)}}{|S(\omega)|^2 + \mu},\]

where mu is the damping term. Smaller damping improves resolution for clean data but amplifies noise and spectral notches. Larger damping produces a smoother, more stable estimate.

WaterLevelDecon

Frequency-domain deconvolution with a water-level floor on the source spectrum. Raising the water level improves stability but reduces resolution.

MultiTaperPowerXcor and MultiTaperPowerSpecDiv

Multitaper frequency-domain operators that use discrete prolate spheroidal sequence (DPSS) tapers to stabilize source-power estimates. The final receiver-function estimate is formed with the untapered data window so delayed converted phases are not weighted by the taper value at their arrival time. In receiver-function workflows, the noise vector should normally describe the source/wavelet component, not the radial or transverse target component.

These MsPASS operators are not paper-faithful Park and Levin (2000) multitaper correlation estimators. They use DPSS spectra to stabilize the source-power denominator, but apply the final inverse operator to the untapered wavelet phase and untapered data spectrum.

NoiseStableDecon

Noise-aware scalar inverse used by ns_gid and exposed as a standalone validation operator. It applies gain limits and frequency-dependent regularization from the source/noise estimate:

\[G(f)=B(f)\frac{\overline{S(f)}}{|S(f)|^2+\mu(f)}.\]

The operator does not pick sparse spikes and does not apply the GID output shaping wavelet.

TimeDomainLeastSquareDecon

Time-domain least-squares deconvolution. This operator builds a linear convolution system for the requested lag window and solves the regularized normal equations. It does not build a circular convolution matrix.

Three-component CNR operators#

CNRDeconEngine and CNR3CDecon implement colored-noise-ratio receiver-function deconvolution. They estimate or load a noise spectrum and regularize the source spectrum as a function of colored noise level and SNR. colored_noise_damping adds a frequency-dependent damping term to the normal-equation denominator. generalized_water_level raises low-SNR source amplitudes before division.

CNRDeconEngine is the current engine used by the Python wrappers and by the CNR inverse mode inside GID. CNR3CDecon is the older three-component prototype kept for compatibility.

The legacy CNRDeconEngine.process(d, noise_spectrum, fl, fh) entry point continues to derive an output shaper from fl and fh. GID instead uses the configured-shaping path, so its inverse-band choices no longer overwrite a configured Ricker frequency. This keeps CNR output shaping independent from the colored-noise inverse regularization band and makes its 0.125/0.25 Hz products comparable with the other deconvolution methods.

Generalized iterative deconvolution (GID)#

TimeDomainGIDDecon and FrequencyDomainGIDDecon estimate a sparse three-component impulse response. Accepted spikes are convolved with the configured output shaping wavelet to produce the shaped receiver function returned by getresult. The raw sparse impulse response is available through sparse_output for diagnostics and validation.

Greedy GID workflow#

Greedy GID is the classic one-spike-at-a-time sparse method. At each iteration the engine applies a configured inverse operator to the current residual to form a detection function,

\[a(t) = g(t) * r(t),\]

where r is the current residual and g is the configured inverse operator for the selected GID mode. The largest acceptable three-component detection-function peak becomes a candidate spike. The engine subtracts that candidate’s predicted contribution from the data-domain residual and keeps the spike only if the residual decreases.

For legacy greedy modes, residual_fractional_improvement_floor implements the Wang & Pavlis (2016) Eq. (15) post-acceptance state test. The highest-ranked candidate is accepted whenever it decreases the residual; the updated residual state then stops the iteration when its fractional L2 improvement reaches the configured floor. It is not a candidate-level rejection or a lower-ranked-candidate scan. ns_gid instead applies its noise-significance and fractional-improvement gate before committing the highest-ranked candidate. A residual-decreasing NS-GID candidate that is at or below that floor terminates the epoch without accepting a lower-ranked candidate.

The inverse operator is therefore used to choose candidate spike locations and amplitudes. It is not the final receiver-function representation. The reported GID receiver function is the sparse impulse response convolved with the output shaping wavelet.

Time-domain and frequency-domain engines#

The two GID engines implement the same sparse model with different internal machinery. The time-domain engine performs the iterative residual update in the time domain. The frequency-domain engine builds inverse and resolution quantities with frequency-domain operations. Both engines update residuals with a compact, peak-normalized residual-update kernel derived from the inverse operator’s actual output. The public actual_output() method returns the full inverse-operator resolution diagnostic, not the compact residual-update kernel.

Because scalar leaves intentionally normalize their public resolution diagnostic, GID measures the corresponding unnormalized source-through-inverse zero-lag gain and applies its reciprocal to the inverse-domain data and noise quantities before sparse selection. This keeps the data, residual, candidate threshold, and zero-lag-normalized update kernel in one amplitude domain. It is especially important for heavily damped least-squares leaves; public actual_output() behavior is unchanged.

The multi_taper inverse mode is still a frequency-domain inverse operator, even when selected inside TimeDomainGIDDecon. In that case “time domain” describes the sparse iteration and residual update, not the spectral method used by the inverse operator.

`NS-GID` inverse mode#

NS-GID is selected with deconvolution_type ns_gid in either GID engine. It is a noise-stable inverse mode used inside the GID iteration, not a separate preprocessing workflow. The inverse operator is

\[G(f) = B(f)\frac{\overline{S(f)}}{|S(f)|^2 + \mu(f)},\]

where S is the source wavelet spectrum, B is the optional band-limiting or reliability taper, and mu(f) is a frequency-dependent damping term. The damping combines a relative spectral floor, noise-spectrum information, and an explicit maximum-gain constraint controlled by ns_gid_gain_max.

Use ns_gid when noisy or spectrally weak source wavelets make greedy GID with a less noise-stable inverse pick noise-generated spikes. External wavelets can be supplied through the GID loadwavelet APIs or wrapper external_wavelet argument. External TimeSeries noise can support both inverse regularization and residual-domain stopping when loaded as residual noise. External PowerSpectrum noise only regularizes the inverse operator; a residual-domain noise window is still needed for spike significance and stopping. In group_sparse mode, that inverse regularization can also change the automatically selected group_sparse_lambda because the automatic lambda uses the NS-GID inverse-filtered noise threshold when it is available.

Convergence and stopping#

GID stops when the configured iteration limit is reached or when a stopping test says another spike is not useful. Common stopping tests include small fractional residual improvement, residual energy reaching a noise floor, no acceptable candidate lag, a lag-weight penalty exhausting valid candidates, and, when ns_gid is selected, optional ns_gid_max_spikes and noise-significance thresholds.

Time-domain and frequency-domain GID use the same iteration-cap behavior. Reaching maximum_iterations stops the iteration, returns the best accepted sparse model after final amplitude refit, and records gid_stop_reason="max_iterations" with gid_converged=false. Stops caused by residual or lag-weight floors are reported as converged because the engine found a configured reason not to add more spikes.

QC metadata#

QCMetrics records both processing status and diagnostic context. Useful first-pass GID fields include decon_operator, decon_processed, gid_converged, gid_stop_reason, gid_iterations, gid_number_spikes, residual_rms_final_fraction, and residual_rms_reduction_fraction. residual_L2_initial and residual_L2_final are retained for compatibility, but are unnormalised Euclidean norms and must not be compared across records or analysis-window lengths. residual_rms_initial and residual_rms_final remove only the analysis-window-length factor; use the fractional fields for comparisons across records with different amplitudes. residual_energy_final_fraction and residual_energy_reduction_fraction are the corresponding squared-norm fractions, useful when the QC policy is expressed in explained residual energy. residual_3c_rms_initial and residual_3c_rms_final are component RMS values computed with denominator sqrt(3*npts) and remove only record length and component-count dependence; they remain amplitude-, bandwidth-, and leaf-dependent and are not cross-record quality measures by themselves. Use the dimensionless residual fractions or the noise-normalized NS-GID fields for cross-record QC. The explicit residual_L2_metric label identifies these as legacy unnormalized Euclidean norms in the new gain-normalized inverse domain; it is not a raw-input-amplitude metric. GID also records gid_leaf_raw_zero_lag_gain, gid_inverse_domain_amplitude_scale, and gid_inverse_domain_scaling_policy. Legacy Eq. (15) audits are available as gid_legacy_eq15_candidates_tested, gid_legacy_eq15_candidates_rejected, and gid_legacy_eq15_rejected_lag_seconds. These legacy names are compatibility pre-trial scan counters, not Eq. (15) state-test counts; the identical, explicit gid_legacy_eq15_pretrial_scan_* fields make that scope clear. Use gid_legacy_eq15_post_acceptance_state_tests and gid_legacy_eq15_post_acceptance_floor_stops for post-acceptance Eq. (15) QC. Legacy leaf modes preserve their historical greedy rule: the highest-ranked candidate is accepted whenever it reduces the residual, then Eq. (15) tests the resulting state improvement. It is not a candidate-level gate and the engine never scans lower-ranked candidates to bypass a sub-floor best candidate. The legacy audit fields are therefore diagnostic only; a post-acceptance floor stop is recorded as gid_stop_reason="fractional_improvement_floor" and gid_legacy_eq15_stop_detail="post_acceptance_fractional_improvement_floor". For ns_gid the same fractional floor is evaluated before committing its highest-ranked noise-significant candidate; a decreasing but sub-floor trial is recorded as fractional_improvement_floor_rejected and terminates that epoch without accepting a lower-ranked candidate.

After a terminal NS-GID refit, ns_gid_final_scan_best_trial_* reports the best trial across the complete read-only candidate audit, while ns_gid_final_scan_acceptable_candidate_remaining and ns_gid_final_scan_global_acceptable_candidate_count report whether any audited candidate clears the floor. These are deliberately distinct from the greedy decision: ns_gid_final_scan_decision_candidate_* and ns_gid_final_scan_decision record the first decreasing ordered candidate, which alone controls resume versus the fractional-floor stop. This prevents a later candidate from overriding a sub-floor greedy decision while retaining complete diagnostic visibility.

For compatibility, the frequency-domain legacy Eq. (15) floor stops at improvement <= floor; the time-domain legacy implementation retains its historical strict improvement < floor comparison. New configurations should not depend on an exact floating-point equality at this boundary.

Greedy lag-weight penalty runs also record gid_penalty_function, gid_penalty_scale_factor, gid_penalty_width, gid_penalty_effective_width, lag_weight_Linf_final, and lag_weight_L2_final. Adaptive-memory runs add fields such as gid_adaptive_penalty_enabled, gid_penalty_noise_amplitude, gid_penalty_last_confidence, gid_penalty_last_decay_factor, gid_penalty_memory_Linf_final, and gid_penalty_memory_L2_final.

gid_noise_samples_loaded and gid_noise_samples_used always describe the GID residual/stopping noise record, regardless of the selected inverse mode. Leaf regularization inputs are reported separately as gid_leaf_noise_samples_loaded, gid_leaf_noise_samples_used, and gid_leaf_noise_truncated. gid_residual_external_noise_used indicates that an externally loaded TimeSeries supplied the GID residual/stopping record, while gid_leaf_external_noise_used indicates that it was consumed by the leaf inverse operator. The legacy-compatible gid_external_noise_used is true if either role used external noise.

When deconvolution_type ns_gid is active, additional ns_gid_* fields record inverse stability and stopping diagnostics, including ns_gid_stop_reason, ns_gid_converged, ns_gid_peak_threshold, ns_gid_noise_amplitude_rms, ns_gid_noise_component_sigma_rms, ns_gid_gain_max_actual, and ns_gid_effective_bandwidth_fraction. The Python wrappers store these QC values in receiver-function metadata subdocuments, so they are saved with normal database records.

For NS-GID, ns_gid_peak_threshold is an amplitude threshold on the three-component vector norm at one candidate lag. With the shipped configuration it is the larger of the configured empirical quantile (ns_gid_peak_threshold_empirical) and ns_gid_peak_sigma_threshold times ns_gid_noise_component_sigma_rms (ns_gid_peak_threshold_sigma). The requested robust reference, ns_gid_noise_component_sigma_rms_robust, is the RMS of three median-centered, signed-component MAD sigma estimates. It is an explicitly defined scalar-component noise sigma, not a statistic of nonnegative vector amplitudes. The shipped multiplier of 3.0 consequently means three scalar-component sigma for this robust reference. ns_gid_empirical_peak_threshold and ns_gid_sigma_peak_threshold are literal compatibility aliases for ns_gid_peak_threshold_empirical and ns_gid_peak_threshold_sigma, respectively.

Sparse or quantized nonzero noise can have a zero (or roundoff-level) MAD. When that occurs, the engine uses the ordinary signed-component RMS aggregate ns_gid_noise_component_rms_aggregate instead and records both the used scale in ns_gid_noise_component_sigma_rms and ns_gid_noise_component_sigma_rms_fallback_used=true. If neither the robust scale nor that ordinary scalar-component RMS reference is finite and positive, processing fails rather than emitting a zero sigma threshold. Thus disabling the empirical guard cannot silently bypass NS-GID filtering.

ns_gid_noise_amplitude_rms remains the ordinary inverse-filtered 3C vector-amplitude RMS, and ns_gid_noise_amplitude_robust remains the quadrature combination of the same three MAD sigma estimates. They are retained as vector-noise diagnostics but do not define the sigma multiplier. ns_gid_component_noise_rms_0 through _2 continue to expose ordinary component RMS contributions. The empirical threshold is still a quantile of the 3C vector-amplitude distribution, and the max policy between it and the sigma threshold is unchanged.

group_sparse uses the same inverse-filtered noise-threshold computation for its automatic group_sparse_lambda (when that parameter is zero), so the scalar-component MAD policy also defines its automatic lambda reference. Its QC records group_sparse_noise_threshold, group_sparse_peak_threshold_empirical, group_sparse_peak_threshold_sigma, group_sparse_use_empirical_noise_threshold, and group_sparse_peak_threshold_controlling_term. The latter has the same meaning as the NS-GID controlling term: empirical, sigma, or tie when the empirical guard is enabled, and sigma_empirical_disabled when it is disabled. QC also records the corresponding group_sparse_noise_component_* requested/used-scale fields. Set an explicit positive group_sparse_lambda to make the regularization strength independent of this noise-scale policy.

Historical real-data diagnostic and scope#

An earlier diagnostic screen used 332 real RF records from 10 station folds (31–39 records per fold), using analysis, wavelet, and noise windows of (-10, 160), (-5, 20), and (-200, -10) seconds at dt=0.05. It was run with the then-current implementation on Vista ARM64; it is not an execution of the corrected current engine. Each parameter sweep reused these same 332 records; the resulting TD/FD sweep rows are repeated measurements, not thousands of independent observations.

This historical screen must not be used to validate or retain shipped defaults: it included the now-corrected frequency-domain threshold-domain error and showed low iteration/high residual-energy symptoms. It is retained only as a diagnostic record. Default selection remains provisional pending a common-event EarthScope benchmark, comparison with other MsPASS deconvolution methods, and waveform/stack QC.

Table 11 Active NS-GID controls and historical diagnostic notes#
Control	Shipped value	Historical note; not a default-validation decision
`maximum_iterations`	100	200/400 reduced cap exits, but increased negative-time energy, IQR, and incoherent content. Historical observation only; no current default decision follows.
`residual_fractional_improvement_floor`	`1e-4`	This is the 0.01% criterion of Wang and Pavlis equation 15. `1e-5` produced no waveform-quality gain in that screen; no current default decision follows.
`ns_gid_refit_interval`	5	Intervals 1, 5, and 20 were stable in the screen. Historical observation only.
`ns_gid_residual_noise_ratio_floor`	1.0	1.2 produced about 21% one-spike solutions, consistent with underfit; 0.8 produced no coherent gain in that screen. Historical observation only.
`ns_gid_peak_sigma_threshold`	3.0	This remains exactly three robust scalar-component MAD sigmas (with the documented fallback), not three vector-amplitude sigmas. This defines the parameter but does not recommend a current value.
empirical threshold	`q=.995`, on	Lowering to .9925 reduced the under-5-spike fraction from 18% to 12% and raised median iterations from 19/20 to 23.5/24.5 (TD/FD), without increasing moveout-corrected P410/P660 coherence; negative-time/IQR spread broadened. .990/.975 and empirical-off showed stronger noise-overfit signatures. This pre-fix observation is not a recommendation for .995.
`ns_gid_max_spikes`	0 (unlimited)	A hard cap of 50 degraded waveform consistency in that screen. Historical observation only.
`ns_gid_ridge_beta`	`1e-10`	`1e-8` degraded waveform consistency and `1e-6` had no stable advantage in that screen. Historical observation only.
FD `residual_ratio_floor`	0.01	This is active in frequency-domain NS-GID, but was not reached before NS criteria in this sample. Historical observation only.
TD lag/residual-probability controls	as shipped	The TD lag-weight and residual-probability controls are not NS-GID stop criteria (though penalty configuration can still affect candidate scoring). They were not retuned as NS stopping defaults.

The final operational peak threshold is the maximum of the scalar-sigma and empirical vector-amplitude thresholds, so it need not equal exactly the stated number of scalar sigmas when the empirical guard dominates. ns_gid_peak_threshold_controlling_term records empirical, sigma, or tie for that max operation. When ns_gid_use_empirical_noise_threshold is false, the operational threshold is exactly the sigma threshold and the controlling term is explicitly sigma_empirical_disabled; that boolean is also published in QC. ns_gid_peak_threshold_scope is always pointwise_candidate_lag. Thus ns_gid_peak_probability_threshold is a single-lag 3C-amplitude quantile, not a global false-alarm probability over a receiver-function search. The three-component norm has a Maxwell, rather than scalar-Gaussian, tail, so a multiplier of three component sigmas is not a conventional one-component 3 sigma exceedance rule. The initial ns_gid_initial_stationary_null_expected_noise_exceedances QC field makes the resulting multiple-lag exposure explicit; it is a diagnostic estimate, not a calibrated family-wise-error guarantee.

ns_gid_residual_noise_rms_ratio is the canonical, dimensionless stopping ratio: final ordinary three-component residual-amplitude RMS divided by the ordinary inverse-filtered three-component noise-amplitude RMS. It is the ratio used for the NS-GID residual-noise stop and is independent of noise-window length. ns_gid_residual_rms_ratio is a compatibility alias with the identical value. ns_gid_residual_rms aliases ns_gid_residual_rms_final (RMS amplitude units), while ns_gid_noise_rms aliases ns_gid_noise_amplitude_rms (the same units). ns_gid_residual_noise_ratio and ns_gid_residual_l2_ratio_legacy retain the earlier raw-L2 ratio for backward-compatible QC interpretation, but no longer drive stopping.

For NS-GID, residual_fractional_improvement_floor is evaluated from the raw-L2 reduction caused by the accepted candidate before a periodic ridge amplitude refit. A ridge refit is retained as the next iteration state, but its residual adjustment cannot by itself cause candidate-fractional convergence. ns_gid_fractional_improvement_final therefore means the final accepted candidate increment; the separately reported ns_gid_fractional_improvement_state_final describes the corresponding post-refit state change. Each ns_gid_iteration_N trace also records residual_l2_before_candidate, residual_l2_trial_pre_refit, residual_l2_post_refit, candidate and state fractional improvements, and periodic_refit_applied/final_refit_applied flags. Candidate rows are historical records: a later terminal refit never rewrites their local stop condition. Final-refit audit fields are published separately in the global ns_gid_final_scan_* metadata. trial_evaluated and metric_available distinguish a residual trial from an early significance or spike-cap exit. Every row reports the measured residual_l2_before_candidate. Trial, post-refit, and improvement fields are present only when a residual trial was actually evaluated; a rejected trial reports its restored pre-candidate L2 state and zero state improvement.

GID constructs a union preprocessing grid spanning the full source and analysis intervals, then maps the leaf response back onto the unchanged analysis grid; it does not relabel the output after vector-only leaf deconvolution. The QC keys gid_analysis_t0, gid_wavelet_t0, gid_analysis_samples, gid_wavelet_samples, and gid_wavelet_alignment_offset_samples record that mapping. For NS-GID, ns_gid_iteration_N_candidate_lag_samples, _candidate_lag_time, _candidate_amplitude, _threshold, _significance, and _accepted make each candidate decision inspectable. The matching _post_residual_rms_ratio is the residual/noise RMS ratio after a candidate is accepted, and _stop_condition is initially continue or the local rejection reason. After final ridge refitting, the final trace row is updated with the final residual ratio for QC, but its recorded stop condition remains the decision made at that iteration. A dirty support is amplitude-refit once; the refit then scans only off-support positive-weight lags. A raw-significant lag resumes NS-GID only if its trial residual is strictly decreasing and clears residual_fractional_improvement_floor. Hard iteration/spike caps still terminate. ns_gid_refit_epochs and ns_gid_refit_resume_count record these transitions.

The repository includes deterministic synthetic checks of the shipped 3.0/ 0.995 profile. The historical screen above does not provide workflow evidence for retaining it, and no current real-data result establishes a calibration. Fixed-lag stacks and IASP91 ray-parameter/moveout checks in that screen had coherence near a random 1/N baseline. They do not constitute P410s/P660s detection or depth inference, and the waveform set is not a representative mantle-conversion benchmark.

Warning

ScalarDecon public virtual diagnostics now return TimeSeries rather than CoreTimeSeries. Downstream custom overrides/callers and prebuilt extensions must update to that signature and be rebuilt. This is an intentional API migration so diagnostic waveforms retain metadata.

For bounded synthetic controls only, a conservative profile is ns_gid_peak_sigma_threshold 4.0 and ns_gid_peak_probability_threshold 0.999 (with the empirical threshold enabled). Its 4.0 multiplier is applied to the robust scalar-component sigma defined above, while the 0.999 empirical vector-amplitude quantile remains the guard for non-Gaussian noise. This profile retained false-positive control in the bounded short-window synthetic regression after the scalar-component-sigma change (signal (-8, 22), noise (-35, -8), with synthetic arrivals at 0, 3, 8, and 18 seconds); it is not a shipped default, a real-data recommendation, or evidence for overriding any pending EarthScope/cross-method calibration.

In a Figure-8-style downstream QC pass, inspect gid_number_spikes, the stop reason, and event-stack consistency together. Very-low-spike solutions are candidates for editing or exclusion after waveform review; the engine does not auto-kill them.

For a beginner-oriented guide to reading these fields after a plot run or database save, see Validation and QC workflow.

Lag-weight penalty framework#

Lag-weight penalties apply only to greedy GID. They modify the score used to choose the next candidate lag after a spike has already been accepted. They do not shrink accepted amplitudes, and they are separate from the grouped-sparsity solve used by deconvolution_type group_sparse (see Group-sparse regularized GID). In all greedy modes the data-domain residual decrease test remains the final acceptance rule.

Recommended use#

The easiest way to choose a penalty is by failure mode:

Situation	Suggested setting
Noisy receiver functions repeatedly pick one strong arrival	`deconvolution_type ns_gid` with `lag_weight_penalty_function adaptive_memory`
Closely spaced true arrivals are the scientific target	`lag_weight_penalty_function none`
Visual diagnostics	compare `none` and the active penalty
Fixed-width controlled tests	`cosine_taper` or `boxcar`
Resolution-footprint studies	`resolution_kernel` or `shaping_wavelet`

The shipped GID parameter files set deconvolution_type ns_gid, lag_weight_penalty_function adaptive_memory, and lag_weight_penalty_scale_factor=0.35. Validation sweeps across moderate and large noise levels selected this as the default because the legacy least_square inverse can generate noise-driven false picks even with the same lag penalty. Treat it as a robust starting point, not a universal optimum.

How the greedy penalty enters GID#

At iteration \(k\), let \(\mathbf{a}^{(k)}_j\) be the three-component detection-function vector at lag sample \(j\), and let \(\ell^{(k)}_j \in [0, 1]\) be the accumulated lag weight. For NS-GID, the raw-amplitude threshold \(T\) is a per-valid-lag test, not a full-search false-alarm probability. The picker first restricts its eligible set, then applies the lag penalty:

\[\mathcal{E}^{(k)} = \{j:\ell^{(k)}_j > 0, \left\|\mathbf{a}^{(k)}_j\right\|_2 \ge T\}, \qquad j_* = \operatorname*{arg\,max}_{j\in\mathcal{E}^{(k)}} \ell^{(k)}_j\left\|\mathbf{a}^{(k)}_j\right\|_2 .\]

If \(\mathcal{E}^{(k)}\) is empty, NS-GID reports candidate_not_significant (subject to the final-refit audit below). Other greedy GID modes retain their legacy unconstrained weighted-energy selection.

After accepting a spike at \(j_*\), the engine multiplies nearby lag weights by a compact penalty kernel \(p_m\):

\[\ell^{(k+1)}_j = \operatorname{clip}\left(\ell^{(k)}_j p_{j-j_*}, 0, 1\right).\]

Outside the penalty support, \(p_m=1\). Repeated hits near the same arrival multiply down the same neighborhood, making the penalty a soft anti-cycling cost. The penalty changes which lag is tried next; the data-domain residual decrease test still decides whether the candidate is accepted.

Current lag-weight penalty methods#

none: Leaves the greedy selector unchanged. It is the safest support-recovery choice when real arrivals may be tightly clustered.
boxcar: Applies fixed-width hard local suppression. It is easy to interpret but can move false picks to the first unpenalized sample outside the window.
cosine_taper: Applies fixed-width smooth compact suppression. The accepted sample receives the strongest penalty and nearby samples are downweighted gradually.
shaping_wavelet: Derives the applied penalty footprint from the requested output shaping wavelet. This is useful for controlled comparisons where the requested output pulse is the desired resolution reference.
resolution_kernel: Derives the footprint from the compact actual residual-update kernel. This reflects inverse-operator regularization, trimming, and peak normalization.
adaptive_memory: Derives the footprint from the resolution kernel and changes memory strength and retention from event to event. Strong, well-localized arrivals create stronger local memory. Broad or low-confidence arrivals create weaker memory that decays quickly.

Adaptive memory theory#

adaptive_memory implements a finite-memory local penalty. After each accepted spike it asks three questions:

Which lags are in the same ambiguous state as the accepted spike?
How reliable is the accepted spike as an explanation of the residual?
How long should that local penalty persist?

Let \(M^{(k)}_j=-\log \ell^{(k)}_j\) be the local penalty cost. The picker can be written as

\[q^{(k)}_j = \exp\left[-2M^{(k)}_j\right] \left\|\mathbf{a}^{(k)}_j\right\|_2^2 .\]

An accepted spike adds local memory shaped by the resolution-kernel coherence, scaled by a noise-normalized confidence, and retained according to the specificity of the footprint. In the practical extremes, high-SNR sharp arrivals create durable local memory, high-SNR broad arrivals create strong but quickly forgotten memory, low-SNR sharp arrivals create cautious memory, and low-SNR broad arrivals create little durable memory. This discourages immediate cycling without permanently forbidding nearby lags.

The algorithm estimates confidence from a search-adjusted three-component detection statistic. If \(\Sigma_s\) is the local detection-function noise covariance at the accepted sample \(s\), a natural single-lag statistic is

\[Z_s^2 = \mathbf{a}_s^T \Sigma_s^{-1}\mathbf{a}_s .\]

Because GID picks the maximum over all currently valid lags, a pure-noise maximum also grows with search size. With \(N_\mathrm{valid}\) valid lag samples, the helper normalizes by a three-component upper-tail bound,

\[E_\mathrm{search} = 1 + 2\sqrt{\frac{2\log N_\mathrm{valid}}{3}} + \frac{4\log N_\mathrm{valid}}{3},\]

using \(E_\mathrm{search}=1\) when there is only one valid lag. The bounded confidence is then

\[\begin{split}B_s = \begin{cases} 0, & Z_{*,s} \le 1, \\ 1 - Z_{*,s}^{-2}, & Z_{*,s} > 1 , \end{cases} \qquad Z_{*,s}^2 = Z_s^2/E_\mathrm{search}.\end{split}\]

The ambiguity footprint comes from the coherence of shifted resolution-kernel atoms. A broad footprint has low information concentration and should be forgotten faster than a sharp footprint. The implemented specificity is

\[S_s = 1 - \frac{\log N_{\mathrm{eff},s}} {\log N_\mathrm{valid}}, \qquad N_{\mathrm{eff},s} = \frac{\left(\sum_j A_s(j)\right)^2}{\sum_j A_s(j)^2},\]

where \(A_s(j)\) is the applied footprint weight. The retention strength is \(\gamma_s=B_sS_s\). Conceptually, the memory update is

\[M^{(k+1)}_j = R^{(k)}_j M^{(k)}_j - \log\left[1-\alpha B_s A_s(j)\right],\]

where \(\alpha\) is lag_weight_penalty_scale_factor and \(R^{(k)}_j\) is the retention already stored at lag \(j\). This is a temporary local penalty, not a permanent exclusion rule.

Benchmark interpretation#

The validation results should be read as evidence for guarded anti-cycling behavior, not as proof that adaptive_memory is universally best. In the moderate-noise synthetic sweep, adaptive_memory reduced some repeated-pick failures and false positives. In dense close-arrival cases, none can remain better because any local anti-cycling memory can suppress the next true arrival. In high-noise sweeps, search-adjusted confidence can make adaptive_memory collapse toward none once selected peaks are no longer strong relative to the full lag search.

Use the penalty plots and QC fields to decide whether the penalty helped the specific data and noise setting being processed.

Serialization and distributed use#

Direct TimeDomainGIDDecon and FrequencyDomainGIDDecon engine objects are pickleable for distributed wrapper use. Their pickle state preserves the parameter-file configuration, successful changeparameter() updates to the underlying inverse operator, and any externally loaded wavelet, TimeSeries noise, or PowerSpectrum noise. It does not preserve loaded seismograms, processed receiver functions, residuals, sparse spike trains, or runtime QC state.

For Dask or Spark jobs, build and configure the reusable deconvolution processor on the driver, then scatter or close over that configured object. Avoid loading per-datum scalar data into an RFdeconProcessor before scattering it, because scalar compatibility mode intentionally preserves cached input vectors for post-processing diagnostics.

Group-sparse regularized GID#

group_sparse is a separate GID mode from the greedy lag-weight penalties described above. Lag-weight penalties modify which candidate lag the greedy iteration tries next. group_sparse instead solves one regularized model for the full sparse impulse response, the lag-indexed spike train returned by sparse_output. Its group-lasso penalty favors shared support, meaning the same nonzero lag samples, across the three components while still allowing different component amplitudes.

Relation to classic greedy GID#

group_sparse estimates the full sparse impulse response by solving an explicit regularized objective instead of selecting one spike per iteration. This makes it a different estimator, not a strict upgrade. Classic greedy GID remains useful for diagnostics, direct per-spike stopping behavior, and cases where close arrivals need to remain separate.

Both time-domain and frequency-domain engines still build a GID inverse operator and compact actual-output/resolution kernel. In the current implementation deconvolution_type group_sparse uses the ns_gid inverse branch internally, so the inverse-filtered data and resolution kernel inherit NS-GID’s gain cap, noise-spectrum damping, and optional reliability taper. After those quantities are built, the two engines use the same group-sparse objective.

Model and solver#

Let \(\mathbf{y}_c\) be component \(c\) of the inverse-filtered data, and let \(R\) be the compact, peak-normalized resolution-kernel convolution operator used by GID residual updates. The unknown sparse impulse response is a matrix \(H\), where each lag sample \(j\) has a three-component group and each \(\mathbf{h}_c\) denotes the lag series for component \(c\).

\[\mathbf{h}_j = \left(h_{E,j}, h_{N,j}, h_{Z,j}\right)^T .\]

The implemented model is

\[\hat{H} = \operatorname*{arg\,min}_{H} \frac{1}{2}\sum_{c \in \{E,N,Z\}} \left\|R \mathbf{h}_c - \mathbf{y}_c\right\|_2^2 + \lambda \sum_j \left\|\mathbf{h}_j\right\|_2 .\]

The first term asks the sparse impulse response, convolved with the actual GID resolution kernel, to explain the inverse-filtered data. The second term is an \(\ell_{2,1}\) group-lasso penalty. It shrinks all three component amplitudes at a lag as one group, favoring shared arrival support.

The solver uses proximal gradient iterations. Each iteration takes a gradient step on the data-misfit term and then applies group soft thresholding,

\[\mathbf{h}^{(k+1)}_j = \max\left(0, 1 - \frac{\tau \lambda}{\|\mathbf{z}_j\|_2}\right) \mathbf{z}_j .\]

The multiplier is interpreted as zero when \(\|\mathbf{z}_j\|_2 = 0\).

The implementation uses a conservative step bound derived from the finite impulse response (FIR) resolution kernel and stops when the objective’s relative change is below group_sparse_tolerance or when group_sparse_max_iterations is reached. After support is selected, the engine runs the same final amplitude refit used by classic greedy GID and recomputes the reported residual metrics.

At convergence, the group-lasso Karush-Kuhn-Tucker conditions clarify the role of group_sparse_lambda. For inactive lags,

\[\left\|R_j^T \mathbf{r}\right\|_2 \le \lambda ,\]

where \(\mathbf{r}=R\hat{H}-\mathbf{y}\) stacks the residuals for all three components and \(R_j\) denotes the resolution-kernel operator block associated with lag \(j\). For active lags,

\[R_j^T \mathbf{r} = -\lambda \frac{\hat{\mathbf{h}}_j} {\left\|\hat{\mathbf{h}}_j\right\|_2}.\]

Thus group_sparse_lambda is the evidence threshold inside the regularized objective. It is not, by itself, the threshold used to decide which nonzero coefficients are exported as sparse arrivals.

Regularization and support reporting#

group_sparse_lambda controls shrinkage inside the objective. A positive value is used directly. The default value 0.0 selects an automatic noise-scaled value. Let \(\lambda_{\mathrm{scale}}\) denote group_sparse_lambda_scale and let \(\eta\) denote the automatic noise reference. Then

\[\begin{split}\lambda = \lambda_{\mathrm{scale}}\eta, \qquad \eta = \begin{cases} \eta_{\mathrm{noise}}, & \mathrm{available}, \\ 0.02\,\|\mathbf{y}\|_\infty, & \mathrm{otherwise}. \end{cases}\end{split}\]

Here \(\eta_{\mathrm{noise}}\) is the NS-GID inverse-filtered noise threshold when that estimate is available. It is the same final threshold reported as group_sparse_noise_threshold: empirical vector quantile when enabled, guarded by the scalar-component sigma threshold described above.

Since the zero-lag normalization policy raw_zero_lag_normalized_v2, both an explicit group_sparse_lambda and a fixed group_sparse_active_threshold are interpreted in this normalized inverse domain. They are not numerically interchangeable with values calibrated against pre-normalization leaf output. Automatic lambda and quantile-based support thresholds migrate with the noise/coefficient scale; users carrying a fixed legacy value should retune it and retain the reported *_used QC fields with the result.

The exported sparse support has a second adaptive decision rule. The regularized coefficient field can contain tiny numerical coefficients or small clustered sidelobe coefficients. Let \(a_0\) denote group_sparse_active_threshold, \(s_a\) denote group_sparse_active_threshold_scale, and \(q\) denote group_sparse_active_threshold_quantile. The default support threshold is

\[\tau_{\mathrm{support}} = \max\left( a_0, s_a Q_q\left(\left\{\|\hat{\mathbf{h}}_j\|_2\right\}\right) \right),\]

with distributed defaults group_sparse_active_threshold=0.02, group_sparse_active_threshold_scale=1.0, and group_sparse_active_threshold_quantile=0.90.

This threshold is the group-sparse analogue of the adaptive part of adaptive_memory. The adaptation is global and distributional rather than a per-spike memory update, but it lets the support decision follow the solved coefficient field instead of a hand-picked constant. group_sparse_lambda is therefore the objective shrinkage strength, while group_sparse_active_threshold_used is the actual cutoff applied to sparse_output and the final refit model.

If the proximal solve is already sparse, most lag groups are zero and the absolute floor controls the exported support. If coherent leakage fills many lag groups, the upper-tail quantile rises and prunes the leaked coefficient field. Benchmark threshold sweeps are therefore important: disabling the adaptive upper-tail cutoff can leave many retained coefficients, while an overly large fixed floor can remove weak converted phases.

QC and interpretation#

QCMetrics records group_sparse_lambda_requested, group_sparse_lambda_scale, group_sparse_lambda_used, group_sparse_noise_threshold, group_sparse_peak_threshold_empirical, group_sparse_peak_threshold_sigma, group_sparse_use_empirical_noise_threshold, group_sparse_peak_threshold_controlling_term, and the group_sparse_noise_component_* requested/used-scale fields, group_sparse_active_threshold, group_sparse_active_threshold_scale, group_sparse_active_threshold_quantile, group_sparse_active_threshold_quantile_value, group_sparse_active_threshold_used, and group_sparse_active_groups. It also records group_sparse_iterations and group_sparse_converged for the proximal solve, plus group_sparse_debiased_objective_final and group_sparse_debiased_fractional_improvement_final for the amplitude refit on the selected support. These fields let downstream QC distinguish the regularized solve from the support-reporting and debiased-amplitude layers.

These are not optional display fields. They are part of the QC subdocument that should be saved with the deconvolved datum and used by downstream screening jobs to reject unstable or poorly constrained receiver functions. When deconvolution_type group_sparse is active, the group-sparse fields describe the solve that produced the output. When another GID method is active, the group_sparse_* fields are absent. The QC document is not a fixed table; method-specific keys are recorded only when that method or internal method component actually ran.

group_sparse_active_groups is a count of retained coefficient groups, not a resolved geologic arrival count. Validation plots and tests use separate detection metrics for arrival-level comparisons.

Because this mode uses the NS-GID inverse branch internally, it also records group-sparse-prefixed inverse QC fields such as group_sparse_inverse_gain_max_actual, group_sparse_inverse_noise_amplification, and group_sparse_inverse_effective_bandwidth_fraction.

The careful benchmark claim is narrow: in the supplied synthetic validation fixtures, group-sparse solving with the adaptive support rule can improve sparse-support recovery and false-positive control. It should not be presented as universally better than classic greedy GID, including adaptive_memory penalty runs. Adaptive-memory GID may still produce a lower final residual norm or preserve dense close arrivals better.

Validation and QC workflow#

The validation tests can write optional figures that are useful for learning how the deconvolution operators behave on controlled synthetic examples. These figures are an audit aid, not a field-data benchmark. The pytest assertions remain the source of pass/fail behavior.

Run the optional validation plots#

To write the default diagnostic figures, run:

python -m pytest \
    python/tests/algorithms/test_decon_algorithm_validation.py \
    --decon-validation-plots \
    --decon-validation-plot-dir /tmp/mspass-decon-validation-plots

Use --decon-validation-noise-scale when you want extra stress plots at a different synthetic noise amplitude:

python -m pytest \
    python/tests/algorithms/test_decon_algorithm_validation.py \
    --decon-validation-plots \
    --decon-validation-plot-dir /tmp/mspass-decon-validation-plots-noise-003 \
    --decon-validation-noise-scale 0.03

--decon-validation-noise-scale is the Gaussian noise amplitude before the synthetic coloring filter. It is not an SNR value and it is not a percent-noise setting. Changing it changes the optional figures used for visual inspection; the core numerical checks still use fixed reproducible validation cases.

What the figures show#

complex_colored_validation_wavelet.png shows the notched source wavelet, the noisy convolved three-component data, and the known sparse impulse response used as synthetic truth.
complex_colored_scalar_methods.png and scalar_noise_<scale>_stress_spike_results.png compare scalar and CNR receiver-function outputs. These traces are shaped receiver functions, not sparse impulse responses.
TimeDomainGIDDecon_* and FrequencyDomainGIDDecon_* overlays compare GID outputs for the time-domain and frequency-domain engines. Files ending in _sparse_results.png show raw sparse impulse responses from sparse_output and are the most direct visual comparison to the known sparse truth.
external_wavelet_all_methods.png compares methods when a prepared external wavelet is supplied. The companion external_wavelet_all_methods_display_filtered.png applies a common plotting-only display filter; that filter is not part of the algorithms or the pass/fail assertions.

Reading noise-scale plots#

Start with the setup figure. The data panel is normalized for display, so use the noise-scale label, pre-event RMS, signal RMS, and waveform shape together. A larger noise scale should make pre-event and inter-arrival fluctuations more visible, but it should not be read as a calibrated prediction of field-data performance.

If a high-noise diagnostic run marks a method as failed or omits that method’s output trace, treat that as stress-run context. It is a reason to inspect the matching QC metadata, not by itself a benchmark result.

Reading sparse-support plots#

For GID methods, compare the raw sparse plots to the true impulse response. The main questions are whether expected arrivals are present near the correct lag, whether signs are plausible, and whether the result contains many extra isolated or clustered spikes.

For group_sparse labels, lam is the regularization value used, thr is the support threshold used after the solve, k is the number of retained coefficient groups, one three-component lag group per lag sample, and it is the number of proximal iterations. k is not automatically a geologic arrival count.

Reading penalty comparison plots#

The penalty comparison figures should be read as a three-part diagnostic:

The shaped-result plot shows the receiver-function traces users would normally inspect.
The sparse-result plot shows how each penalty changed the accepted sparse support.
The lag-weight plot shows the final penalty state. Values near 1.0 are effectively unpenalized; lower values mark lag neighborhoods that were downweighted after earlier accepted spikes.

The labels report the penalty function, scale factor, effective width, and, for adaptive_memory, the last confidence and decay values. These plots can show whether a penalty reduced repeated picking near the same arrival, but they should not be summarized as one penalty universally beating another. Prefer wording such as “in this synthetic stress case” or “for this configured noise level.”

Inspecting QC metadata after deconvolution#

After running deconvolution, inspect the QC subdocument saved on the output datum. RFdecon uses RFdecon_properties by default. The direct GID wrappers use TimeDomainGIDDecon_properties and FrequencyDomainGIDDecon_properties by default. If you set a custom QCdocument_key or QCdata_key, inspect that key instead.

Useful first-pass fields are:

decon_operator and decon_processed confirm which engine ran and whether it processed the datum.
decon_window_start, decon_window_end, noise_window_start, and noise_window_end confirm that the analysis and noise windows match the intended workflow.
residual_rms_final_fraction and residual_rms_reduction_fraction are within-configuration fit metrics: final/initial residual RMS and its complement. Compare them across records only when the method, leaf operator, analysis window, preprocessing, and output shaping are the same. Different leaf bandwidths and inverse domains make them unsuitable as a standalone cross-algorithm ranking. Smaller final fraction is a larger fit reduction, but is not automatically better if the sparse support becomes physically implausible. Consult residual_rms_fraction_valid first: an unprocessed result or a zero initial residual has undefined fractional metrics (stored as NaN), not a perfect fit. residual_rms_initial and residual_rms_final normalize only for analysis-window length and therefore retain amplitude units. residual_energy_final_fraction and residual_energy_reduction_fraction are respectively the square of the RMS fraction and one minus that square; use them when a policy needs energy, rather than amplitude, reduction. residual_L2_initial and residual_L2_final remain legacy raw norms for within-record debugging only.
gid_converged, gid_stop_reason, gid_iterations, and gid_number_spikes summarize the GID iteration outcome.
gid_penalty_function, gid_penalty_effective_width, lag_weight_L2_final, and the gid_penalty_* adaptive-memory fields explain how the lag-weight penalty affected candidate selection. These fields are present for greedy GID modes, not for deconvolution_type group_sparse.
group_sparse_enabled, group_sparse_converged, group_sparse_iterations, group_sparse_lambda_used, group_sparse_active_threshold_used, group_sparse_active_groups, group_sparse_objective_initial, and group_sparse_objective_final summarize the regularized group-sparse solve and the exported sparse support decision. group_sparse_objective_final is the solver’s final regularized objective after thresholding the proximal solution. The amplitude refit is reported separately as group_sparse_debiased_objective_final and group_sparse_debiased_fractional_improvement_final because it is the same regularized objective evaluated after the debiased amplitude update on the selected support. These fields are present only for deconvolution_type group_sparse.
ns_gid_gain_max_actual, ns_gid_noise_amplification, and ns_gid_effective_bandwidth_fraction help audit inverse-operator stability for deconvolution_type ns_gid.
For ns_gid, ns_gid_last_selected_candidate_lag_samples, ns_gid_last_selected_candidate_lag_weight, and ns_gid_last_selected_candidate_weighted_amplitude describe the candidate selected after applying the lag penalty. The ns_gid_max_raw_candidate_* and ns_gid_last_scan_raw_significant_candidate_remaining fields describe the final iteration scan before amplitude refitting, not necessarily the terminal residual. The ns_gid_final_scan_* fields distinguish the off-support raw maximum from the penalized ordered trial list and expose significant_candidate_count, the best trial, and whether an off-support candidate clears both the significance and strict fractional-improvement gates. Existing-support amplitude is reported separately, so an unstable ridge refit cannot be hidden by excluding that support from candidate insertion. ns_gid_provisional_stop_reason_before_final_refit preserves the pre-refit decision; ns_gid_stop_reason is reclassified from the final residual or resumes the bounded iteration. Iteration-trace stop conditions remain historical and are not overwritten by that audit. ns_gid_noise_samples_at_or_above_peak_threshold and ns_gid_initial_stationary_null_expected_noise_exceedances are an initial stationary-null plug-in estimate, not a terminal QC metric. The latter is the observed noise exceedance rate multiplied by ns_gid_initial_stationary_null_search_lag_count (the actual initial positive-weight, complete-kernel search lags), not by the full analysis window length. These are tuning diagnostics, not hard QC cutoffs.
group_sparse_inverse_gain_max_actual, group_sparse_inverse_noise_amplification, and group_sparse_inverse_effective_bandwidth_fraction help audit the NS-GID inverse branch used internally by deconvolution_type group_sparse. These fields are absent for other GID modes.

Implementation and compatibility notes#

The C++/pybind multitaper classes are still registered as MultiTaperXcorDecon and MultiTaperSpecDivDecon for ABI and pickle compatibility. MultiTaperPowerXcorDecon and MultiTaperPowerSpecDivDecon are Python-level class aliases to those classes, not distinct runtime types. The old RF processor algorithm names MultiTaperXcor and MultiTaperSpecDiv remain available as deprecated compatibility aliases for the power-stabilized operators.

MultiTaperSpecDivDecon still exposes legacy all_inverse_wavelets, all_rfestimates, and all_actual_outputs methods. In the current implementation those compatibility methods normally contain one combined multitaper estimate rather than one independent product per taper.

Parameter files tuned against older tapered-numerator multitaper behavior should be retuned and revalidated. damping_factor is now interpreted in the power-domain denominator, so values should not be assumed numerically equivalent to historical Park-Levin-style implementations.