Abstract
The thickness of sediments above bedrock controls seismic site response and groundmotion amplification, but is difficult to map at high resolution in urban areas with sparse boreholes. This control is captured by the shallow equivalent thickness—the effective sediment-layer thickness that reproduces the observed near-surface resonance, rather than the purely geometric depth to bedrock. Conventional horizontal-to-vertical (H/V) inversions assume a single homogeneous layer, while purely data-driven models lack physical constraints and may produce physically inconsistent predictions under limited samples. This study proposes a physics-constrained random forest (RF) machine learning model for spatial prediction of the shallow equivalent thickness from threecomponent microtremor data, guided by quarter-wavelength theory. At 101 sites and 11 co-located boreholes in the Youjiang river-valley basin, the H/V fundamental frequency and amplification factor, three-component energy ratios, a waveform factor and site coordinates form a multi-source feature vector. A site-specific fundamental-frequency–thickness power law is calibrated against boreholes (R2 = 0.96; mean error 5.32%) and embedded as the dominant feature of the random forest. Under 5- fold cross-validation the RF achieves R2 = 0.857, RMSE = 1.81 m and MAE = 1.02 m, outperforming support vector regression, radial basis function and k-nearest neighbour methods (R2 = 0.726–0.778) and traditional interpolation (R2 < 0.16). Predicted thickness lies between 10 and 26 m and varies inversely with the fundamental frequency (2.29–6.71 Hz), consistent with quarter-wavelength theory; removing it reduces R2 by 0.283 and inflates RMSE and MAE by 72.4% and 122.5%. The framework validates that integrating physical constraints with data-driven learning enables high-precision, interpretable urban-scale subsurface characterization under limited sample conditions.
1 Introduction
The rapid acceleration of urbanization has led to increasing spatial heterogeneity in shallow subsurface structures in urban areas. The configuration of near-surface geological media exerts a fundamental control on site dynamic characteristics, as well as on the propagation and amplification of seismic waves, making it a key factor in determining seismic response. Previous studies have demonstrated that the site fundamental frequency and its associated amplification effect are critical parameters for characterizing site response, both of which are closely related to the thickness of shallow sedimentary layers and the shear-wave velocity structure. Therefore, accurate estimation of shallow structural parameters that reflect site dynamic behavior is essential for urban seismic hazard assessment and earthquake-resistant engineering design. Among these parameters, the shallow equivalent thickness has emerged as an effective proxy for comprehensively characterizing the influence of near-surface sedimentary structures on ground motion (Zhong et al., 2022; Xu et al., 2023; ; Zhu et al., 2024).
Currently, the determination of shallow subsurface structural parameters primarily relies on methods such as borehole and seismic explorations and estimations using empirical formulas. While borehole data can provide reliable stratigraphic information, it is costly and has limited spatial coverage, making it difficult to meet the requirements for detailed mapping at the urban scale; Seismic exploration methods, by contrast, can provide continuous information on subsurface structures; however, they are highly susceptible to environmental noise, often requiring advanced machine learning and signal-processing techniques for noise suppression (; Zhan et al., 2024). In addition, their application in complex urban environments is frequently constrained by construction conditions (); Empirical methods based on the quarter-wavelength theory offer the advantage of simple calculations, but they typically rely on the assumption of a homogeneous medium and are difficult to apply to complex geological conditions with significant spatial heterogeneity (; ).
Relevant studies indicate that in complex urban sites, shallow subsurface structures often exhibit significant lateral heterogeneity, which further limits the applicability of traditional methods. In recent years, the H/V method, based on microtremor observations, has been widely applied to extract site fundamental frequencies and investigate site effects, providing an effective means for indirectly estimating subsurface structural parameters. However, traditional approaches often rely on a single parameter or simple empirical relationships, making it difficult to fully characterize nonlinear features under complex conditions, thereby limiting their accuracy and stability in urban-scale applications (; ; ; Xu and Wang, 2021; ; ; ).
With the rapid development of machine learning methods in the field of geophysics, these techniques have demonstrated significant advantages in modeling complex nonlinear relationships and integrating multi-source data. They have been widely applied in areas such as site classification, seismic motion prediction, and subsurface structure inversion, and have demonstrated potential to outperform traditional methods in modeling complex systems (Yuan, 2026). Among these, ensemble learning models, such as random forests, possess excellent robustness and generalization capabilities, offering a new technical approach for the high-precision modeling of complex subsurface structure parameters. However, existing research still has certain limitations: most methods rely on a single feature or lack physical constraints, resulting in insufficient model interpretability and limited generalization capabilities of the model. Simultaneously, traditional methods for acquiring subsurface structure data rely on dense and costly borehole data, making it difficult to meet the demand for rapid assessments at the urban scale (; ; ; ; ; ; ; Yang et al., 2024; ; ; Zou et al., 2021; Zuo and Carranza, 2023; ).
Furthermore, commonly used spatial interpolation methods are essentially linear or locally weighted models, making it difficult to capture nonlinear spatial variation characteristics under complex geological conditions, which limits their effectiveness. Against this backdrop, urban site issues in typical complex geological environments have gradually attracted increasing research attention. Taking the city of Baise in Guangxi, China, a typical river valley city, as an example, the influence of fluvial deposition has led to the development of shallow sedimentary structures with significant spatial variability. The subsurface structure exhibits pronounced lateral heterogeneity and is characterized by a certain degree of seismic hazards. Therefore, conducting research in such environments is of great significance for revealing the dynamic response characteristics of sites under complex geological conditions and verifying the applicability of relevant methods (Zhou et al., 2022; Zhang et al., 2025).
The remainder of this paper is organized as follows. Section 2 describes the study area, including its geographical, geological, and tectonic settings. Section 3 details the microtremor data acquisition, the H/V spectral ratio analysis, the multi-source feature construction, the physics-constrained random forest regression model, and the model evaluation framework. Section 4 presents the spatial distribution of the dynamic site parameters, the calibration of the empirical relationship between the fundamental frequency and the equivalent thickness, the borehole verification, and the predictive performance of the proposed model. Section 5 discusses the synergistic mechanism underlying the superior performance of the proposed model, the methodological innovations, the convergence and uncertainty of the prediction results, and the engineering applications and future perspectives. Section 6 concludes the paper.
2 Study area
2.1 Geographic setting
The study area is located in the Youjiang District of Baise City, Guangxi, China. Multiple observation points were established within the study area to conduct microseismic observations and fundamental frequency analyses. The spatial distributions of the variables are shown in Figure 1. Geological data indicate that the study area is situated on the first terrace of the Youjiang River Valley at an elevation of approximately 125 m above sea level, with generally flat topography. The Quaternary sedimentary layers in the study area are well developed, and the shallow strata consist primarily of plant fill, alluvial clay, silty clay, and gravel layers overlying the Tertiary mudstone bedrock. This stratigraphic sequence reflects the characteristics of a typical fluvial-alluvial depositional environment. This area belongs to a typical river valley urban landform unit controlled by fluvial alluvial processes. The thickness and physical properties of the sedimentary layers vary significantly across the study area, providing a geological context for the formation of spatial heterogeneity in the shallow subsurface structure.
FIGURE 1
According to the “Seismic Motion Parameter Zoning Map of China” (GB-18306-2015), the characteristic period of the design response spectrum for the study area is with a peak ground acceleration of 0.10 g. The corresponding basic seismic intensity was , classifying the area as a region with moderate seismic hazard. In river valley sedimentary environments, shallow subsurface structures typically exhibit significant spatial heterogeneity owing to spatial variations in sediment thickness and physical properties. This heterogeneity in the subsurface structure may significantly affect the site dynamic characteristics and lead to substantial spatial variations in the fundamental frequency. Therefore, conducting dense microseismic observations in this region is of great significance for revealing the spatial distribution characteristics of the fundamental frequency.
2.2 Urban environment
The study area is a typical urban built environment characterized by a dense distribution of buildings, with urban functions primarily centered on residential, commercial, and public-service facilities. The buildings in the area are predominantly reinforced concrete structures, mostly mid- and high-rise buildings. The area is also equipped with a relatively well-developed urban road network and public infrastructure, which creates a complex urban spatial environment. The spatial distribution of the microseismic observation points in the study area is shown in Figure 2.
FIGURE 2
In urban areas, human activities such as road traffic, pedestrian movement, and machinery operation continuously generate environmental vibration signals. These anthropogenic vibration sources constitute a major source of microseismic signals in urban environments and form a relatively stable background vibration field. This type of persistent environmental vibration provides stable data conditions for passive-source observations. The relatively consistent environmental vibration conditions offer a solid data foundation for passive-source microseismic observations, making it possible to extract the site-fundamental frequency using the H/V spectral ratio method. Building on this foundation, combined with dense observation data and data-driven methods, further analyses of site dynamic parameters can be conducted.
Tectonically, the study area is located within the Baise Basin, a NW–SE-elongated Cenozoic continental rift basin developed along the southwestern margin of the South China Block, within the western Youjiang Fold Belt that was shaped by the Indosinian Orogeny during the convergence of the South China and Indochina blocks. The basin is structurally controlled by NW–SE-striking boundary faults, with the Baise–Bama Fault Zone (part of the broader Youjiang Fault System) serving as the principal basin-controlling and seismogenic structure, while subordinate NE-trending faults further partition the basin into secondary sub-blocks. This tectonic framework, together with fluvial channel migration and differential subsidence during the Quaternary, accounts for the substantial lateral variation in sediment thickness and physical properties across the study area, providing the geological basis for the pronounced shallow subsurface heterogeneity investigated in this study.
3 Data and methods
This study employed a three-component broadband intelligent seismometer (IGU-BD3C-5) to conduct ambient microtremor observations. A total of 101 observation points were deployed over an area of approximately 2 km × 3 km, with a relatively uniform grid distribution and spacing of about 200–300 m between adjacent points. Uniform acquisition parameters were applied across all observation points, with a sampling interval of 0.25 m (4,000 Hz). GNSS positioning and GPS time synchronization were used to ensure temporal consistency among different instruments, and all sensors were oriented northward. During the observation period, the instruments were placed directly on the ground with good coupling conditions to ensure stable signal acquisition. The seismometer exhibits excellent low-frequency response characteristics, with a natural frequency range of 0.16–0.24 Hz, sensitivity of 160–240 and impedance of 1,665–2,035 Ω. The gain of all channels was uniformly set to 6 dB.
Each observation lasted 20–30 min to ensure stable and representative ambient vibration signals. Observations were conducted during periods of relatively stable ambient conditions to minimize transient disturbances. The collected data underwent initial quality control to remove anomalous records, after which the valid data were segmented into multiple stable time windows. Signal processing was then performed for each window, including zero-padding, bandpass filtering, Fast Fourier Transform (FFT), energy normalization, and spectral smoothing. Considering that the original sampling frequency (4,000 Hz) was much higher than the frequency range of interest, bandpass filtering (0.05–30 Hz) using a zero-phase second-order Butterworth filter was applied prior to spectral analysis to extract the target frequency band and suppress out-of-band noise, thereby meeting the requirements for microtremor analysis in the frequency range of 0.5–20 Hz.
To minimize the influence of transient noise, this study adopted a multi-stage screening procedure. Observations were conducted during periods of low anthropogenic activity, and records with obvious disturbances were re-acquired in the field. Each raw three-component record was detrended, demeaned, and visually inspected to discard transient spikes or instrument anomalies. The remaining records were segmented into non-overlapping 40-s windows with 0% overlap, and a short-term-average to long-term-average (STA/LTA) criterion was applied; windows with STA/LTA ratios exceeding 2.0 on any component were excluded. The H/V curves of the retained windows were then checked for consistency following the SESAME (2004) criteria, and windows deviating significantly from the site average were further rejected. Only the stable windows passing all the above criteria were used for FFT analysis with a 5% cosine (Tukey) taper applied to each window and bandpass filtering (0.05–30 Hz), Konno–Ohmachi spectral smoothing (b = 40), and final H/V averaging. The complete processing workflow is summarized in Figure 3.
FIGURE 3
3.1 H/V spectral ratio analysis
To extract the fundamental frequencies at each observation point, this study employed the H/V spectral ratio method to analyze the microseismic data. First, the continuous microseismic records were divided into multiple time windows, and a fast Fourier transform was applied to each component signal to obtain the spectrum. The H/V spectral ratio is defined as the ratio of the spectral amplitude of the horizontal component to that of the vertical component, where the horizontal spectral amplitude is composed of two horizontal components, north–south and east–west (), and is calculated as follows (Equation 1):where , and represent the spectral amplitudes of the north-south, east-west, and vertical components at a frequency .
The fundamental frequency f0, defined as the frequency at which the H/V spectral ratio reaches its maximum (Equation 2), reflects the resonance characteristics of the site and serves as a key parameter governing seismic amplification. The amplification factor is defined as the H/V spectral amplitude at the fundamental frequency as follows (Equation 3): This characterizes the vibration response of the site at the resonance frequency and can be used to reflect the amplification effect of shallow strata on seismic waves.
The quarter-wavelength theory is commonly used to describe the relationship between the fundamental frequency and the thickness of sedimentary layers. According to classical theory, for SH waves, the fundamental resonance frequency of the sedimentary layer can be expressed as : represents the equivalent shear wave velocity of the deposit layer, and represents the equivalent thickness calculated based on the quarter-wavelength theory.
Urban subsurface structures typically exhibit significant spatial heterogeneity, making it challenging to accurately reflect true geological characteristics with a single, constant shear wave velocity. Therefore, an equivalent shear wave velocity constrained by borehole data is introduced (Equation 5):
Here, represents the actual stratigraphic thickness revealed by the borehole, while denotes the equivalent shear wave velocity derived by inverting the fundamental frequency and the borehole-revealed thickness. This parameter is not a direct measurement but rather an equivalent parameter that reflects the integrated dynamic characteristics of the local site. Through this calibration process, shear wave velocities more consistent with actual geological conditions are obtained, thereby enhancing the reliability and applicability of sedimentary layer thickness estimates based on the fundamental frequency. In this step, a coupling constraint is established between the physical model and the observed data.
Under complex subsurface conditions, the thickness of sedimentary layers is difficult to accurately represent using an ideal homogeneous model. Therefore, the concept of equivalent thickness is introduced to characterize the effective thickness of sedimentary layers governing the site’s dynamic response. It should be noted that the equivalent thickness does not strictly correspond to the geometric thickness of the strata, but rather represents an equivalent parameter reflecting the site’s dynamic characteristics (). On this basis, considering that the site’s fundamental frequency and the equivalent thickness typically exhibit a power-law relationship in observations, this relationship can be expressed as:
In this equation, represents the equivalent site thickness, while and are empirical constants closely related to the subsurface conditions of the study area. Within the study area, a finite number of key borehole locations are selected to spatially cover the range of fundamental frequency variations. Three-component microtremor observation points are typically deployed in proximity to borehole control points. By calculating the H/V spectral ratio curves at these locations and extracting the fundamental peak frequencies, a data-driven relationship between the equivalent thickness and the site’s fundamental frequency is established.
Taking the natural logarithm of both sides of Equation 6 yields a linear form (Equation 7):
Taking as the independent variable and as the dependent variable, linear regression is performed using the least squares method. The regression slope corresponds to parameter , while the intercept corresponds to , thereby determining the specific values of the empirical parameters and . Under the constraint of borehole data, this method can establish an empirical relationship between sedimentary layer thickness and the fundamental frequency in the study area.
3.2 Multi-source feature construction
In this study, to comprehensively characterize the site dynamic response under complex urban conditions and improve the model’s capability to represent shallow subsurface structures, a multi-source feature system was constructed based on three-component microtremor observations. In addition to the derived from the H/V spectral ratio, three-component energy features were introduced as auxiliary variables to reflect the relative contributions of vibrational energy in different directions.
To improve the accuracy of equivalent thickness prediction in complex urban environments, multi-source features were extracted from the three-component microtremor data to construct an input feature vector (Equation 8):
In this equation, denotes the multidimensional input feature vector used for training the machine learning model, which integrates the spectral characteristics, energy distribution features, and spatial location information of the observation points. and represent the longitude and latitude of the observation points, respectively, and are used to characterize spatial location features, thereby assisting the model in capturing the spatial distribution patterns of site parameters. and , characterize the site’s resonance properties. and denote the normalized proportions of vibrational energy in the vertical, north–south, and east–west directions, respectively, as defined in Equation 10.
The energy characteristics of the three components are obtained by integrating the Fourier amplitude spectra of each component signal, and are defined as Equation 9:
Here, represents the energy of the -th component within the specified frequency band, where denotes the vertical, north–south, and east–west components, respectively. denotes the Fourier amplitude spectrum of the corresponding component, while and denote the lower and upper bounds of the analysis frequency band. By integrating the squared spectral amplitude over the specified frequency band, the energy of each component can be obtained, thereby characterizing the contributions of vibrations in different directions.
Based on this, the three-component energy is normalized to yield:
These normalized energies satisfy the following constraint (Equation 11):
In the following, and continue to denote the normalized three-component energies. A statistical analysis of the normalized results was performed, and their distribution characteristics will be discussed further below.
To comprehensively characterize the resonance and amplification effects of the site, we introduce the waveform factor proposed by , which is defined as (Equation 12):
This index incorporates both and , making it a key parameter for reflecting the site’s sensitivity to seismic motion.
3.3 Random forest regression model
In this study, the term physics-constrained refers to how physical knowledge is embedded into the data-driven framework rather than to any modification of the Random Forest algorithm itself. The physical constraint is implemented through physics-guided feature engineering–adopting the fundamental frequency derived from the H/V spectral ratio and the quarter-wavelength theory as the dominant input feature–and through empirical calibration of the – relationship using borehole data, while the Random Forest algorithm is used in its standard form.
To predict , a random forest (RF) regression model was established based on the previously constructed feature vector to characterize the nonlinear relationship between the input features and . The model is a typical ensemble learning approach that improves predictive accuracy and stability by constructing multiple decision trees and aggregating their outputs. For regression tasks, the prediction can be expressed as Equation 13:
Here, represents the prediction of the ith regression tree for the input feature vector, where i denotes the number of regression trees, and denotes the thickness predicted by the model.
In this study, the input variable is denoted as , and the output is the model-predicted equivalent thickness . To evaluate model performance, the dataset is randomly divided into a training set and a validation set, where the training set is used for model training and the validation set is used to assess generalization capability. After training, the RF model is applied to regular grid points across the study area to perform spatial prediction of , thereby obtaining its continuous spatial distribution. The specific hyperparameters used for the Random Forest model in this study are summarized in Table 1.
TABLE 1
| Hyperparameter | Value | Setting |
|---|---|---|
| Task type | Regression | Manually specified |
| Number of trees | 200 | Manually specified |
| Minimum leaf size | 5 | Manually specified |
| Number of predictors per split | ⌈m/3⌉ ≈ 3 (for m = 8 features) | Default |
| In-bag sampling fraction | 1.0 (bootstrap with replacement) | Default |
| Random seed | 42 | Manually specified |
| Feature standardization | Z-score (mean = 0, std = 1) | Applied before training |
Hyperparameters of the Random Forest regression model.
3.4 Spatial interpolation methods and data-driven methods
To evaluate the predictive performance of the proposed model, several commonly used spatial interpolation and data-driven methods are selected for comparative analysis, including Support Vector Regression (SVR), Inverse Distance Weighting (IDW), Natural Neighbor, K-Nearest Neighbors (KNN), and Radial Basis Function (RBF).
SVR is a supervised learning method based on statistical learning theory. By introducing kernel functions, the input features are mapped into a high-dimensional feature space, enabling the modeling of nonlinear relationships. While controlling model complexity, the generalization ability is improved based on the principle of structural risk minimization. The regression function can be expressed as Equation 14: denotes the training samples, is a kernel function (such as a radial basis function), and are Lagrangian multipliers, and is the bias term.
IDW is a classical deterministic spatial interpolation method based on the assumption that observations closer to the estimation point have a greater influence than those farther away. The method assigns weights using a power function of distance to perform a locally weighted average. Its expression can be written as Equation 15:where denotes the observed value at the -th sample point, represents the predicted value at the target location , is the distance between the estimation point and the sample point , is the corresponding weight, and is the distance decay exponent.
Natural Neighbor is a spatial interpolation method based on Voronoi diagrams. The method introduces the estimation point into the Voronoi diagram constructed from the observation points, and determines the weights according to the area overlap between the newly generated polygon and the neighboring polygons of the observation points, thereby achieving locally weighted interpolation. Its expression can be written as Equation 16: is a weighting coefficient determined by the proportion of the Voronoi polygon’s area, such that .
KNN is an instance-based learning method that calculates the distance between the sample to be predicted and the training samples, selects the nearest neighbors in the feature space, and uses the observations of these samples to make a prediction. The formula is as follows (Equation 17): denotes the observation of the -th nearest neighbor, and denotes the number of neighbors.
RBF is a spatial interpolation method based on radial basis functions. It constructs radially symmetric functions centered at observation points to fit and interpolate spatial variables. The predicted value is expressed as a linear combination of multiple radial basis functions, and can be written as Equation 18: represents the weight coefficient, denotes a radial basis function (such as a Gaussian function or a multiquadratic function), and denotes the Euclidean distance.
3.5 Model evaluation
To ensure the reliability and robustness of model evaluation, 5-fold cross-validation is adopted to assess model performance. In this scheme, the dataset is randomly partitioned into five non-overlapping folds of approximately equal size; the model is trained on four folds and evaluated on the remaining one, and the procedure is repeated five times so that every sample is used for validation exactly once. The random seed for the partition is fixed at 42, and the predictions from all five folds are aggregated to compute the cross-validation metrics. This scheme provides a balanced compromise between training-set size and validation coverage, and is widely adopted for machine-learning models with sample sizes around one hundred.
The coefficient of determination measures the model’s ability to explain the variance of the observed data. To further distinguish between fitting capability and generalization performance, the fitting coefficient of determination and prediction coefficient of determination are defined, representing model performance on the training dataset and validation dataset, respectively. Their expressions are given as Equation 19: represents the observed equivalent thickness of the -th sample, is the average of the observed equivalent thicknesses, and is the corresponding model prediction. The root mean square error indicates the overall magnitude of the prediction error (Equation 20):
The mean absolute error represents the average magnitude of the prediction error (Equation 21):where denotes the absolute value. By comparing the performance of RF, SVR, IDW, Natural Neighbor, KNN, RBF, and Kriging across the above evaluation metrics, the differences in accuracy and stability among the methods for spatial prediction of equivalent thickness can be systematically assessed.
Under the 5-fold cross-validation framework, the prediction error of the model can be expressed as Equation 22:Here, denotes the observed equivalent thickness of the ith sample, denotes the equivalent thickness predicted by the model when the ith sample was assigned to the validation fold, and denotes the total number of samples.
4 Results
This study conducted a statistical analysis of the energy-normalized results (see Supplementary Figure A1). The analysis validated the validity of the normalization process and the reliability of the data. Although there were some variations in energy distribution among different observation points, the overall trend was stable, with no significant outliers observed. A total of 101 valid observation points were included in the analysis.
4.1 H/V spectral ratio results
The H/V spectral ratio curves obtained from different observation points exhibit distinct peak characteristics, and the identified fundamental frequencies () from the main peaks show significant spatial variability across the study area, indicating pronounced heterogeneity in the shallow subsurface structure (see Supplementary Table A1 for the f0 values at all 101 observation points) Figure 4 presents the H/V spectral ratio curves for two representative observation points. The horizontal axis denotes frequency, while the vertical axis represents the H/V spectral ratio. Clear peaks are observed at both locations, with peak frequencies of approximately 3.60 Hz and 6.66 Hz, respectively. The peak frequency corresponds to , whereas the peak amplitude reflects the site amplification effect at that frequency.
FIGURE 4
The recorded microtremor signals exhibited distinct characteristics of ambient vibrations across all three components. Time-series analysis showed that the amplitudes of the preprocessed signals remained generally stable, with no significant transient disturbances. Spectral analysis indicated that most of the microtremor energy was concentrated within specific frequency ranges associated with the dynamic response of shallow subsurface media.
4.2 Spatial distribution characteristics of site parameters
Based on H/V spectral ratio results extracted from all observation points, spatial interpolation analysis of dynamic site parameters was carried out, and their spatial distributions are shown in Figure 5.
FIGURE 5
The results show that ranges from 2.29 to 6.71 Hz, generally exhibiting low-to mid-frequency characteristics. Low-frequency zones (approximately 2.3–3 Hz) are mainly distributed in the central and southeastern parts of the study area, whereas high-frequency zones (approximately 5–6.5 Hz) are located in the northeastern and southwestern local areas. ranges from 8.44 to 32.37 m, showing a pronounced spatial zonation pattern. The central and southeastern parts are characterized by relatively high values, with local thickness exceeding 25–30 m, while the northern and western parts are dominated by thinner layers, with thickness generally ranging from 10 to 18 m. ranges from 1.80 to 3.78, and its spatial distribution is relatively smooth, although local high-value anomalies occur in the southern and eastern parts of the study area. ranges from 0.59 to 3.55, and its spatial distribution is generally consistent with that of the amplification factor, with high values mainly concentrated in the southern and southeastern areas.
A comparison between the spatial distributions of and shows a good spatial correspondence. Low-frequency areas generally correspond to greater depths, whereas high-frequency areas correspond to thinner layers. This spatial pattern of “low frequency–thick layer and high frequency–thin layer” is generally continuous and consistent, which directly confirms the fundamental physical relationship described by the quarter-wavelength theory.
4.3 Establishment of empirical relationships and drilling verification
Based on the established quantitative relationship between and in the study area, an empirical model relating to was constructed by integrating borehole data with three-component ambient microtremor observations at adjacent locations. A total of 11 boreholes were collected in the study area, and their spatial distribution is shown in Figure 5d. For ease of comparison, the boreholes were sequentially labeled from west to east as ZK01–ZK11.
Through regression analysis, a power-law relationship between and in the study area is obtained as follows:
The empirical relationship indicates a pronounced inverse power-law relationship between and : as decreases, increases; conversely, as increases, decreases, indicating thinner sediment layers or a shallower bedrock depth.
As shown in Figure 6a, the observed data points are generally distributed along the fitted curve, indicating a clear negative correlation between and . The fitting result yields a coefficient of determination of 0.96, suggesting that the empirical model can effectively characterize the quantitative relationship between the two variables. It is worth noting that although the number of boreholes used for calibration is relatively limited (11 in total), their spatial distribution covers the main variations in geological conditions, and the established empirical relationship still exhibits high fitting accuracy and reliability.
FIGURE 6
The comparison results among different empirical models are shown in Figure 6b. It can be observed that the empirical relationship obtained in this study generally falls within the range of existing models and shows good consistency in the low-frequency range. This indicates that the proposed model has good applicability in the study area, while also reflecting the influence of regional geological differences on the parameters of empirical relationships.
To further validate the validity of the empirical relationships established, the results of this study were compared with empirical models reported in the existing literature. The mathematical expressions, scope of application, and relevant characteristic parameters of each empirical model are summarized in Table 2.
TABLE 2
| References | Empirical relationship | Study area | Frequency (Hz) | Remarks |
|---|---|---|---|---|
| This study | Right River District, Baise, China; 106 sites, 6 boreholes | 2.29–6.71 | Sedimentary layer thickness ranges from 9 to 34 m | |
| Bajo Segura Basin, Spain; 23 stations | 1–10 | Sedimentary thickness less than 100 m; shear wave velocity <250 m/s | ||
| Nanning City, Guangxi, China; 561 points compared boreholes | 1–10 | Most thicknesses between 10 and 40 m | ||
| Sanhe–Pinggu area, China; 3 boreholes and 4 arrays | 0.2–10 | Shallow shear wave velocity <180 m/s | ||
| Pearl River Delta, China; 52 boreholes | 1–10 | Sedimentary thickness 7.9–39.6 m | ||
| Po Plain region of Italy | 0.2–1 | The thickness of the sedimentary layer is less than 500 m | ||
| Lower Rhine Embayment, Germany; 34 | 0.14–4.5 | Cover thickness 30–1,600 m; shear-wave velocity varies significantly | ||
| Data from 20 boreholes and sites in Harbin, China | 1.23–4.89 | Shear wave velocity greater than 500 m/s The thickness of the overburden layer ranges from 41 to 84.5 m |
Comparison of empirical relationships between site dominant frequency and sediment thickness and their applicable conditions.
To further verify the applicability of the empirical relationship in the study area, the values obtained from microtremor observations near the boreholes were substituted into Equation 23 to calculate , and the results were compared with the soil–rock interface depths revealed by the boreholes. The calculated values and relative errors at each borehole are listed in Table 3, and the comparison profiles are shown in Figure 6c.
TABLE 3
| Drill hole | Drilling depth(m) | Observation point | Frequency (Hz) | Estimated equivalent thickness(m) | Relative error (%) |
|---|---|---|---|---|---|
| ZK01 | 15.66 | 94 | 4.32 | 14.59 | 6.83% |
| ZK02 | 24.01 | 50 | 2.77 | 25.53 | 6.33% |
| ZK03 | 21.98 | 84 | 3.22 | 21.15 | 3.78% |
| ZK04 | 14.33 | 68 | 4.32 | 14.59 | 1.81% |
| ZK05 | 21.06 | 36 | 3.08 | 22.34 | 6.08% |
| ZK06 | 19.10 | 21 | 3.51 | 19.93 | 4.35% |
| ZK07 | 15.80 | 01 | 4.02 | 16.89 | 6.90% |
| ZK08 | 20.80 | 20 | 3.43 | 20.44 | 1.73% |
| ZK09 | 17.50 | 05 | 3.73 | 18.53 | 5.89% |
| ZK10 | 11.80 | 03 | 5.11 | 12.56 | 6.44% |
| ZK11 | 24.30 | 06 | 2.81 | 25.36 | 4.36% |
Comparison between calculated equivalent thickness and borehole measurements.
The comparison results (Figure 6c; Table 3) show that the calculated is generally in good agreement with the borehole-derived depths . The relative errors at the borehole locations range from 1.73% to 6.90%, with an average error of approximately 5.32%. These results indicate that, given a properly designed methodology, a limited but representative set of borehole data can effectively calibrate a physically based empirical relationship applicable to the local area, thereby providing a reliable physical basis for subsequent large-scale spatial prediction.
4.4 Model performance and interpretation
4.4.1 Comparison of different methods
To evaluate the suitability of different spatial forecasting methods for estimating , this paper selects six methods—RF, SVR, IDW, Natural Neighbor, KNN, and RBF—for comparative analysis. ,RMSE, and MAE were used as evaluation metrics; the results are shown in Table 4.
TABLE 4
| Method | Input features | RMSE (m) | MAE (m) | |
|---|---|---|---|---|
| Full | 0.857 | 1.81 | 1.02 | |
| RBF | 0.778 | 2.25 | 1.82 | |
| SVR | 0.739 | 2.44 | 1.65 | |
| KNN | 0.726 | 2.50 | 1.90 | |
| IDW | 0.151 | 4.40 | 3.52 | |
| Natural Neighbor | 0.113 | 4.50 | 3.41 |
Performance comparison of different prediction methods.
Table 4 presents the predictive performance of different methods. In terms of the core evaluation metrics, the RF model performs best among all methods, with = 0.857, RMSE = 1.81 m, and MAE = 1.02 m, indicating superior fitting capability and predictive accuracy. In comparison, nonlinear models based on the full feature set (RBF, SVR, and KNN) show lower performance, with ranging from 0.726 to 0.778, RMSE from 2.25 to 2.50 m, and MAE from 1.65 to 1.90 m, suggesting that RF is more effective in capturing complex nonlinear relationships between multi-source seismic parameters and spatial features. Furthermore, the RF model explains more than 85% of the spatial variability in equivalent thickness and achieves higher predictive accuracy than the other models, providing a reliable basis for subsequent uncertainty analysis.
For interpolation methods that rely solely on spatial coordinates, such as IDW and Natural Neighbor, the predictive performance is the poorest, with values all below 0.16, RMSE values exceeding 4.40 m, and MAE values exceeding 3.41 m. These results indicate that relying only on spatial information is insufficient to achieve high-accuracy prediction. In contrast, incorporating signal features such as , ,energy-related features, and provides key information reflecting subsurface properties and seismic wave propagation processes, thereby significantly improving model predictive performance.
Further examination of the scatter plots of predicted versus observed values shown in Figure 7 reveals that the predictions of the RF model are highly consistent with the 1:1 reference line, indicating excellent fitting performance. The SVR model reflects the overall data trend, but exhibits noticeable dispersion in local regions. In contrast, the scatter distributions of IDW, Natural Neighbor, KNN, and RBF deviate markedly from the reference line, with predicted values concentrated within a relatively narrow range, exhibiting a pronounced smoothing effect and failing to capture the actual spatial heterogeneity. In addition, the residual distributions of the six methods are presented in Supplementary Figure A2, which provide a visual assessment of the error distribution characteristics and stability of each model.
FIGURE 7
After obtaining the predicted values at discrete observation points, natural neighbor interpolation is applied to achieve a continuous spatial representation. Based on the Voronoi tessellation principle, this method preserves the local characteristics of the original data without introducing noticeable oscillations, making it suitable for spatial reconstruction under irregularly distributed sampling conditions.
During the interpolation process, a regular grid (approximately 50 × 50 cells) is constructed over the study area to map discrete data onto a continuous spatial field. To avoid unrealistic extrapolation beyond the study area boundary, a boundary constraint is further introduced. By extracting the outer boundary of the observation points and constructing a polygonal domain, only the interpolated results within the polygon are retained, while the external region is masked. This ensures the geological plausibility of the spatial distribution results (see Figure 8a).
FIGURE 8
Based on the spatial consistency validation, the model prediction accuracy was further quantitatively evaluated using statistical metrics and residual analysis. As shown in Figure 8b, a good linear relationship is observed between the observed and predicted, with a high coefficient of determination ( = 0.952), indicating that the model has strong explanatory capability for . Figure 8c presents the residual distribution between observed and predicted values. It can be seen that the residuals are generally randomly distributed around zero, without obvious trends or clustering patterns, indicating that no systematic bias exists across different thickness ranges.
Combined with the quantitative evaluation metrics (RMSE = 1.22 m, MAE = 0.61 m), the overall prediction error of the model falls within an acceptable range for engineering applications. In addition, the uniform distribution of residuals indicates that the model has strong adaptability under different geological unit conditions. It should be noted that this error analysis is based on the predictions of the RF model for discrete samples, while the spatial continuous distribution is obtained through subsequent interpolation and gridding. Therefore, the high goodness-of-fit and residual randomness reflected in Figures 8b, 7c jointly indicate that the model maintains high accuracy in both discrete prediction and spatial representation, exhibiting high precision, low bias, and good generalization ability.
It should be noted that the range of values shown in Figure 8a (approximately 10–26 m) is smaller than that of (approximately 8–32 m). This is mainly due to the smoothing effect of the RF model and the influence of the spatial interpolation process, whereby the predicted results represent a spatially averaged response, leading to a certain reduction of extreme values. Therefore, the results reflect a relatively smooth and stable spatial distribution rather than the full variability of individual observations.
4.4.2 Feature ablation analysis
To quantitatively evaluate the contribution of each input feature to the model predictive performance, a systematic feature ablation experiment was designed and conducted in this study. By sequentially removing key features (or feature groups) from the model, the resulting model performance was compared with that of the Full model including all features. The detailed results are shown in Table 5; Figure 9.
TABLE 5
| Method | Input features | RMSE (m) | MAE (m) | |
|---|---|---|---|---|
| Full | 0.857 | 1.81 | 1.02 | |
| M 1 | Without | 0.574 | 3.12 | 2.27 |
| M 2 | Without | 0.877 | 1.68 | 0.96 |
| M 3 | Without | 0.825 | 2.00 | 1.19 |
| M 4 | Excluding three-component energy features () | 0.874 | 1.70 | 0.93 |
| M 5 | Excluding spatial coordinates () | 0.836 | 1.94 | 1.15 |
Feature ablation results of the RF model.
FIGURE 9
The Full model achieved good predictive performance when all features were included. When (M1) was removed, the model performance decreased significantly. Compared with the Full model, decreased by approximately 0.283, while RMSE and MAE increased by 72.4% and 122.5%, respectively. This indicates that is the key controlling feature for achieving accurate prediction.
After removing (M2) or the three-component energy feature group (M4), all performance metrics of the model were better than those of the Full model. The of M2 increased to 0.877, and RMSE and MAE decreased to 1.68 and 0.96 m, respectively; M4 achieved better (0.874) and MAE (0.93 m). This indicates that these features may introduce additional noise in the Full model or exhibit a certain degree of redundancy and multicollinearity with other features.
Finally, after removing (M3) or spatial coordinates (M5), the model performance showed a slight but consistent decrease. For example, the of M3 decreased to 0.825, and that of M5 decreased to 0.836, with corresponding increases in errors. This indicates that and spatial coordinates provide useful supplementary information for the prediction task, but their contributions are far less significant than that of the core feature .
Figure 9 shows the residual scatter plots of the RF model under different feature ablation scenarios. It can be seen that the residuals of the M1 model (without ) are the most dispersed and clearly deviate from zero, indicating a significant increase in prediction error and a substantial reduction in the model’s ability to characterize the target variable. In contrast, the residuals of the other models (M2–M5) are mainly concentrated around zero and exhibit relatively compact patterns, indicating lower overall error levels. Among them, M2 and M4 show the most concentrated residual patterns, further indicating that the removal of these features does not weaken model performance. Overall, the differences in residual patterns clearly demonstrate that the absence of leads to a significant amplification of prediction errors and is a key controlling factor affecting prediction accuracy, whereas the effects of other features are relatively limited. The residual scatter plots of all ablation models (M1–M5) are provided in Supplementary Figure A3.
4.4.3 Feature importance and partial dependence analysis
To further reveal the prediction mechanism of the model, feature importance analysis and partial dependence analysis were conducted (see Figure 10). The results show that has the most significant contribution to the prediction of , followed by and , whose contributions are relatively weaker (Figure 10a). This indicates that plays a dominant role in characterizing the subsurface structure, which is consistent with its high sensitivity to the structure.
FIGURE 10
Partial dependence analysis further reveals the nonlinear relationships between each feature and (Figures 10b–d). As increases, decreases significantly (Figure 10b), exhibiting a clear negative correlation, which is consistent with the theoretical understanding that higher frequencies correspond to shallower subsurface structures. In contrast, the influence of is more complex: its response curve shows noticeable fluctuations with an overall limited contribution (Figure 10c), suggesting that its effect may be influenced by interactions with other variables or local nonlinear effects.
shows a positive correlation with (Figure 9d), with a pronounced variation in the range of 1.0–2.0, followed by a gradual leveling off, indicating a marginal diminishing effect. To emphasize the contribution of physically meaningful parameters, spatial coordinates (latitude and longitude) were excluded during feature interpretation to prevent them from dominating the model explanation due to spatial autocorrelation, thereby improving the physical interpretability and reliability of the results. This result further confirms that is the decisive factor in predicting , provides supplementary explanatory information, and the independent influence of is relatively weak.
5 Discussion
This study proposes and validates a machine learning framework that integrates physical mechanisms with data-driven approaches to achieve high-resolution spatial prediction of shallow in complex urban environments. The proposed framework not only significantly improves prediction accuracy but also achieves a good balance among physical consistency, model interpretability, and engineering applicability. The following sections systematically discuss the sources of model performance, interpretability value, methodological innovations, uncertainty analysis, as well as engineering applications and future perspectives.
5.1 Synergistic mechanism underlying superior model performance
The RF model in this study demonstrates superior predictive performance compared with traditional spatial interpolation methods (IDW, Natural Neighbor) and other data-driven methods (SVR, RBF, KNN), achieving an of 0.857. This superiority is rooted in the deep synergy between physics-guided feature engineering and the intrinsic advantages of the algorithm. The model is able to explain more than 85% of the spatial variability in the data, with prediction accuracy significantly exceeding that of conventional methods. Meanwhile, the approximately 15% of unexplained variability may be attributed to microtremor observation noise, insufficiently characterized local geological anomalies, and the inherent simplifications of the concept when applied to complex stratigraphic conditions.
First, the input feature system has a solid physical foundation. The fundamental frequenc is directly derived from the quarter-wavelength resonance theory and serves as a key physical parameter linking wavefield observations to subsurface structures. The parameter integratesand , while the three-component energy features characterize the anisotropy of the wavefield, and spatial coordinates implicitly reflect variations in the geological background. This feature construction ensures that the model learns relationships driven by geophysical mechanisms rather than spurious statistical correlations.
Second, the random forest algorithm is well suited for capturing high-dimensional nonlinearity and interactions among features. Urban subsurface structures exhibit pronounced spatial heterogeneity and parameter coupling. By integrating multiple decision trees, RF can effectively model such complex relationships. Moreover, its built-in random subspace sampling mechanism enhances model robustness under limited sample conditions, thereby mitigating the risk of overfitting.
Feature ablation experiments (Table 5) provide empirical support for this finding: model performance is highly sensitive to the removal of the core physical feature , whereas it remains stable or even improves when redundant features, such as or three-component energy features, are removed. This indicates that the model can automatically identify and rely on key physical variables while effectively suppressing the influence of noise or collinear features, thereby achieving high prediction accuracy under small-sample conditions.
5.2 Interpretability analysis: from “black-box” prediction to “glass-box” mechanism
This study systematically conducts interpretability analysis to transform the machine learning model from a traditional “black-box” prediction tool into a “glass-box” model with clear decision logic and consistency with physical principles. The analysis not only verifies the consistency between the model’s decision-making process and geophysical theory, but also reveals the complex interactions among features and their dominant controlling mechanisms, thereby establishing a reliable interpretive bridge between data-driven approaches and physical theory.
The feature importance analysis (Figure 10A) indicates that is the dominant controlling factor for predicting , with an importance significantly higher than that of other input variables. This result is not incidental, but rather provides strong data-driven support for the physical basis of the H/V spectral ratio method. It quantitatively validates the core concepts of the quarter-wavelength theory expressed in Equations 4, 6: within a given range of , serves as the most sensitive and direct physical proxy controlling . The model’s strong reliance on this feature fundamentally ensures the physical plausibility and theoretical consistency of the predictions.
The partial dependence analysis (Figures 10b–d) further reveals complex and nonlinear functional relationships between key features and the target variable, offering mechanistic insights beyond simple statistical correlations. The most prominent pattern is the stable monotonic decreasing relationship between and (Figure 10b), which is highly consistent with the physical law described by the quarter-wavelength theory—namely, higher frequencies correspond to smaller equivalent thicknesses. This indicates that the model is not performing unconstrained empirical fitting, but rather successfully reconstructs and explicitly represents this key geophysical mechanism within a data-driven framework.
This study, through systematic interpretability analysis, transforms the machine learning model from a purely “black-box” prediction tool into a “glass-box” model with transparent decision logic and clear physical consistency. This analysis not only verifies the consistency between the model’s decision-making process and geophysical principles, but also reveals the complex interactions and dominant mechanisms among features, ultimately establishing a robust interpretability bridge between data-driven approaches and physical theory.
The role of parameter further reflects the model’s capability to learn more complex site response patterns. The partial dependence results show a positive correlation between and (Figure 10d), indicating that the predicted equivalent thickness increases with increasing values. From a physical perspective, an increase in may arise from two mechanisms: either stronger site amplification (larger ) or lower resonance frequency (smaller ), both of which are typically associated with thicker or softer sedimentary layers. Therefore, the model’s positive dependence on suggests that it successfully captures a coupled and composite site response mechanism, in which thick and soft sediments are often characterized by both more pronounced resonance amplification and lower fundamental frequencies. This highlights the ability of the machine learning model to identify and exploit complex feature interactions without relying on explicitly defined physical equations.
In contrast, shows a relatively low contribution, and its partial dependence curve exhibits stronger fluctuations and nonlinearity (Figure 10c), further indicating a more complex and less stable relationship with . This observation is consistent with the feature ablation experiment (M2, Table 5), where removing led to a slight improvement in overall model performance ( = 0.877). This suggests that, under the specific site conditions and data quality of the study area, the information carried by that is independent of and effective for thickness prediction is relatively limited. On the one hand, its signal may exhibit strong collinearity with the core feature ; on the other hand, it may contain substantial local noise unrelated to thickness prediction (e.g., shallow surface heterogeneity, sensor coupling differences).
This insight should not be regarded as a limitation of the model, but rather as a finding of practical significance. It indicates that, in future observation design and feature engineering, priority should be given to ensuring the accuracy and stability of extraction, while may be more appropriately used as a supplementary reference parameter or for other site evaluation purposes.
5.3 Methodological innovation: from fixed-parameter empirical models to adaptive physics-informed machine learning
The core contribution of this study lies in establishing a modeling paradigm that deeply integrates physical mechanisms with data-driven approaches. Compared with traditional methods, this paradigm achieves three key methodological transformations, thereby significantly enhancing the model’s representational capacity and generalization ability. This improvement ultimately leads to the superior predictive performance (
= 0.857):
From global equivalent assumptions to spatially variable parameter modeling: enabling adaptive correction of lateral heterogeneity. Traditional empirical formulations based on the quarter-wavelength theory (e.g., Equation 6) or regionally calibrated approaches essentially adopt a spatially constant to establish the relationship between frequency and thickness, thereby neglecting the lateral heterogeneity that commonly exists in urban subsurface structures. The RF model constructed in this study, by introducing spatial coordinates and multi-source physical features as joint inputs, implicitly learns and characterizes the continuous spatial variation of . In other words, when performing thickness prediction, the model does not rely on a global average parameter, but instead adaptively selects an optimal equivalent shear-wave velocity based on local feature combinations. The feature ablation experiment (M5, Table 5) shows that when spatial coordinates are removed, the model performance decreases from = 0.857 to 0.836, thereby inversely verifying the critical role of spatial information in capturing local velocity variations and improving prediction accuracy. This mechanism constitutes the fundamental reason why the model outperforms traditional empirical formulations (Figure 6A).
From discrete point-based calibration to continuous spatial field learning: integrating sparse constraints with dense spatial observations. Traditional methods rely heavily on dense and costly borehole data for calibration and validation, which essentially represent an interpolation process of “using points to approximate a surface.” In this study, within a machine learning framework, the sparse calibration information provided by 11 boreholes is integrated with dense spatial observational data obtained from 101 microtremor measurement points. The model not only accurately learns the – relationship at borehole locations (Table 3, with an average error of 5.32%), but more importantly, effectively captures the spatial co-variation characteristics and distribution patterns of site dynamic parameters (, , ) across the study area through the dense observation network (Figure 8). Ultimately, the model generates a high-resolution continuous prediction field covering the entire study area with physical consistency (Figure 8a), rather than a simple interpolation result based on discrete points, thereby significantly improving exploration efficiency and prediction accuracy.
From result-oriented prediction to interpretable modeling: constructing a physically consistent and trustworthy decision logic. Traditional machine learning models are often regarded as black-box approaches, whose predictions are difficult to gain full trust from domain experts. In this study, through systematic feature importance analysis and partial dependence analysis (Figure 10), the decision-making process of the model is expressed in an interpretable manner. The study not only achieves high-accuracy predictions of , but further reveals the underlying decision basis of the model: it identifies as the dominant controlling factor (Figure 10a), follows the negative physical relationship between and (Figure 10b), and simultaneously reflects the composite site effects represented by parameter (Figure 10d).
This extension from “prediction results” to “mechanistic interpretation” elevates the model from a simple fitting tool to an analytical framework capable of revealing underlying physical mechanisms and providing interpretable decision logic, thereby significantly enhancing the credibility and practical value of the results in both scientific understanding and engineering applications.
5.4 Convergence and uncertainty of prediction results and their implications for the methodological paradigm
The physics-constrained, data-driven modeling framework constructed in this study produces not only a high-accuracy spatial prediction field, but also prediction results whose value range and spatial distribution characteristics contain important physical and geological information. The values predicted by the RF model converge within the range of 10–26 m (Figure 8a), which is more concentrated than the theoretical range (8–32 m) derived from traditional physical inversion based on broad velocity assumptions (e.g., = 150–500 m/s).
It should be emphasized that this convergence does not reflect a limitation of the model’s representational capacity, but rather a direct manifestation of the data-driven model adaptively learning regional geological conditions under physical constraints. Specifically, by integrating a limited but representative set of borehole calibration data (11 boreholes) with dense microtremor observations (101 measurement points), the model effectively filters out “thickness–velocity–frequency” combinations that, although consistent with ideal physical relationships, have a low geological probability in the sedimentary environment of the study area (the Youjiang River alluvial setting in Baise, China).
Therefore, the convergence of the predicted range signifies a transition from a theoretically feasible solution space to a geologically plausible solution space, which constitutes an important foundation for its engineering applicability.
However, a certain degree of uncertainty still exists within the converged solution space. The systematic spatial overestimation and underestimation in the predictions (as revealed by the borehole validation results in Figure 6c and the residual distribution in Figure 8b) indicate that the uncertainty mainly arises from the following aspects. First, local noise in the feature extraction process may lead to deviations in key parameters such as , particularly in areas with strong interference. Second, the Heq parameter represents a simplified expression of complex stratigraphic structures; in the presence of significant vertical heterogeneity (e.g., gravel lenses), a single Heq is insufficient to fully characterize the actual site response. Third, the limited coverage of extreme geological conditions in the training samples leads to a certain degree of conservatism when the model is extrapolated spatially. These uncertainties do not represent deficiencies of the method; rather, they clearly delineate the current applicability boundaries of the model and point to directions for future improvement.
5.5 Engineering applications and future perspectives
The core value of the RF model lies in providing a novel technical paradigm for rapid, high-resolution, and low-cost investigation of shallow urban subsurface structures. Under the constraint of limited borehole data, this method integrates a dense ambient microtremor observation network with a random forest model, successfully achieving high-accuracy spatial prediction of . The resulting high-resolution distribution map can directly support engineering site classification, the determination of input parameters for ground motion simulation, and refined assessment of urban seismic hazard, demonstrating clear practical engineering value. This approach is particularly suitable for built-up areas where conventional exploration methods are difficult to implement, as well as for emergency scenarios requiring rapid surveys, exhibiting significant advantages in terms of efficiency, cost-effectiveness, and environmental friendliness.
Building upon the current research outcomes, the further development of this method can be advanced from multiple perspectives. Future studies may focus on incorporating stronger physical constraints, such as integrating or dispersion curve data obtained from active-source surface wave exploration, to directly constrain the shallow shear-wave velocity structure, thereby effectively reducing the inherent “thickness–velocity” non-uniqueness in current inversion.
At the same time, it is crucial to promote the model from deterministic prediction to probabilistic prediction, in order to quantify the unexplained variance (approximately 15%) identified in this study. In the future, uncertainty in predictions can be quantified by incorporating Bootstrap or Bayesian frameworks, allowing the output of confidence intervals and thus providing more reliable inputs for risk-based probabilistic seismic hazard analysis. In addition, improving the generalization ability of the model is crucial for transforming it into a universally applicable tool. Exploring transfer learning techniques, by leveraging multi-regional datasets to construct a pre-trained model, would enable rapid adaptation to new areas with only a small number of local samples, thereby significantly enhancing its cross-regional applicability.
From a longer-term perspective, advancing the inversion from the integrated parameter toward a “refined velocity structure”is an inevitable direction. Future studies may explore the integration of one-dimensional physical models with physics-informed neural networks to achieve direct inversion of continuous profiles with depth from microtremor data, ultimately providing a truly transparent three-dimensional physical property model of the urban subsurface.
6 Conclusion
To address the challenges of sparse samples, insufficient physical constraints, and limited model interpretability in the spatial prediction of shallow
in complex urban sites, this study proposes a physics-informed data-driven machine learning framework. Based on dense three-component ambient microtremor observations and limited borehole calibration data, the framework achieves high-accuracy and high-resolution spatial prediction of
in the study area. The main contributions and key findings of this study are as follows:
Validation of methodological effectiveness: The constructed RF model significantly outperforms traditional spatial interpolation methods (IDW, Natural Neighbor) and other comparative data-driven methods (SVR, RBF, KNN) in predictive performance, achieving an of 0.857. This indicates that integrating physically derived features from the H/V spectral ratio method (, , , three-component energy, and spatial coordinates) with the nonlinear modeling capability of machine learning can effectively characterize the spatial heterogeneity of complex subsurface structures.
Physical interpretability of model decision-making: Through systematic feature importance analysis and partial dependence analysis, the “black box” of the model is effectively opened. The results indicate that is the most critical feature for predicting and exhibits a clear negative correlation with , strongly supporting the fundamental physical principle of the quarter-wavelength theory. This analysis not only enhances the reliability of the model results but also provides data-driven insights into the mechanisms of site dynamic response.
Synergistic advantage of physics–data integration: Feature ablation experiments show that model performance is highly sensitive to the core physical feature , while remaining relatively robust when certain features (e.g., and ) are removed. This reveals the intrinsic advantage of the proposed framework: physics-based mechanisms guide feature engineering to ensure a reasonable learning direction, while the data-driven capability of RF enables adaptive capture of complex nonlinear interactions and spatial variability among features, thereby improving prediction accuracy. Meanwhile, the predicted results exhibit clear convergence and a significantly reduced range compared with traditional theoretical estimates, indicating a transition from the theoretically feasible solution space to a geologically plausible solution space. In essence, the framework achieves an implicit and adaptive modeling of the spatial variability function of in the study area.
This work not only represents a successful application of the proposed methodology, but also further demonstrates that the deep integration of physical mechanisms and data-driven approaches is an effective pathway for advancing near-surface geophysical exploration and interpretation toward higher accuracy and stronger interpretability.
Statements
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Author contributions
YO: Conceptualization, Data curation, Formal Analysis, Methodology, Writing – original draft.
Funding
The author(s) declared that financial support was not received for this work and/or its publication.
Conflict of interest
The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2026.1841221/full#supplementary-material
References
1
AhmedK. A.KhanS.NisarU. B.MughalM. R.SultanM. (2021). Machine seismic: an automatic approach for the identification of subsurface structural models. Soft Comput.25 (13), 8169–8176. 10.1007/s00500-021-05740-2
2
ChenQ.GaoG.-Y.YangJ. (2011). Dynamic response of deep soft soil deposits under multidirectional earthquake loading. Eng. Geol.121 (1-2), 55–65. 10.1016/j.enggeo.2011.04.013
3
CultreraG.De RubeisV.TheodoulidisN.CadetH.BardP.-Y. (2014). Statistical correlation of earthquake and ambient noise spectral ratios. Bull. Earthq. Eng.12 (4), 1493–1514. 10.1007/s10518-013-9576-7
4
DelgadoJ.López CasadoC.GinerJ.EstévezA.CuencaA.MolinaS. (2000). Microtremors as a geophysical exploration tool: applications and limitations. Pure Applied Geophysics157 (9), 1445–1462. 10.1007/PL00001128
5
DuanY.ShenY.CanbulatI.LuoX.SiG. (2021). Classification of clustered microseismic events in a coal mine using machine learning. J. Rock Mech. Geotechnical Eng.13 (6), 1256–1273. 10.1016/j.jrmge.2021.09.002
6
GuoJ.ZhengY.LiuZ.WangX.ZhangJ.ZhangX. (2025). Pattern-based multiple-point geostatistics for 3D automatic geological modeling of borehole data. Nat. Resour. Res.34 (1), 149–169. 10.1007/s11053-024-10405-6
7
HaghshenasE.BardP. Y.TheodulidisN.TeamS. W. (2008). Empirical evaluation of microtremor H/V spectral ratio. Bull. Earthq. Eng.6 (1), 75–108. 10.1007/s10518-007-9058-x
8
HarsukoM. R. C.ZulfakrizaZ.NugrahaA. D.SarjanA. F. N.WidiyantoroS.RosaliaS.et al (2020). Investigation of hilbert–huang transform and fourier transform for horizontal-to-vertical spectral ratio analysis: understanding the shallow structure in Mataram city, Lombok, Indonesia. Front. Earth Sci.8, 334. 10.3389/feart.2020.00334
9
HuynhN. N. T.MartinR.OberlinT.PlazollesB. (2023). Near-surface seismic arrival time picking with transfer and semi-supervised learning. Surv. Geophys.44 (6), 1837–1861. 10.1007/s10712-023-09783-y
10
Ibs-von SehtM.WohlenbergJ. (1999). Microtremor measurements used to map thickness of soft sediments. Bull. Seismol. Soc. Am.89 (1), 250–259. 10.1785/bssa0890010250
11
KalininaA. V.AmmosovS. M.TatevossianR. E.BykovaV. V. (2022). Applicability of the H/V method in the seismic microzoning problem. Seism. Instrum.58 (1), S79–S88. 10.3103/S0747923922070052
12
LehmannF.GattiF.BertinM.ClouteauD. (2022). Machine learning opportunities to conduct high-fidelity earthquake simulations in multi-scale heterogeneous geology. Front. Earth Sci.10, 10–2022. 10.3389/feart.2022.1029160
13
LiR.ProzziJ. A.HongF. (2025). Quantification of post-rainfall moisture content in pavement unbound layers using long-term pavement performance data. Transp. Res. Rec.2680, 03611981251380276–03611981251380672. 10.1177/03611981251380276
14
LiangD.GanF.ZhangW.JiaL. (2018). The application of HVSR method in detecting sediment thickness in karst collapse area of pearl river Delta, China. Environ. Earth Sci.77 (6), 259. 10.1007/s12665-018-7439-x
15
LiangR.ZhangC.HuangC.LiB.SaydamS.CanbulatI.et al (2024). Multimodal data fusion for geo-hazard prediction in underground mining operation. Comput. and Industrial Eng.193, 110268. 10.1016/j.cie.2024.110268
16
LiuY.ShiL. (2018). Site characteristic parameters’ quick measurement based on micro-tremor's H/V spectra. Zhendong Yu Chongji/Journal Vib. Shock.37 (13), 235–242. 10.13465/j.cnki.jvs.2018.13.037
17
LiuP.-C.TsaiC.-C. (2022). Influence of local site condition on vertical-to-horizontal spectrum ratio – insight from site response analysis. J. Earthq. Eng.26 (5), 2283–2300. 10.1080/13632469.2020.1759473
18
LvX.WangG. (2024). GIS-Based mineral prospectivity mapping using machine learning methods: a case study from Duobaoshan Ore district, Northeastern China. Ore Geol. Rev.175, 106352. 10.1016/j.oregeorev.2024.106352
19
LvP.ChenW.LiH.SongW. (2024). SsL-VGMM: a semisupervised machine learning model of multisource data fusion for lithology prediction. Nat. Resour. Res.33 (5), 1993–2007. 10.1007/s11053-024-10375-9
20
MaH.LiaoB.-Y. (2024). Study on site effect in Nanning city using HVSR method. Civ. Eng. Archit.12 (2), 1165–1179. 10.13189/cea.2024.120235
21
MacauA.BenjumeaB.GabàsA.FiguerasS.VilàM. (2015). The effect of shallow Quaternary deposits on the shape of the H/V spectral ratio. Surv. Geophys.36 (1), 185–208. 10.1007/s10712-014-9305-z
22
MaoH.MaJ.WangX.JiangL.WangL.ZhangY.et al (2020). Harmonic noise suppression of vibroseis data based on adaptive dictionary learning. Geophys. Prospect. Petroleum59 (5), 725–735. 10.3969/j.issn.1000-1441.2020.05.006
23
MarzánI.MartíD.LoboA.AlcaldeJ.RuizM.Alvarez-MarrónJ.et al (2021). Joint interpretation of geophysical data: applying machine learning to the modeling of an evaporitic sequence in Villar de Cañas (Spain). Eng. Geol.288, 106126. 10.1016/j.enggeo.2021.106126
24
MascandolaC.MassaM.BaraniS.AlbarelloD.LovatiS.MartelliL.et al (2019). Mapping the seismic bedrock of the Po Plain (Italy) through ambient‐vibration monitoring. Bull. Seismol. Soc. Am.109 (1), 164–177. 10.1785/0120180193
25
MokhberiM. (2015). Vulnerability evaluation of the urban area using the H/V spectral ratio of microtremors. Int. J. Disaster Risk Reduct.13, 369–374. 10.1016/j.ijdrr.2015.06.012
26
MolnarS.SiroheyA.AssafJ.BardP.-Y.CastellaroS.CornouC.et al (2022). A review of the microtremor horizontal-to-vertical spectral ratio (MHVSR) method. J. Seismol.26 (4), 653–685. 10.1007/s10950-021-10062-9
27
NakamuraY. (1989). A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface. Q. Rep. RTRI (Railw. Tech. Res. Inst.). 30 (1), 25–33.
28
PengC.WangC.LiZ. (2025). Review of geophysical data acquisition methods for underground feature detection and future trends. Tunn. Undergr. Space Technol.163, 106731. 10.1016/j.tust.2025.106731
29
RuoHanZ.PeiFenX.SuQunL.YaNanD.ZhiWeiY.ZhiHuiW.et al (2020). Detection of the soil-rock interface based on microtremor H/V spectral ratio method: a case study of the Jinan urban area. Chin. J. Geophys. (In Chinese)63 (1), 339–350. 10.6038/cjg2020M0678
30
SignaniniP.TorreseP. (2004). Application of high resolution shear-wave seismic methods to a geotechnical problem. Bulletin of Engineering Geology and the Environment63 (4), 329–336. 10.1007/s10064-004-0252-7
31
WQ. T.WeijunW.HuadongK. (2020). H/V spectral analysis of ground vibration in the Sanhe-Pinggu area: plateau response, shallow sedimentary structure and its reflected fault activity. Chinese Journal of Geophysics (In Chinese)63 (10), 16. 10.19975/i.dqyxx.2021-003
32
WangZ.ZuoR.JingL. (2021). Fusion of geochemical and remote-sensing data for lithological mapping using random forest metric learning. Mathematical Geosciences53 (6), 1125–1145. 10.1007/s11004-020-09897-8
33
XuR.WangL. (2021). The horizontal-to-vertical spectral ratio and its applications. EURASIP Journal on Advances in Signal Processing.2021 (1), 75. 10.1186/s13634-021-00765-z
34
XuZ.ZhuangH.XiaZ.YangJ.BuX. (2023). Study on the effect of burial depth on seismic response and seismic intensity measure of underground structures. Soil Dynamics and Earthquake Engineering166, 107782. 10.1016/j.soildyn.2023.107782
35
YangX.LiS.CaoA.WangC.LiuY.BaiX.et al (2024). Deep transfer learning for P-wave arrival identification and automatic seismic source location in underground mines. International Journal of Rock Mechanics and Mining Sciences182, 105888. 10.1016/j.ijrmms.2024.105888
36
YuanF. (2026). Reconstructing the evolution of fault activity in Feidong, Anhui province by integrating multimodal deep learning and geological big data. Frontiers in Earth Science14, 14–2026. 10.3389/feart.2026.1753000
37
ZhanW.ChenY.LiuQ.LiJ.SacchiM. D.ZhuangM.et al (2024). Simultaneous prediction of petrophysical properties and formation layered thickness from acoustic logging data using a modular cascading residual neural network (MCARNN) with physical constraints. Journal of Applied Geophysics224, 105362. 10.1016/j.jappgeo.2024.105362
38
ZhangY.HeL.-L.JiaoY.-Y.PengH.-F.LiuS.-C.ZhangQ.-B. (2025). Multiscale progressive 3D geological modeling based on isochronous stratigraphy identification in urban underground space. Bulletin of Engineering Geology and the Environment84 (4), 171. 10.1007/s10064-025-04185-3
39
ZhongZ.ShenY.ZhaoM.LiL.DuX. (2022). Seismic performance evaluation of two-story and three-span subway station in different engineering sites. Journal of Earthquake Engineering26 (14), 7505–7535. 10.1080/13632469.2021.1964647
40
ZhouF.LiM.HuangC.LiangH.LiuY.ZhangJ.et al (2022). Lithology-based 3D modeling of urban geological attributes and their engineering application: a case study of Guang’an city, SW China. Frontiers in Earth Science10, 10. 10.3389/feart.2022.918285
41
ZhuX.PeiX.YangS.WangW.DongY.FangM.et al (2024). Spatial prediction of ground substrate thickness in shallow Mountain area based on machine learning model. Frontiers in Earth Science12, 1455124. 10.3389/feart.2024.1455124
42
ZouC.ZhaoL.XuM.ChenY.GengJ. (2021). Porosity prediction with uncertainty quantification from multiple seismic attributes using random forest. Journal of Geophysical Research Solid Earth126 (7), e2021JB021826. 10.1029/2021jb021826
43
ZuoR.CarranzaE. J. M. (2023). Machine learning-based mapping for mineral exploration. Mathematical Geosciences55 (7), 891–895. 10.1007/s11004-023-10097-3
Summary
Keywords
fundamental frequency, H/V, machine learning, random forest, shallow equivalent thickness, spatial prediction
Citation
Ou Y (2026) Physics-constrained random forest for spatial prediction of shallow equivalent thickness from three-component microtremor data. Front. Earth Sci. 14:1841221. doi: 10.3389/feart.2026.1841221
Received
28 March 2026
Revised
22 May 2026
Accepted
25 May 2026
Published
15 June 2026
Volume
14 - 2026
Edited by
Ashutosh Chamoli, Indian Institute of Technology Roorkee, India
Reviewed by
R. B. S. Yadav, Kurukshetra University, India
Weichen Zhan, The University of Texas at Austin, United States
Updates
Copyright
© 2026 Ou.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yangkai Ou, ouyangkai.research@gmail.com
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.