Modified Stockwell Frequency-Wavenumber Workflow for Surface-Wave Dispersion Representation

Abstract

# MSFK MATLAB implementation and reproducible validation This package accompanies **Modified Stockwell Frequency-Wavenumber Workflow for Surface-Wave Dispersion Representation** by Hanbing Ai, Jiangtao Li, Yunus Levent Ekinci, Yang Hu, Ling Ning, and Yingwei Yan. The MSFK algorithm was developed by **Hanbing Ai**. It selects the dominant complex coefficient of the standard Stockwell transform at each frequency and receiver and passes those coefficients to a spatial Fourier transform. The standard Stockwell transform is unchanged. MSFK requires neither group-velocity steering nor repeated slant stacking. ## Version and citation status - Prepared package version: **1.1.0**.- Preparation and validation date: **2 October 2026**.- Primary developer: Hanbing Ai, hanbingai@foxmail.com.- Validation environment: MATLAB R2024b, PCWIN64.- Version 1.1.0 DOI: https://doi.org/10.5281/zenodo.23121000.- Original version DOI: https://doi.org/10.5281/zenodo.22108446. The authors reserved the version 1.1.0 DOI in Zenodo on 3 October 2026. Zenodo registers this DOI when the record is published. The original version DOI remains the identifier of the earlier archive. Original source headers retain that historical DOI; use the version 1.1.0 citation below for the complete package and added validation. No publication date is asserted in this package. Recommended citation: Ai, H., Li, J., Ekinci, Y. L., Hu, Y., Ning, L., and Yan, Y. (2026). Modified Stockwell Frequency-Wavenumber Workflow for Surface-Wave Dispersion Representation (Version 1.1.0) [Computer software]. Zenodo. https://doi.org/10.5281/zenodo.23121000. ## Files and quick start The original implementation and example are retained without changes: - `msfk_dispersion.m`: MSFK dispersion-image implementation.- `msfk_stockwell_transform.m`: documented complex Stockwell transform.- `Stockwell_FK.m`: compatibility wrapper.- `test.m`: original synthetic example.- `OUTPUT_SG_DWM_20260323_300Hz.mat`: original synthetic input. Run the original example from the package directory: ```matlabtest``` The second-round validation is in `validation_R2/`. Its copies of the MSFK core and synthetic input are byte-identical to the original files above. This duplication allows the validation folder to run independently as journal supplementary code. ```matlabaddpath('validation_R2')run_round2_testscalculate_band_RMSE``` These commands write into `validation_R2/results/`. Example outputs are already supplied there; copy the package to a writable directory before rerunning it. `validation_R2/README.txt` describes the tests and their parameters. ## Additional validation and interpretation The package adds known-phase single-mode controls, overlapping and separated two-mode controls, an unnormalized input-spectrum and Parseval audit, Gaussian/Hann coefficient-window ablation, and five paired noise realizations at each noise level. The noise results compare each method with its own noise-free image. They measure normalized-image sensitivity, not velocity RMSE, repeated-human-pick uncertainty, or modal detection probability. The point audit retains the supplied picks and frozen mode assignments. Theoretical velocities are linearly interpolated from saved Shi et al. SASEM branches at each pick frequency, without extrapolation or bridging missing segments. No new solver call is made at each point. The calculator independently reconstructs those residuals and applies the same 5-80 Hz eligibility band to every method, with no residual-based rejection. Verification results: - Both 40 Hz controls peak at the prescribed 320 m/s on the spatial FFT lattice; maximum inter-receiver phase residuals are approximately 4.12e-14 and 1.35e-14 rad.- The independent Parseval check has relative error approximately 1.34e-16.- Independent re-interpolation agrees with cached residuals to a maximum difference of approximately 4.12e-13 m/s.- Separate portable control and eight-method ensemble runs reproduce the validation statistics; the band-table difference is below 1e-12 m/s. The controls support phase consistency for the stated narrowband cases and retention of concurrent contributions; weaker temporally separated arrivals can be suppressed. The supplied Figure 2 M5-M6 picks coexist at 99-102 Hz in a weak-input-energy interval. At 100 Hz, their velocities are 284.8 and 301.2 m/s, compared with SASEM reference values of 285.36 and 301.45 m/s. This supports local ridge separation and reference agreement, without validating every high-frequency feature. The M5 labels within the supplied mixed-mode file are nearest-admissible-branch assignments; the M6 labels are explicitly supplied. Assignments alone are not independent modal ground truth. The common 5-80 Hz table band is a conservative evaluation choice for these records, not a universal MSFK frequency limit. No inversion experiment is included. ### Publication palette `validation_R2/publication_colormap.mat` stores the exact 254-row RGB palette used by manuscript Figure 2. `format_round2_publication.m` applies it explicitly to each normalized spectral axes, with colour limits [0,1] and colour-bar ticks at 0, 0.5, and 1. This keeps the supplementary spectra consistent with the synthetic and field comparisons. The presentation correction does not change spectral arrays, theoretical curves, supplied picks, or numerical metrics. Grayscale gathers and line/error-bar plots retain their distinct encodings; the original example script is preserved unchanged. ## External comparator implementations To rerun the eight-method noise diagnostic: ```matlabrun_round2_ensemble(benchmarkDir, sfkDir)``` `benchmarkDir` must contain the original `phase_kai.m`, `FK_Dispersion.m`, `Radon_conjgrad.m`, and helpers; `sfkDir` must contain SFK and its dependencies. These third-party functions are not redistributed here. F-J and F-H use Kai Zhang's MATLAB `phase_kai.m` from [PassiveFW](https://github.com/jordankai/PassiveFW/tree/main/active%20seismic%20modeling%20by%20DWM). The supplied wrapper records the original calls: time-by-receiver data; offsets 3:1:62 m; velocities 50:2:600 m/s; dt = 0.0005 s; fmax = 150 Hz; df_factor = 1; vertical component `v`; kernels `bessel` and `hankel`. The additional ensemble uses fmax = 80 Hz. [CC-FJpy](https://github.com/ColinLii/CC-FJpy) provides related Python implementations; it was not timed in this study. SFK is obtained from the source associated with Serdyukov et al. (2019). Obtain comparator code from its original providers and comply with their terms. Restricted legacy `st.m` and `seidisp.m` are excluded and are not needed by the self-contained controls. ## Attribution and license Original source headers, authorship, creation dates, and copyright notices are preserved. Algorithm development is credited to Hanbing Ai; the added validation scripts document tests prepared for this study on 2 October 2026. See `LICENSE` and `THIRD_PARTY_NOTICES.txt` for the software-license scope and exclusions. No new ownership or license is asserted over third-party software or field datasets. No field seismic records or published field curves are redistributed here. ## References Ai, H., Li, J., Ekinci, Y. L., Hu, Y., Ning, L., and Yan, Y. (2026). Modified Stockwell Transform-Driven Frequency-Wavenumber Method for Reliable Surface-Wave Dispersion Energy Representation [Computer software]. Zenodo. https://doi.org/10.5281/zenodo.22108446. Original archive title, retained for historical citation. Stockwell RG, Mansinha L and Lowe RP (1996). Localization of the complex spectrum: The S transform. IEEE Trans Signal Process 44(4):998-1001. https://doi.org/10.1109/78.492555. Shi C, Ren H, Li Z and Chen X (2022). Calculation of normal and leaky modes for horizontal stratified models based on a semi-analytical spectral element method. Geophys J Int 230(3):1928-1947. https://doi.org/10.1093/gji/ggac163. Wang J, Wu G and Chen X (2019). Frequency-Bessel transform method for effective imaging of higher-mode Rayleigh dispersion curves from ambient seismic noise data. J Geophys Res Solid Earth 124(4):3708-3723. https://doi.org/10.1029/2018JB016595. Yang Z, Sun YC, Zhang D, Han P and Chen X (2024). A frequency-Hankel transform method to extract multimodal Rayleigh wave dispersion spectra from active and passive source surface wave data. Geophysics 89(2):KS69-KS81. https://doi.org/10.1190/geo2023-0189.1. Xi C, Xia J, Mi B, Dai T, Liu Y and Ning L (2021). Modified frequency-Bessel transform method for dispersion imaging of Rayleigh waves from ambient seismic noise. Geophys J Int 225(2):1271-1280. https://doi.org/10.1093/gji/ggab008. Li Z, Zhou J, Wu G, Wang J, Zhang G, Dong S, Pan L, Yang Z, Gao L, Ma Q, Ren H and Chen X (2021). CC-FJpy: A Python Package for Extracting Overtone Surface-Wave Dispersion from Seismic Ambient-Noise Cross Correlation. Seismol Res Lett 92(5):3179-3186. https://doi.org/10.1785/0220210042. Serdyukov AS, Yablokov AV, Duchkov AA, Azarov AA and Baranov VD (2019). Slant f-k transform of multichannel seismic surface wave data. Geophysics 84(1):A19-A24. https://doi.org/10.1190/geo2018-0430.1. Zhang K (2021). Active and passive seismic full-wavefield modeling based on discrete wavenumber method (Version 1.0) [Computer software]. Zenodo. https://doi.org/10.5281/zenodo.5559375.