Large seismic waveform models can efficiently support earthquake detection, phase picking, and event classification by transforming waveforms into high-dimensional feature vectors. However, the seismological interpretation of these intermediate representations remains difficult to establish owing to the black-box nature of deep neural networks. This study investigates whether feature vectors extracted from the pre-trained DiTingSFM-0.1B model retain measurable waveform characteristics, catalog parameters, event-related information, and spatial sampling patterns. Unsupervised clustering is adopted as a probe for representations rather than a tool to define intrinsic earthquake categories.

The dataset comprises 3 473 three-component local-earthquake waveform samples selected from DiTing Dataset 2.0, corresponding to 3 333 independent natural earthquakes recorded at 714 stations between February 2022 and February 2023. All samples have epicentral distances less than 200 km. Epicenters lie between approximately （18°N−52°N, 73°E−132°E）, with local magnitudes ranging from M_L0.0−5.7. Each sample includes vertical, north–south, and east–west components sampled at 100 Hz. A 60 s waveform window is aligned relative to the P-wave arrival time. Prior to model inference, each trace is demeaned and normalized by its standard deviation computed over the full window. A 512-dimensional feature vector is extracted after global average pooling and prior to the task-specific classification head of the pre-trained DiTingSFM-0.1B model.

Four feature-input schemes were compared: raw model outputs, z-score-standardized features, L₂-normalized features, and z-score standardization followed by L₂ normalization. K-means＋＋ clustering was implemented for candidate numbers of clusters ranging from 2 to 10. The silhouette coefficient, Calinski–Harabasz index, and Davies–Bouldin index were adopted to assess clustering structure. The median Hopkins statistic was 0.757 （95% interval: 0.754−0.761）, indicating non-random aggregation within the feature space. Nevertheless, the gap statistic failed to identify a unique optimal number of clusters. The silhouette coefficient and Calinski–Harabasz index favoured a coarse two-cluster structure, while the Davies–Bouldin index supported further subdivision at k＝4. Therefore, k＝4 was adopted as an operational partition to compare waveform, spectral, catalog, and spatial characteristics, rather than interpreted as the unique count of intrinsic seismological classes.

The reproducibility of the k＝4 partition was assessed via repeated initialization, independent-event resampling, and alternative clustering algorithms. Across 100 independent initializations, the median adjusted Rand index and normalized mutual information were 0.992 and 0.985, respectively. Event-level resampling yielded corresponding median values of 0.979 and 0.963, with cluster-specific Jaccard retention rates ranging from 0.981 to 0.987. Principal-component K-means and spherical K-means yielded results broadly consistent with the original partition, whereas the Gaussian mixture model and HDBSCAN failed to recover the same structure. Accordingly, the k＝4 partition is reproducible within K-means and related distance-partitioning frameworks yet partially dependent on the choice of clustering algorithm. The four clusters contained 433, 1604, 748, and 688 waveform samples, respectively.

A station-aware nearest-neighbour analysis was performed to evaluate geographic consistency within the feature space. After removing additional records belonging to the same earthquake, 38.1% of all feature-space nearest neighbours originated from the same station, substantially higher than the median value of 0.37% for random pairs. This finding demonstrates pronounced station dependence within the feature space. The analysis was repeated after excluding records from both the same earthquake and the same station. The observed median epicentral distance between nearest neighbours was 881.7 km. A station-blocked permutation test, which preserved each station’s sample size and geographic sampling extent, yielded a null-distribution median of 889.1 km （95% interval: 873.7−902.1 km）. The observed distance was not significantly smaller than the null distribution （p＝0.203）.

A covariate-residualization analysis was further carried out to control for station identity, magnitude, epicentral distance, back azimuth, and recording time. The nearest-neighbour distance computed on residualized features was 1 197.3 km, compared with a station-blocked null median of 1 191.7 km; the difference was not statistically significant （p＝0.868）. Between station identity and cluster labels, the bias-corrected Cramér’s V was 0.540 and normalized mutual information was 0.182, indicating a non-negligible yet partial association. For 134 earthquakes recorded at multiple stations, 62.3% of the 146 cross-station record pairs were assigned to the same cluster, significantly exceeding the random baseline median of 28.8%. Station-grouped holdout validation also yielded an adjusted Rand index of 0.971±0.013 and normalized mutual information of 0.954±0.016. These results show that feature-space structure is strongly modulated by station and sampling conditions yet not solely determined by station identity. The representations retain measurable event-related waveform information, though nearest-neighbour analysis alone cannot independently prove that they encode propagation paths or regional medium structure.

Waveform envelopes, power spectral densities, catalog parameters, and waveform statistics were compared across the four clusters. This analysis adopted a processing workflow independent of the standardization applied to model inputs. Each trace of the waveform in digital counts was centered by subtracting its median and normalized by its maximum absolute value. Relative root-mean-square （RMS） values were computed for the full 60 s window, the P-S window, and the 10 s post-S-wave window. These dimensionless quantities describe waveform fluctuations relative to the trace peak and do not represent absolute ground-motion amplitude or seismic wave energy. Power spectral density was estimated via Welch’s method, and dominant frequency was defined as the frequency corresponding to the maximum average three-component power spectral density.

For statistical comparisons, a single record was retained per independent earthquake, and station-blocked permutations were used to construct null distributions. Back azimuth was parameterized using its sine and cosine components, and the Holm correction was applied across all statistical tests. Differences in back azimuth were minor, whereas magnitude and epicentral distance exhibited significant inter-cluster differences, with effect sizes of 0.336 and 0.069, respectively. Focal depth showed no significant inter-cluster difference. The full-window, P-S-window, and post-S-wave-window relative RMS yielded effect sizes of 0.361, 0.300, and 0.176, while dominant frequency had an effect size of 0.435.

Clusters C₀ and C₁ generally featured lower magnitudes, intermediate epicentral distances, relatively low dominant frequencies, and larger relative RMS values. Cluster C₂ was characterized by shorter epicentral distances, moderately larger magnitudes, smaller relative RMS values, and substantially higher dominant frequencies. Cluster C₃ typically exhibited larger epicentral distances and magnitudes, intermediate relative RMS values, and stronger spectral decay. These results demonstrate that this operational feature-space partition corresponds to measurable differences in peak-normalized waveform envelopes, relative temporal fluctuations, spectral structure, magnitude, and epicentral distance. Since instrument responses and station gains were not uniformly corrected, these statistics are interpreted only as relative waveform differences under the present processing conditions.

Six low-dimensional baseline models were constructed to examine whether the k＝4 partition could be reproduced using a small set of variables. These baselines comprised three single-variable cases （feature-vector norm, full-window relative RMS, dominant frequency）, alongside simple catalog-parameter, waveform-statistic, and combined baselines. Agreement with the 512-dimensional clustering remained limited: the adjusted Rand index did not exceed 0.226, and normalized mutual information did not exceed 0.292. Hence, none of the tested individual variables or low-dimensional combinations directly reproduced the four-cluster partition. Nonetheless, this low agreement does not prove that high-dimensional representations contain richer or more physically meaningful information; discrepancies may alternatively arise from unmeasured waveform attributes, station and instrument characteristics, high-dimensional noise, or feature-space geometry.

Cluster epicenters were further compared with topography and horizontal gradients derived from WGM2012 Bouguer gravity anomalies. The transition between Cluster C₀ and the remaining clusters lies broadly near the north-northeast-trending Greater Khingan–Taihang–Wuling gravity-gradient belt. This pattern suggests a plausible spatial correspondence between cluster distribution and regional geophysical background. However, station-blocked permutation and covariate-residualization analyses revealed no significant epicentral proximity after accounting for station-specific sampling biases. The association between clustering results, topography, and Bouguer gravity gradients should thus be treated as an exploratory spatial correspondence rather than direct evidence that the model identifies tectonic boundaries.

Overall, unsupervised clustering offers a powerful approach for probing feature vectors generated by large seismic waveform models. Representations from the pre-trained DiTingSFM-0.1B model reflect observable differences in magnitude, epicentral distance, peak-normalized waveform features, relative RMS, and dominant frequency. Their spatial organisation is also strongly linked to station composition and geographic sampling, and cross-station consistency indicates partial retention of event-related information. Since the records used in this study overlap in source properties and time window with a subset of the model’s pre-training corpus, these results characterise representations formed under the present data distribution and should not be interpreted as proof of generalisation to fully unseen data. Further validation using independent time periods, regions, and stations, combined with uniformly corrected instrument responses and explicit event–station propagation paths, is required to disentangle the respective contributions of earthquake source, propagation path, station response, instrumental characteristics, and sampling structure.