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:

  1. Skim Essential terms so the method descriptions use familiar words.

  2. Read Wrappers and engines, Common conventions, and Choosing an operator to understand the processing workflow.

  3. Read the method section for the operator you plan to use.

  4. 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 as the interval to analyze and return. 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.

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.

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.

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. Their default parameter files do not apply a GID-style output shaping wavelet (shaping_wavelet_type none). If an older scalar parameter file configures an output filter, treat it as a legacy display or bandlimiting filter.

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.

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.

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.

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_L2_initial, and residual_L2_final.

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.

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_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 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.

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. The implemented selection statistic is

\[q^{(k)}_j = \left(\ell^{(k)}_j\right)^2 \left\|\mathbf{a}^{(k)}_j\right\|_2^2, \qquad j_* = \operatorname*{arg\,max}_j q^{(k)}_j .\]

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:

  1. Which lags are in the same ambiguous state as the accepted spike?

  2. How reliable is the accepted spike as an explanation of the residual?

  3. 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.

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_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_L2_initial and residual_L2_final check whether the fitted model reduced residual energy. Smaller is not automatically better if the sparse support becomes physically implausible.

  • 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.

  • 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.