Three-Dimensional Stratigraphic Structure and Property Collaborative Modeling in Urban Engineering Construction

Abstract

In urban engineering construction, ensuring the stability and safety of subsurface geological structures is as crucial as surface planning and aesthetics. This study proposes a novel multivariate radial basis function (MRBF) interpolant for the three-dimensional (3D) modeling of engineering geological properties, constrained by the stratigraphic structural model. A key innovation is the incorporation of a well-sampled geological stratigraphical potential field (SPF) as an ancillary variable, which enhances the interpolation of geological properties in areas with sparse and uneven sampling points. The proposed MRBF method outperforms traditional interpolation techniques by showing reduced dependency on the distribution of sampling points. Furthermore, the study calculates the bearing capacity of individual pile foundations based on precise stratigraphic thicknesses, yielding more accurate results compared to conventional methods that average these values across the entire site. Additionally, the integration of 3D geological models with urban planning facilitates the development of comprehensive urban digital twins, optimizing resource management, improving decision-making processes, and contributing to the realization of smart cities through more efficient data-driven urban management strategies.

1. Introduction

In urban engineering construction, the development of a fine-grained engineering geological model is critical not only for enhancing engineering geological work and digital city frameworks but also for the quantitative evaluation of the urban environment, urban planning, subsurface space development, and rapid disaster response [1,2]. The three-dimensional (3D) visualization and characterization of geologic bodies offer a more realistic and intuitive depiction of subsurface geologic structures and properties [3]. Since Houlding (1994) introduced the concept of 3D geological modeling [4], it has been a vital tool in geological information processing and scientific research. Three-dimensional geological modeling comprises two key components, i.e., structural modeling and property modeling. Structural modeling involves constructing models by combining geometric entities, i.e., points, lines, surfaces, and bodies, with a focus on the shape of geological structures and their interrelationships. Property modeling, on the other hand, divides the 3D geological space into regular or irregular voxels, focusing on the spatial distribution of physical and chemical properties. Property models can be seamlessly integrated with geophysical fields [5,6,7], geochemical fields [8,9,10], structural analysis [11,12,13], metallogenic analysis [14,15], and reservoir characterization [16,17,18].

Spatial interpolation is a fundamental technique for property modeling, wherein the properties of sampling points are used to infer the properties of points to be estimated. Commonly used interpolation methods include kriging [19,20,21], discrete smooth interpolation (DSI) [22,23], inverse distance weighted (IDW) [24,25], and radial basis function (RBF) [26,27,28]. The performance of spatial interpolation is significantly influenced by data density and the spatial distribution of sample points [29]. Generally, well-sampled ancillary variables that provide useful information about the under-sampled primary variable should be considered. Wang et al. (2005) compared cokriging (CK) and ordinary kriging (OK) methods to estimate the quantity and distribution of soil chemicals from a limited number of available samples, showing that cokriging of half of the observations, using the same number of samples for the secondary variable (total salt), resulted in more accurate results than the OK of all the samples [30]. Guastaldi et al. (2013) demonstrated that multivariate cokriging spatial interpolation of the under-sampled primary airborne gamma-ray data, using a well-sampled geologic map as ancillary variable, exhibited a distinct correlation between geological formations and radioactivity contents [31]. Zhu et al. (2013) transformed qualitative geological constraints into quantitative control parameters during data preprocessing and applied different property interpolation schemes to different types of geo-objects, resulting in geologically reasonable property models controlled by geological structures [32]. Foehn et al. (2018) explored the regression cokriging approach to efficiently combine weather radar data with data from two heated rain gauge networks of different quality, demonstrating the added value of the radar information compared to using only ground station data [33]. Liu et al. (2024) developed a multivariate spatial interpolation method based on Yang Chizhong filtering (CoYangCZ), integrating geometry and statistics-based strategies to improve the accuracy of a variable of air pollution using ancillary meteorological datasets [34]. However, the accuracy of these methods largely depends on the sufficient number and proper distribution of sampling points. In practice, the number of geological property sampling points is often limited and unevenly distributed due to high costs. Additionally, geological property interpolation that only considers the spatial distance and direction effects between scattered points is insufficient. In 3D geological space, a stratum deposited during the same geological era will exhibit strong similarity in its physical and chemical properties, e.g., lithology and materials, even when spatial distances are maximum.

Geological property is strata-bound, and its interpolation should be conducted under the constraints of the geological structures [32]. A 3D stratigraphic potential field (SPF), using an implicit mathematical function, describes the continuous stratigraphic distribution and represents geological boundaries [35]. Using an SPF to express a set of conformable strata in 3D space is convenient for spatial analysis, statistics, and simulation. This study proposes the multivariate RBF (MRBF), incorporating the SPF as an ancillary variable into the engineering geological property interpolation algorithm to achieve strata-bound results. The property model is constructed using a hexahedral grid, with the properties of each grid cell interpolated using the MRBF algorithm. Three properties, i.e., natural water content, natural pore ratio, and water saturation, are modeled with the constraints of the 3D SPF model. The proposed method provides a more precise tool for modeling subsurface geological conditions, which is crucial for designing safe and stable foundations in urban engineering construction. Subsequently, the bearing capacity of the pile foundation for each building is designed and calculated based on detailed stratigraphic models, which minimizes geotechnical risks and enhances construction safety. Finally, the surface building inclined photogrammetric models are integrated with the subsurface 3D engineering geological models.

2. Methods and Materials

2.1. Multivariate Radial Basis Function (MRBF)

2.1.1. RBF Interpolation

RBF interpolation involves constructing a linear combination of RBFs from sampling points, where the RBF value typically depends on the distance between the point to be interpolated and the sampling points [36]. Given discrete property points ${\left\{\left({\mathit{u}}_i,f_i\right)\right\}}_{i=1}^N\in R^n\times R$, where ${\mathit{u}}_i$ denotes the position of the point i and f_i denotes the property value at the point i, the interpolation function can be expressed as follows:

$$\widehat{f}\left(\mathit{u}\right)=\sum _{i=1}^N{w}_i\phi \left(\mathit{u}-{\mathit{u}}_{\mathit{i}}\right)$$

(1)

where N is the total number of RBFs used for interpolation, $\phi \left(‖\mathit{u}-{\mathit{u}}_i‖\right)$ represents the generalized form of the RBFs used, $‖\mathit{u}-{\mathit{u}}_i‖$ denotes the Euclidean distance between the position vectors of $\mathit{u}$ and ${\mathit{u}}_i$, and $w_i$ is the corresponding weight coefficient. The commonly used RBFs are listed in Table S1.

By substituting the coordinates ${\mathit{u}}_i$ and the corresponding property values $f_i$, the equations can be expressed in matrix form as follows:

$$\mathit{\Phi }·\mathit{W}=\mathit{F}$$

(2)

$$\left(\begin{array}{cccc}{\phi }_{11}& {\phi }_{12}& \cdots & {\phi }_{1n}\\ {\phi }_{21}& {\phi }_{22}& \cdots & {\phi }_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ {\phi }_{n1}& {\phi }_{n2}& \cdots & {\phi }_{nn}\end{array}\right)\cdot \left(\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_n\end{array}\right)=\left(\begin{array}{c}{f}_1\\ f_2\\ ⋮\\ f_n\end{array}\right)$$

(3)

where ${\phi }_{ij}=\phi \left(‖{\mathit{u}}_i-{\mathit{u}}_j‖\right)$. The weight coefficients can be calculated by solving the linear system of equations. Once the coefficients are determined, the interpolation function $\widehat{f}\left(\mathit{u}\right)$ can be obtained, and the property values can be calculated by substituting the coordinates of the points to be interpolated.

2.1.2. MRBF Interpolation

When a property filled exhibits linear characteristics, these linear characteristics can also be incorporated into the RBF interpolant [37,38]. Inspired by this, when interpolating the properties of the strata, especially when the number of sampling points is limited, we can introduce an SPF model, which reflects the depositional order of the strata, as an ancillary variable to enhance the property interpolation. This approach helps control the overall strata-bound trend of the properties, thereby improving the validity of the interpolation results.

The stratigraphic interface is represented by a specific equipotential surface of the SPF $P\left(\mathit{u}\right)$. The SPF values with a stratum change gradually from bottom to top, reflecting the temporal order of stratum deposition [35,39]. When the number of sampling points for property values is limited, an interpolation function can be established, assisted by the known SPF values as ancillary variables. The SPF can be estimated using an RBF interpolant, allowing the final interpolation function to be expressed as follows:

$$\widehat{f}\left(\mathit{u}\right)=\sum _{i=1}^n{w}_i\phi \left(‖\mathit{u}-{\mathit{u}}_i‖\right)+\alpha \widehat{P}\left(\mathit{u}\right)+\beta$$

(4)

where $\phi \left(‖\mathit{u}-{\mathit{u}}_i‖\right)$ is the RBF and $\widehat{P}\left(\mathit{u}\right)$ is the estimated SPF function. There are $n+2$ parameters to be solved; i.e., $w_i$(i = 1, 2, …, n) are the weighting coefficients of the RBF, α is the weighting coefficient of the SPF, and β is a constant term to reduce interpolation error when the point to be interpolated is far from the sample points. With n coordinates and corresponding property values from sampling points, to obtain n + 2 weighting coefficients requires the orthogonality conditions of the basis functions, as follows:

$${\sum }_{i=1}^n{w}_{i}P\left({\mathit{u}}_{\mathit{i}}\right)=0,{\sum }_{i=1}^n{w}_i=0$$

(5)

Thus, there are $n+2$ equations, which can be expressed in matrix form. The solution for the $n+2$ parameters can be determined, allowing the interpolation function to be specified as follows:

$$\left[\begin{array}{cc}\mathit{A}& \mathit{B}\\ {\mathit{B}}^{\mathrm{T}}& \mathbf{0}\end{array}\right]·\left[\begin{array}{c}\mathit{W}\\ \mathit{\Theta }\end{array}\right]=\left[\begin{array}{c}\mathit{F}\\ \mathbf{0}\end{array}\right]$$

(6)

$$\left[\begin{array}{c}\mathit{W}\\ \mathit{\Theta }\end{array}\right]={\left[\begin{array}{cc}\mathit{A}& \mathit{B}\\ {\mathit{B}}^{\mathrm{T}}& \mathbf{0}\end{array}\right]}^{-1}·\left[\begin{array}{c}\mathit{F}\\ \mathbf{0}\end{array}\right]$$

(7)

where $\mathit{A}=\left[\begin{array}{cc}\begin{array}{cc}0& {\phi }_{21}\\ {\phi }_{12}& 0\end{array}& \begin{array}{cc}\cdots & {\phi }_{n1}\\ \cdots & {\phi }_{n2}\end{array}\\ \begin{array}{cc}⋮& ⋮\\ {\phi }_{1n}& {\phi }_{2n}\end{array}& \begin{array}{cc}\ddots & ⋮\\ \cdots & 0\end{array}\end{array}\right],\mathit{B}=\left[\begin{array}{c}\begin{array}{cc}P\left({\mathit{u}}_1\right)& 1\\ P\left({\mathit{u}}_2\right)& 1\end{array}\\ \begin{array}{cc}⋮& ⋮\\ P\left({\mathit{u}}_n\right)& 1\end{array}\end{array}\right],\mathit{W}=\left[\begin{array}{c}{w}_1\\ ⋮\\ w_n\end{array}\right],\mathit{\Theta }=\left[\begin{array}{c}\alpha \\ \beta \end{array}\right],\mathit{F}=\left[\begin{array}{c}{f}_1\\ ⋮\\ f_n\end{array}\right]$.

2.2. Verification Experiments

2.2.1. Experiment Models

To verify the effectiveness of the proposed MRBF interpolation method, a two-dimensional (2D) experimental model was generated to compare this method with existing common interpolation methods. The property field of the experimental model is defined by the equation $f\left(x,y\right)=\sin (\frac{x}{50})+5\cos (\frac{y}{50})$, where the region scope is $0\le x\le 150 \mathrm{a}\mathrm{n}\mathrm{d} 0\le y\le 100$. The property values range from −2.08 to 6, with an average value of approximately 2.9245. Based on the distribution of the property field, a distance function $P\left(x,y\right)=\sqrt{{\left(x-75\right)}^2+{\left(y+80\right)}^2}$ was selected as the ancillary model for multivariate interpolation. The experimental model and the ancillary model are illustrated in Figure 1. The correlation coefficient r between f and P is −0.9804, indicating a significant negative correlation, validating their use for co-interpolation.

A set of sample points with known property values f_i were randomly generated, as shown in Figure 2a. Additionally, 31 × 21 sampling points for the ancillary variable P_j were obtained using uniform sampling with a 5 m interval, as shown in Figure 2b. Their coordinates and ancillary variable values are then substituted into the RBF interpolation formula to calculate P(u). Grid points served as the points to be interpolated, which were compared to the true values.

Thirty experiments, each with 20 randomly generated sampling points and corresponding property values, were conducted. The MRBF interpolation results using four radial basis functions, i.e., cubic, linear, Gaussian, and thin plate, were compared, as shown in Figure S1. The results clearly demonstrate that the cubic function consistently outperforms the other three RBFs, both in terms of maximum and average errors.

2.2.2. Interpolation Results

By comparing the original model, it indicates that there is no significant difference between the interpolation results and the true values, as shown in Figure 3a. To further evaluate the accuracy, Figure 3b presents a scatter plot with the true values as the horizontal axis and the interpolation results as the vertical axis. The line $\widehat{f}=f$ signifies perfect consistency between true and interpolated values. Most points are distributed near the line $\widehat{f}=f$, indicating a high degree of consistency. The interpolation results tend to be slightly smaller than the true values at lower values and slightly larger at higher values. Figure 3c illustrates the distribution of absolute errors. The maximum absolute error is 0.39, with an average error of 0.041, demonstrating a high degree of interpolation accuracy. Data points with large absolute errors are primarily located in the lower left corner, corresponding to areas with sparse sampling points.

2.2.3. Influence of Number of Sampling Points

In addition to the proposed MRBF interpolation algorithm, several other commonly used interpolation algorithms, i.e., IDW, OK, and DSI, were employed to assess the effectiveness of the MRBF algorithm. To demonstrate the effectiveness of the proposed method, 10, 20, and 30 points were randomly sampled, as shown in Figure S2, and an interpolation method comparison was conducted, as shown in Figure S3.

The performance of these methods was evaluated using three evaluation metrics of maximum absolute bias (MABS), mean absolute error (MAE), and mean square error (MSE).

MABS measures the largest absolute difference between interpolated and true values, indicating the worst-case error.

$$\mathrm{MABS}=\max (\left|f_i-f_{i0}\right|)$$

(8)

MAE computes the average absolute differences between interpolated and true values, reflecting overall accuracy.

$$\mathrm{MAE}=\frac{1}{N}\sum _{i=1}^N\left|f_i-f_{i0}\right|$$

(9)

MSE calculates the average of the squared differences between interpolated and true values, emphasizing larger errors.

$$\mathrm{MSE}=\frac{1}{N}\sum _{i=1}^N{\left(f_i-f_{i0}\right)}^2$$

(10)

Table S2 reveals that the MRBF method proposed in this study consistently outperforms other methods across different numbers of sampling points. Even with a decrease in sampling points, the MRBF method still maintains robust interpolation performance. Overall, the interpolation methods are sorted by performance as follows: MRBF > OK > DSI > IDW. Notably, MRBF demonstrates superior interpolation performances even with a limited number of sampling points. The variogram for OK is obtained by automatically fitting a spherical model and the fitting parameters are provided in Table S3. The kriging interpolation exhibits instability when sampling points are limited or when the variogram fitting is inaccurate, reducing its overall performance. While RBF ranks second after MRBF with 30 and 20 sampling points, it shows poorer performance with only 10 sampling points, with an MAE of 6.529. This discrepancy may stem from RBF’s stringent requirement for uniformly distributed sampling points.

2.2.4. Influence of Distribution of Sampling Points

To investigate whether MRBF faces similar challenges related to the spatial distribution of sampling points, experiments were conducted with 18 sampling points using both random and uniform sampling methods. Figure S4 illustrates the spatial distribution of these sampling points, and Figure S5 presents the interpolation results for each algorithm under both scenarios.

As shown in Table S4, when the sampling points are uniformly distributed, the accuracy of RBF interpolation significantly improves compared to random distribution. This underscores RBF’s requirement for uniformly distributed sampling points. In contrast, the MRBF interpolation shows only marginal improvement, indicating that the ancillary variable integration in MRBF relaxes the uniformity requirement for sampling points and maintains robust interpolation results even with sparse distributions. Among all the methods, kriging’s performance suffers from poor variogram fitting over uniformly distributed sampling points, leading to unreliable results. Under random sampling, although kriging exhibits the smallest MABS value, its MAE and RMSE values are larger than those of MRBF, placing it second after MRBF. IDW and DSI methods exhibit similar performance, showing less sensitivity to sampling point uniformity.

2.3. Dataset

The data source for this case study is the detailed geotechnical investigation dataset of a building site at the stage of construction drawing design. As shown in Figure 4, the construction site is located beside South Zhengyang Road in Hengyang City. Its geomorphology is characterized by the Xiangjiang III terrace landform, relatively flat and open, with ground elevations ranging from 84.45 m to 87.17 m. The site is situated in the South China Fracture Zone, within the Hengyang Basin, which lies in the middle and lower part of the Yangtze River Faulted Depression. This area is a Cretaceous–Tertiary terrestrial stabilized basin. The site’s stratigraphy belongs to the Tertiary inland lake sedimentary, predominantly composed of clastic rocks. The regional tectonics are dominated by the Xishan Stage, with two main groups that are NNE- and NNW-oriented. There are no traces of large faults passing through the site or its neighboring lots, and recent tectonic activity, especially newly active rupture structures, is not evident, indicating the site is in a relatively stable state. According to geological drilling results, the exposed strata within the exploration depth of this construction site include plain fill soil (${\mathrm{Q}}_4^{\mathrm{ml}}$), silty clay (${\mathrm{Q}}_4^{\mathrm{a}\mathrm{l}}$), gravelly sand (${\mathrm{Q}}_4^{\mathrm{a}\mathrm{l}}$), clay (${\mathrm{Q}}_4^{\mathrm{e}\mathrm{l}}$), and the Lower Tertiary (E) mudstone.

The dataset used for modeling includes 37 boreholes (Figure 5a), 16 engineering geological profiles (Figure 5b), a location map of buildings, and an experimental report of physical and mechanical properties of 41 soil samples. A statistical summary of the properties at the sampling points is presented in

Table 1. Statistical summary of the properties at the sampling points.

Stratum	Number of Sampling Points	Natural Moisture Content			Water Saturation			Natural Pore Ratio
		Mean (%)	Std. (%)	Skewness	Mean (%)	Std. (%)	Skewness	Mean	Std.	Skewness
Silty clay	5	27.14	1.54	−0.556	88.0	5.74	−0.566	0.851	1.54	0.785
Clay	18	42.67	8.10	0.187	98.6	2.64	−2.18	1.17	0.222	0.0919
Mudstone	14	34.4	5.61	0.443	97.2	3.55	−0.875	0.952	0.153	0.183

. Five buildings are planned for construction, with exploration completed for buildings 1#, 2#, and 3#, as shown in Table S5. Figure 6 shows an example of an engineering geologic profile 2-2’.

The surface building model is an inclined photogrammetric dataset in OSGB (open scene graph binary) format. Figure 7 shows the surface buildings in the vicinity of the site. The total planned land area is 18,306.1 m², and the total building area is 95,211.6 m².

3. Results

3.1. Stratigraphic Structure Models

Geological profiles combine original engineering geological investigation data with expert interpretation, making them suitable for the 3D modeling of stratigraphic contact surfaces. In this study, the 16 profiles interpreted by engineer geologists were processed to obtain a series of stratum boundaries. These boundaries were then stitched up according to the strata. Subsequently, the surface model with the stratigraphic undulation pattern was constructed based on the DSI algorithm. Finally, the boundaries were stitched together to generate a solid model. The 3D geological structure modelling process is shown in Figure 8. The volume and average thickness of each stratum were calculated, as shown in

Table 2. Volume and thickness of each stratum.

Stratum	Volume	Average Thickness	Thickness Range
Plain fill soil (${\mathrm{Q}}_4^{\mathrm{ml}}$)	10,901.9	0.72	0.5–3.4 m
Silty clay (${\mathrm{Q}}_4^{\mathrm{a}\mathrm{l}}$)	54,567.2	3.60	0.8–12.3 m
Gravelly sand (${\mathrm{Q}}_4^{\mathrm{a}\mathrm{l}}$)	74,792.8	4.94	0.3–16.5 m
Clay (${\mathrm{Q}}_4^{\mathrm{e}\mathrm{l}}$)	135,235	8.93	1.2–17.1 m
Fully weathered mudstone (E)	170,587	11.3	Partially undissected; maximum control stratum thickness: 15.5 m.
Strongly weathered mudstone (E)	298,465	20	Largely undissected.

.

Groundwater at the site consists primarily of loose rock pore water, predominantly stored within the gravel sand strata. It is recharged by the infiltration of atmospheric precipitation and lateral flows from surrounding aquifers, discharging towards lower-lying areas, particularly to the east and northeast border of the site, exhibiting dynamic variability and generally adequate water content. The groundwater levels (GWLs) in the boreholes range between 1.9 m and 20.0 m in depth, with elevations spanning from 66.91 m to 83.45 m. Figure 9 illustrates that the GWL surface aligns with the bottom surface of the gravel–sand stratum, highlighting it as the principal aquifer and flow conduit on the site. Groundwater flow generally trends from southwest to northeast.

3.2. Stratigraphic Property Models

The solid model was gridded into cells with size of 1 m × 1 m × 0.2 m, where the cell properties were filled with the values at the cell center. Three properties, i.e., natural pore ratio, natural moisture content, and water saturation, were selected for modeling according to the experimental report of physical and mechanical properties of soil samples. The SPF values of the stratum surfaces were set to 10, 20, 30, 40, and 50 from top to bottom. The SPF values at points between the stratum surfaces were obtained using the IDW interpolation method (Figure 10). The SPF model was uniformly sampled under the structural constraints of the strata and used as a covariate to assist in the main variable interpolation. The MRBF method was applied to interpolate the properties of the entire model. Since the properties vary across different strata while being similar within the same stratum, the interpolation was conducted under the constraints of the stratum. Specifically, the interpolation was performed separately for the silty clay, clay, and mudstone strata. The interpolation results are visualized in Figure 11, and the histogram illustrates the distribution of the properties values of each stratum, as shown in Figure 12. The distributions of natural moisture content and natural pore ratio are closer to normal, especially for clay and mudstone, which exhibit low skewness (<0.3) and relatively symmetric shapes. In contrast, silty clay shows a slight right skew. On the other hand, water saturation significantly deviates from normality, exhibiting high negative skewness (<−1) with sharp peaks and long left tails, particularly in silty clay and clay. This suggests that water saturation is likely constrained by boundary conditions (i.e., an upper limit of 100%), whereas natural pore ratio and moisture content are influenced by more random processes, resulting in distributions that are closer to a Gaussian shape. The distribution of the moisture content and pore ratio revealing a similar pattern, suggesting a possible correlation between them.

3.3. Pile Foundation Bearing Capacities

Based on the 3D geological models and the engineering geological conditions of the site, combined with the engineering scale of the proposed buildings, it is recommended that a prestressed pile foundation be adopted, using the strongly weathered mudstone stratum as the bearing stratum. The expected length of the piles is 20.0~25.0 m (measured from the basement floor), and the diameter of the piles is 400~500 mm. The pile bearing capacity is determined according to the analysis and statistical results of geotechnical in situ and indoor tests. Drawing insights from Chinese standards [40,41], construction industrial standards [42], and regional past engineering experiences [43], the pile bearing capacity calculation formula is expressed as follows:

$$Q_{uk}=Q_{sk}+Q_{pk}=u\sum q_{sik}{l}_i+\alpha q_{pk}{A}_P$$

(11)

where $Q_{sk}$ is the standard value of total ultimate lateral resistance, $Q_{pk}$ is the standard value of total ultimate end resistance, μ is the perimeter of pile body, $q_{sik}$ is the standard value of ultimate lateral resistance of the ith stratum of soil on the side of the pile, $l_i$ is the thickness of the ith stratum of soil on the side of the pile, $\alpha$ is the correction coefficient of the pile end resistance (taking 0.9), $q_{pk}$ is the standard value of the ultimate end resistance, and $A_P$ is the area of the pile end. The values of $q_{sik}$ and $q_{pk}$ generally derived from local experiences, as shown in Table S5.

The plain fill soil and silty clay strata will be completely excavated during construction, so their lateral resistance is not considered. Except for the strongly weathered mudstone stratum, other strata cannot serve as the bearing stratum due to inadequate capacity. Therefore, only the strongly weathered mudstone stratum is used to calculate the standard value of ultimate end resistance, taken as half of this standard value for estimation purposes. The actual characteristic value of the bearing capacity of a single pile should be determined through pile load tests and adjusted based on the test results if necessary.

The pile length was set to 25.0 m and measured from the basement floor at an elevation of 78m. As shown in Figure 13a, 491 piles were arranged in a full pile configuration, with a diameter of 500 mm and a spacing of 2000 mm between piles. The thickness of the pile foundation through each soil stratum was calculated from the coordinates of the intersection points between the pile foundation and each horizon surface. The standard value of the bearing capacity of each pile is illustrated in Figure 13b. Tables S7–S9 present the characteristic values of bearing capacity for each pile. The characteristic value of the bearing capacity of the building is taken as the minimum value of the pile bearing capacity.

Table 3. The number of pile foundations and the bearing capacity characteristic values of each building.

Building No.	Number of Piles	Characteristic Value of Pile Bearing Capacity (kN)
1#	236	2162.48
2#	166	2191.08
3#	89	2153.81

shows the characteristic values of the pile bearing capacity of each building, all of which exceed the maximum axial force of a single column, 2000 kilonewtons (kN). This indicates that the bearing capacity of the prestressed piles is adequate. Using conventional calculations based on the average thickness of each soil stratum a characteristic value of pile bearing capacity of 2333.11 kN was obtained, which is higher than the values calculated separately. The construction is designed according to the calculated value; if the actual bearing capacity is insufficient, it may lead to potential safety risks. For example, if the maximum shear stress is 2200 kN, the pile foundation bearing capacity calculated using the conventional calculation seems sufficient. In fact, some of the piles’ bearing capacity does not reach 2200 kN, potentially causing safety risks and resource wastage. The proposed method of calculating the bearing capacity of each pile foundation separately mitigates this issue, allowing for a clear identification of locations with weak bearing capacity.

4. Discussions

The proposed MRBF interpolation algorithm incorporates an ancillary variable to control the overall property trends. The results of the verification experiments, conducted using a 2D model, demonstrate improved accuracy compared to the other commonly used interpolation methods. MRBF maintains robust performance with fewer sampling points and mitigates uniformity requirement drawback of RBF, enabling itself as an effective and stable interpolation technique. Moreover, this algorithm allows for the inclusion of any other relevant variable. When the number of sampling points of the primary variable is insufficient, another variable with more sampling points, related to the primary variable, can be selected as an ancillary variable to assist in interpolation. This can improve interpolation accuracy to some extent and make the interpolation of variables with limited sampling points possible.

Engineering geological profiles were used to model the stratigraphic structure and constrain the spatial distribution of properties within the strata. Boreholes and testing samples are crucial as they provide direct and accurate information about stratigraphy and subsurface properties. The influence propagation of borehole data in space is determined by both spatial proximity and stratigraphic relationships. Specifically, the influence of borehole data decreases with increasing distance, with closer points exerting a stronger similarity, as reflected in the primary variable. Additionally, subsurface properties within the same stratum are closely related, exhibiting similar properties across the stratum. This relationship was incorporated using an ancillary variable, which helps capture stratigraphic continuity and ensures more accurate property modeling beyond the borehole locations.

Ordinary kriging is the cornerstone algorithm of geostatistics, offering a robust approach that combines many desirable features of alternative geospatial inference methods [44,45]. Li and Heap (2011) found that the interpolation performance of OK is superior to that of IDW, with kriging being more suitable for spatial interpolation than non-geostatistical methods [46]. Cokriging has the advantages of OK, while also incorporating ancillary information to improve estimates of a primary variable. Stein and Corsten (1991) indicated that cokriging performs better when covariates are highly correlated with the primary variable [47]. Schiavo (2024) demonstrated that cokriging is the most suitable approach for reproducing water table levels when local terrain elevation is used as an ancillary variable, overcoming challenges in the spatial estimation of historical water table levels [48]. However, in some cases, the assumptions underlying kriging interpolation may not hold, or when the sampling points are few and unevenly distributed, the semi-variogram may not be accurately fitted, leading to significant errors in the kriging interpolation results. While the OK variogram function parameters were fitted in our comparison experiments, there was no specific discussion of how their potential changes might affect the resulting variograms and subsequent spatial predictions, as shown in Schiavo (2023) [49]. Given the importance of these parameters in geostatistical modeling, future studies on OK could assess the impact of variogram parameter changes on the model by performing sensitivity analyses. In contrast, the MRBF approach reduces dependence on sampling density, making it a robust tool for subsurface property modeling.

The integrated surface building and subsurface engineering geological models are visualized in Figure 14. Integrating the processed surface buildings with the 3D geological models enables the expression of integrated spatial information in a digital city. The development of urban digital twins relies heavily on accurate subsurface geological models. By combining detailed geological modeling with urban digital twin frameworks, improved subsurface geological analysis and urban planning can be achieved, contributing to the creation of smarter, safer, and more sustainable urban environments [50,51].

5. Conclusions

This study proposes a 3D geological property modeling method that incorporates an SPF model to assist property interpolation with a limited number of sampling points. Property interpolation was performed based on the MRBF method under stratigraphic model constraints. Additionally, the characteristic values of pile foundation bearing capacity for each building were calculated, providing a more accurate reference for subsequent construction. Finally, an integrated representation of surface and subsurface space was achieved, offering enhanced functionality for managing, mining, analyzing, visualizing, and sharing geological models. The main conclusions are as follows:

(1)

The proposed MRBF method effectively utilizes ancillary variable for multivariate interpolation, enabling relatively accurate interpolation with a few property sampling points. Experimental results from a 2D model indicate its high accuracy. The MRBF method outperforms several commonly used methods, i.e., IDW, DSI, and kriging, particularly when the number of sampling points is limited. Additionally, the MRBF method’s performance remains stable regardless of the uniformity of sampling point distribution, making it robust for sparse sampling.

(2)

A geological structural model was constructed based on profiles. The bearing capacity of each pile foundation was calculated individually rather than using an overall estimate based on average stratum thickness, providing a more precise reference for construction. The GWL surface aligns closely with the bottom surface of the gravel sand stratum, indicating that the site’s groundwater is primarily loose rock pore water stored in the gravel sand stratum’s pore spaces.

(3)

By interpolating properties using the MRBF method under structural model constraints, property models, consistent with stratigraphic deposition, were conducted. Integrating the processed surface buildings with the 3D geological models enables the expression of integrated spatial information in a digital city.

Despite these advancements, this study has limitations concerning the MRBF method and the deterministic model. The MRBF algorithm requires that ancillary variables accurately reflect the distributional trends of the main variables for effective multivariate interpolation. Otherwise, the inclusion of an ancillary variable may introduce errors and reduce performance. Furthermore, the profile-based stratigraphic structure modeling and interpolation-based property modeling methods employed in this study are limited in their ability to quantify uncertainty. Future work should focus on developing a probabilistic 3D geological modeling method (e.g., multiple-point statistics method [52]) to effectively quantify uncertainty by simulating multiple possible geological models.