Three-Dimensional Stratigraphic and Attributive Collaborative Modeling in Urban Engineering Construction

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 framework 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 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 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 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. By accurately calculating pile foundation bearing capacity based on detailed stratigraphic models, this study minimizes geotechnical risks and enhances construction safety. Finally, the surface building inclined photogrammetric models are integrated with the subsurface 3D engineering geological models.

2. Method

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 an attribute filed 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 auxiliary 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)+\widehat{P}(\mathit{u})+\mathrm{\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\text{，}{\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}}^T& 0\end{array}\right]·\left[\begin{array}{c}\mathit{W}\\ \Theta \end{array}\right]=\left[\begin{array}{c}\mathit{F}\\ 0\end{array}\right]$$

(6)

$$\left[\begin{array}{c}\mathit{W}\\ \Theta \end{array}\right]={\left[\begin{array}{cc}\mathit{A}& \mathit{B}\\ {\mathit{B}}^T& 0\end{array}\right]}^{-1}·\left[\begin{array}{c}\mathit{F}\\ 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],\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)=\mathit{sin}({\displaystyle \frac{x}{50}})+5\mathit{cos}({\displaystyle \frac{y}{50}})$, where the region scope is $0\le x\le 150and0\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 5m 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, ordinary Kriging (OK), and DSI, were employed to assess the effectiveness of the MRBF algorithm. To demonstrate the superiority of our proposed method, we randomly sampled 10, 20, and 30 points, as shown in Figure S2, and conducted the interpolation method comparison, as shown in Figure S3.

We evaluated these methods 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}=\mathit{max}(\left|f_i-f_{i0}\right|)$$

(8)

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

$$\mathrm{MAE}={\displaystyle \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}={\displaystyle \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 > Kriging > DSI > IDW. Notably, MRBF demonstrates superior interpolation performances even with a limited number of sampling points. 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, we conducted experiments 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 S3, 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 value 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.

In summary, the proposed MRBF interpolation algorithm incorporates an ancillary variable to control the overall trends. The results of the verification experiments, conducted using a 2D model, demonstrate superior accuracy compared to the other commonly used interpolation methods. MRBF maintains robust performance with fewer sampling points, and mitigates RBF’s uniformity requirement drawback, 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 main variable is insufficient, another variable with more sampling points related to it can be selected as an ancillary variable to assist interpolation, which can improve the interpolation accuracy to a certain extent and make the interpolation of variables with limited sampling points possible.

3. Case Study

3.1. Dataset

The data source for this case study is the detailed geotechnical investigation dataset of a building site at the stage of construction drawings 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.45m to 87.17m. 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 of 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 status. According to geological drilling results, the exposed strata within the exploration depth of this construction site include plain fill soil (Q₄^ml), silty clay (Q₄^al), gravelly sand (Q₄^al), clay (Q₄^el), and the Lower Tertiary (E) mudstone.

The dataset used for modeling include 37 borehole histograms, 16 engineering geological cross-sections (Figure 5a), a location map of buildings and exploration points (Figure 5b), and an experimental report of physical and mechanical properties containing 41 soil samples. Five buildings are planned for construction, with exploration completed for buildings 1#, 2#, and 3#, as shown in Table S4. Figure 6 shows an example of an engineering geologic cross-section 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.2. Geological Structural Modeling

Geological cross-sections combine of original engineering geological investigation data with expert interpretation, making them suitable for 3D modeling of stratigraphic contact surfaces. In this study, the 16 cross-sections 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 1.

Groundwater at the site consists primarily of loose rock pore water, predominantly stored within the gravel sand layers. It is recharged by 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 layer, highlighting it as the principal aquifer and flow conduit on the site. Groundwater flow generally trends from southwest to northeast.

3.3. Geological Property Modeling

The solid model was grided into cells with size of 1m×1m×0.2m, 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. We employed the MRBF method to interpolate the properties of the entire model. Since the properties vary among different strata while being similar within the same stratum, this interpolation was conducted under the constraints of the stratum. Specifically, the interpolation was performed separately for the silty clay, clay, and mudstone layers. 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 distribution of the moisture content and pore ratio revealing a similar pattern, suggesting a possible correlation between them.

4. Discussions

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 to adopt prestressed pile foundation, using the strongly weathered mudstone layer as the bearing layer. 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 layer of soil on the side of the pile, $l_i$ is the thickness of the ith layer 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 layers will be completely excavated during construction, so their lateral resistance is not considered. Except for the strongly weathered mudstone layer, other layers cannot serve as the bearing stratum due to inadequate capacity. Therefore, only the strongly weathered mudstone layer is used to calculate the standard value of ultimate end resistance, taken as half of this standard value for estimation purpose. 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 is set to 25.0 m, 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 layer 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 S6–S8 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 2 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 kiloNewton (kN). This indicates that the bearing capacity of the prestressed piles is adequate. Using conventional calculations based on the average thickness of each soil layer, we obtained a characteristic value of pile bearing capacity of 2333.11 kN, 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 reached the 2200 kN, potentially causing safety risks and resource wastage. Our 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.

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.

5. Conclusions

This study proposes a 3D geological property modeling method that incorporates SPF model to assist property interpolation with a limited number of sampling points. We achieved property interpolation based on the MRBF method under stratigraphic model constraints. Additionally, we calculated the characteristic values of pile foundation bearing capacity for each building, providing a more accurate reference for subsequent construction. Finally, we achieved an integrated representation of surface and subsurface space. 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 cross-sections. We calculated the bearing capacity of each pile foundation individually rather than using an overall estimate based on average layer thickness, providing a more precise reference for construction. The GWL surface aligns closely with the bottom surface of the gravel sand layer, indicating that the site’s groundwater is primarily loose rock pore water stored in the gravel sand layer’s pore spaces.

(3) By interpolating property using the MRBF method under structural model constraints, we obtained property models consistent with stratigraphic deposition. 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 integrated model. The MRBF algorithm requires a significant correlation between primary and ancillary variables for meaningful multivariate interpolation. Without a strong correlation, adding an ancillary variable may introduce errors, reducing its performance. Future work should establish a standard for the required correlation strength. The integrated model could also offer more functionalities for managing, mining, analyzing, visualizing, and sharing geological models.