A soil property three-dimensional mapping method, device and storage medium based on integrated nested Bayesian and Laplace approximation

Abstract

The application discloses a kind of based on integrated nested bayesian and laplace approximation soil attribute three-dimensional mapping method, device and storage medium, method includes: constructing three-dimensional unstructured grid;Optimal feature set is screened to environmental covariate;Establish three-dimensional bayesian spatial statistical model based on INLA-SPDE, soil attribute is characterized as Gaussian random field in three-dimensional space and is carried out efficient inference;Realize the three-dimensional space continuous prediction of soil attribute and generate uncertainty distribution chart;Soil carbon storage can be further estimated.The device includes processor and memory, for executing the method.The application breaks through the limitation of traditional two-dimensional and layered three-dimensional prediction method, realizes high-precision, interpretable modeling and embedded uncertainty evaluation in true three-dimensional space, with high calculation efficiency, and provides a powerful technical tool for precision agriculture, environmental management and carbon cycle research.

Description

Technical Field

[0001] This invention belongs to the field of digital soil mapping and spatial information technology, specifically relating to a method, apparatus and storage medium for three-dimensional mapping of soil properties based on integrated nested Bayesian and Laplace approximation. Background Technology

[0002] Soil is a core component of the Earth's critical zone, and its spatial distribution information, especially soil organic carbon, is crucial for interpreting soil formation processes, assessing soil quality, simulating the global carbon cycle, and formulating land management policies to address climate change. Traditional digital soil mapping methods are primarily based on a two-dimensional plane, making spatial predictions by establishing statistical relationships between soil properties and surface environmental factors. However, soil is a typical three-dimensional continuum, and its properties often exhibit much greater variability in the vertical direction than in the horizontal direction. This vertical heterogeneity is a key indicator for assessing soil development, classifying soil types, and simulating soil processes. The literature (Liang, Z., Chen, S., Yang, Y., Zhou, Y., & Shi, Z. (2019). High-resolution three-dimensional mapping of soilorganic carbon in China: Effects of SoilGrids products on national modeling. Science of the total environment, 685, 480-489; Liu, F., Zhang, GL, Song, X., Li, D., Zhao, Y., Yang, J., ... & Yang, F. (2020). High-resolution and three-dimensional mapping of soil texture of China. Geoderma, 361, 114061; Xie, Xianli, Xia, Chengye, Yin, Biao, Li, Anbo, Li, Kaili, & Pan, Xianzhang. Research progress on three-dimensional spatial soil inference and soil model construction. Acta Pedologica Sinica, 62(1), 14-28.) points out that with the rapid development of three-dimensional geoscience information technology, three-dimensional soil mapping has become a research frontier.

[0003] Literature (Liu, Y., Qian, Y., Zhu, Y., Xu, W., Wei, G., Huang, J., ... &Ma, Q. (2025). Spatial estimation of large-scale soil salinity using enhanced inverse distance weighting method and identifying its driving factors. Agricultural Water Management, 317, 109645; Li, Z., Tao, H., Zhao, D., & Li, H. (2022). Three-dimensional empirical Bayesian kriging for soil PAHsinterpolation considering the vertical soil lithology. Catena, 212, 106098; Gao, Z., Peña-Arancibia, JL, & Siyal, AA (2025). Three-dimensional soilsalinity mapping with uncertainty using Bayesian Hierarchical Modelling, Random Forest Regression and remote sensing data. Agricultural WaterManagement, 309, (109318) discloses representative existing 3D spatial prediction methods, including: 3D geostatistical methods such as 3D Inverse Distance Weighted Method (3D-IDW), 3D Kriging, and 3D Empirical Bayesian methods (such as 3D-EBK). While 3D-IDW is simple and easy to use, it assumes that geographic objects are isotropic, leading to a significant decrease in prediction accuracy when data distribution is uneven or outliers are present, and it fails to provide uncertainty information. Methods like 3D-EBK can handle spatial autocorrelation and provide uncertainty metrics, but their computational complexity increases dramatically with the amount of data, making them inefficient for large-scale 3D data processing and sensitive to initial model parameters. Hierarchical 3D machine learning methods, such as hierarchical 3D Random Forest (3D-RF), are also discussed. This method establishes prediction models for different soil layers at different depths, enabling it to handle high-dimensional data and complex nonlinear relationships.However, this method is essentially a stack of multiple two-dimensional models, which fails to achieve true three-dimensional continuous modeling. It is difficult to characterize the continuous change and interaction of attributes in the vertical direction, and there is a risk of overfitting. The model has poor interpretability and is also difficult to provide reliable spatial uncertainty quantification.

[0004] In summary, existing three-dimensional soil mapping technologies suffer from the following common defects: (1) difficulty in achieving true three-dimensional continuous modeling and prediction; (2) low computational efficiency, making it difficult to handle large-scale three-dimensional data; (3) insufficient or missing ability to quantify the uncertainty of prediction results; and (4) complex models with poor interpretability. These limitations severely restrict the depth and breadth of three-dimensional soil mapping in scientific research and practical applications.Literature (Folly, CL, Konstantinoudis, G., Mazzei-Abba, A., Kreis, C., Bucher, B., Furrer, R., & Spycher, BD (2021). Bayesian spatial modeling of terrestrial radiation in Switzerland. Journal of Environmental Radioactivity, 233, 106571; Sun, XL, Minasny, B., Wang, HL, Zhao, YG, Zhang, GL, & Wu, YJ (2021). Spatiotemporal modeling of soil organic matter changes in Jiangsu, China between 1980 and 2006 using INLA-SPDE. Geoderma, 384, 114808; Wright, N., Newell, K., Lam, KBH, Kurmi, O., Chen, Z., & Kartsonaki, C. (2021). Estimating ambient The paper "Air Pollutant Levels in Suzhou through the SPDE approach with R-INLA" (International Journal of Hygiene and Environmental Health, 235, 113766) and "Assessing drivers of estuarine debris using a Bayesian spatial modelling approach (INLA-SPDE)" (Estuarine, Coastal and Shelf Science, 296, 108592) discloses an integrated nested Laplace approximation (INLA) algorithm, which is an emerging and efficient Bayesian inference algorithm that avoids the problems of large computation and poor convergence of traditional Markov chain Monte Carlo methods through a series of Laplace approximations.Stochastic partial differential equations (SPDEs) provide a mathematical framework for discretizing continuous Gaussian fields into Gaussian Markov random fields. Combined with sparse matrix techniques, this significantly reduces computational complexity and effectively improves the efficiency of model parameter fitting. The combination of INLA and SPDE has already demonstrated its advantages in two-dimensional spatial modeling. However, the innovative extension of this INLA-SPDE framework to three-dimensional space and its successful application to true three-dimensional continuous prediction and uncertainty quantification of soil properties is currently lacking. This invention aims to overcome the shortcomings of the aforementioned existing technologies. Summary of the Invention

[0005] This invention provides a three-dimensional mapping method for soil properties based on integrated nested Bayesian and Laplace approximation, in order to solve the technical problem that existing technologies cannot simultaneously achieve continuous prediction in three-dimensional space, high computational efficiency, and accurate uncertainty quantification.

[0006] This invention provides a method for three-dimensional mapping of soil properties based on integrated nested Bayesian and Laplace approximation, including: Multiple soil profile samples from the study area were obtained, along with the three-dimensional spatial coordinates (longitude, latitude, and depth), soil property measurements, and environmental covariates for each soil profile sample, to construct a sample set, which was then divided into a training sample set and a validation sample set. The 3D-INLA-SPDE model is constructed as follows: based on the three-dimensional spatial coordinates of the study area boundary and multiple soil profile samples, the constrained Delaunay triangulation method is used to obtain a three-dimensional unstructured mesh covering the entire study area; the SPDE method is used to discretize the three-dimensional Gaussian random field on the three-dimensional unstructured network to obtain a Gaussian Markov random field; based on the Gaussian Markov random field, the training samples are used as input to perform Bayesian inference based on the three-dimensional Bayesian model using the INLA algorithm, and the posterior distribution parameters of the model parameters are calculated. The 3D-INLA-SPDE model is trained using a training sample set. Based on the environmental covariates and three-dimensional spatial coordinates of the area under study, the model is trained to obtain the predicted values ​​of soil properties in three-dimensional space of the area under study, thereby realizing three-dimensional mapping of soil properties.

[0007] This invention utilizes the stochastic partial differential equation (SPDE) method to process continuous three-dimensional Gaussian random fields. Discretization is performed on the generated 3D triangular mesh, transforming it into a Gaussian Markov random field. This process leverages the sparsity of the triangular mesh to reduce computational complexity, thereby solving the problem of large-scale computational complexity in 3D modeling. NThe problem is that we further employ the INLA (Integrated Nested Laplace Approximation) algorithm for model inference. INLA directly and efficiently calculates model parameters (such as...) through a series of analytical approximations and numerical integrations. , The model outputs the posterior distribution of the variance, spatial correlation range, and hyperparameters, avoiding the convergence problems and computational burdens common in Markov Chain Monte Carlo (MCMC) algorithms. The model can output complete posterior distribution information, including the range parameter of the Matérn covariance function, the variance parameter, and parameters reflecting vertical autocorrelation.

[0008] Preferably, in the constrained Delaunay triangulation method covering the entire three-dimensional unstructured mesh of the study volume, the fineness of the three-dimensional unstructured mesh is controlled by key parameters, including offset parameters, truncation parameters, and longest edge length parameters. The offset parameters are used to define the mesh boundary extension distance, the truncation parameters are used to set the minimum allowable distance between vertices, and the longest edge parameter is used to set the maximum allowable boundary of the triangle.

[0009] More preferably, the fineness of the mesh is controlled by three key parameters: offset (Used to extend the mesh boundary to reduce edge effects, instance value is 6 km) cutoff (Define the minimum distance between vertices to prevent the generation of invalid triangles; instance value is 150 m) and maxedge (The maximum side length of the triangle is limited to control the mesh density; the instance value is 6 km). Optimizing these parameters is key to achieving a balance between computational efficiency and spatial approximation accuracy.

[0010] Preferably, the three-dimensional Bayesian model is: ,in, These are soil property observation values; This is an environmental covariate matrix constructed from environmental covariates; This is the fixed effects coefficient; To characterize a three-dimensional Gaussian random field with spatial dependence, For Gaussian measurement error, the posterior distribution parameters include The variance of S and the spatial correlation range.

[0011] Preferably, after obtaining the environmental covariates, a recursive feature elimination method is used to select features for the environmental covariates.

[0012] More preferably, a recursive feature elimination method is used to select features for environmental covariates, including: A k-fold cross-validation strategy is adopted, with RMSE and other indicators as evaluation criteria. The optimal subset of feature variables with the smallest prediction error for the target soil properties is iteratively selected from the obtained set of environmental covariates. Then, the soil property measurement values ​​and the corresponding values ​​of the optimal feature variable subset at the corresponding spatial locations are associated with the nodes of the constructed three-dimensional unstructured grid to form the structured dataset required for model training, i.e., the sample set.

[0013] Preferably, the three-dimensional Gaussian random field is discretized onto the three-dimensional unstructured network using the SPDE method to obtain a Gaussian Markov random field, including: Based on the parameters of a three-dimensional unstructured mesh (i.e., the mesh cells do not have a fixed arrangement rule or topological structure), a SPDE model is constructed by setting the Matérn covariance function, and the SPDE model is solved by the finite element method to obtain a discretized Gaussian Markov random field.

[0014] Preferably, the method for constructing a three-dimensional Gaussian random field includes: A Gaussian random field is constructed in the horizontal dimension using the stochastic partial differential equation method. The spatial correlation in the horizontal direction is obtained by defining an exponential model or a Gaussian model based on the constructed Gaussian random field in the horizontal dimension. The horizontal correlation and vertical correlation are separated and modeled using anisotropic covariance functions to reflect the characteristic differences of different soil layers or depth ranges. The vertical fluctuation range constraint or correlation decay rule is established through a cross covariance model to capture the indirect influence of deep soil on shallow parameters. The final model can simultaneously characterize the horizontal spatial variation and vertical distribution characteristics.

[0015] Preferably, when the soil property is soil organic carbon content, the soil bulk density and gravel content data are calculated based on the predicted soil property values ​​in three-dimensional space, combined with empirical formulas. The soil organic carbon density and soil organic carbon storage in the study area are then calculated using these formulas. The empirical formula for soil bulk density BD is as follows: (1) The soil organic carbon density (SOCD) of the i-th soil sample is: (2) The soil organic carbon (SOCS) storage of the i-th soil sample is: (3) In the formula, Let be the sampling depth of the i-th soil sample. This represents the side length of a cell in a three-dimensional unstructured mesh. Let be the gravel content of the i-th soil sample.

[0016] Preferably, based on the environmental covariates and three-dimensional spatial coordinates of the area under study, a trained model is used to obtain predicted values ​​of soil properties in three-dimensional space for the area under study, thereby achieving three-dimensional mapping of soil properties, including: This invention also provides three-dimensional spatial prediction and full-domain uncertainty quantification: using the 3D-INLA-SPDE model fitted by the above training, the soil properties of any unsampled point within the study area are spatially predicted. This process relies on inputting the coordinates of grid points throughout the three-dimensional space and their corresponding environmental covariate values ​​into the trained model, thereby generating a true three-dimensional continuous spatial distribution map of soil properties.

[0017] Specifically, the model, through its inherent Bayesian framework, provides posterior predicted distributions of soil properties for each spatial location in three-dimensional space to be predicted (typically defined by finely defined triangular mesh nodes). When generating a continuous spatial distribution map, the model systematically traverses all mesh nodes throughout the entire three-dimensional study domain. Based on spatial random effects and fixed environmental covariate effects, the model can map the continuous variations of soil properties in both the horizontal and vertical directions. This true three-dimensional mapping clearly reveals the spatial differentiation patterns of soil properties such as organic carbon content at different soil depths (e.g., within the range of 0–100 cm).

[0018] The spatial random effects provided by this invention refer to residual variations that cannot be explained by known covariates but are spatially correlated. It captures the effects of unknown, unmeasured, or inherently spatially continuous processes.

[0019] The fixed environmental covariate effects provided by this invention refer to the portions that can be explained linearly or nonlinearly by known environmental variables (such as temperature, altitude, precipitation, population density, etc.). These effects are "fixed," meaning they have a definite, global impact.

[0020] Furthermore, in the 3D-INLA-SPDE model, uncertainty quantification is directly based on the posterior marginal distribution of the latent Gaussian field. For each prediction node in the three-dimensional space, the INLA algorithm not only calculates the mean (or median) of the posterior prediction distribution—representing the best estimate of soil organic carbon content—but also calculates key statistics characterizing prediction accuracy.

[0021] Specifically, the standard deviation (or variance) of the posterior marginal distribution at each three-dimensional spatial location is used to quantify the uncertainty of local prediction results.

[0022] Furthermore, this study extracts 95% confidence intervals (i.e., the 2.5% and 97.5% quantiles of the posterior distribution) based on the posterior prediction results of the 3D-INLA-SPDE model, thus providing a probabilistic measure of the prediction range. These uncertainty indicators are then visualized across the entire three-dimensional space of the wetland to identify areas with lower prediction reliability, typically corresponding to locations with sparse data coverage or complex spatial heterogeneity. This explicit uncertainty output provides a more rigorous basis for risk assessment and decision-making in soil carbon mapping and monitoring.

[0023] On the other hand, the present invention also provides a soil property three-dimensional mapping device based on integrated nested Bayesian and Laplace approximation, comprising: one or more processors; a memory for storing one or more programs; when the one or more programs are executed by the one or more processors, the one or more processors implement the soil property three-dimensional mapping method based on integrated nested Bayesian and Laplace approximation.

[0024] On the other hand, the present invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the aforementioned method for three-dimensional mapping of soil properties based on integrated nested Bayesian and Laplace approximation.

[0025] Compared with the prior art, the beneficial effects of the present invention are as follows: This invention employs the SPDE method to discretize a three-dimensional Gaussian random field onto a three-dimensional unstructured network, obtaining a Gaussian Markov random field. For the first time, it successfully extends the INLA-SPDE framework to three-dimensional space, achieving accurate simulation of continuous and seamless changes in soil properties in both the horizontal and vertical directions. Compared to hierarchical prediction models, this method more realistically reflects the three-dimensional spatial structure of soil properties, and the prediction results are more consistent with natural laws.

[0026] The method incorporates Bayesian statistical uncertainty quantification capabilities: Based on Bayesian theory, this invention provides a complete posterior prediction distribution, capable of generating a comprehensive spatial distribution of uncertainty within the study area. This provides crucial scientific evidence, which is generally lacking in existing methods, for assessing the reliability of prediction results, guiding subsequent sampling design, and making risk decisions.

[0027] The computational efficiency is greatly improved: By discretizing the continuous spatial problem into a sparse matrix problem through SPDE and combining it with the efficient inference algorithm INLA, this invention effectively overcomes the bottleneck of excessive computational cost of traditional three-dimensional geostatistical methods, making it feasible to process large-scale, high-resolution three-dimensional soil data. The model running time is reduced by 37.5% compared with the 3D-RF model.

[0028] The prediction accuracy is significantly better than traditional methods: A case study in Poyang Lake wetland shows that the method of this invention has significantly better prediction accuracy for soil organic carbon than mainstream methods such as 3D-IDW, 3D-EBK, and 3D-RF, where R... 2 Compared with mainstream methods such as 3D-IDW, 3D-EBK and 3D-RF, it improves efficiency by 77.78%~116.22% and RMSE by 40.23%~44.68%, demonstrating its superior performance. Attached Figure Description

[0029] Figure 1 This is a schematic diagram of the study area according to an embodiment of the present invention, showing the relative positions of specific examples of the present invention.

[0030] Figure 2 This is a flowchart illustrating the overall method of an embodiment of the present invention, showing the complete workflow from data preparation to result output.

[0031] Figure 3 The image shows the result of variable selection using the recursive feature elimination algorithm in this embodiment of the invention, displaying the root mean square error of cross-validation corresponding to different numbers of variables.

[0032] Figure 4 This is a schematic diagram of a three-dimensional triangular mesh constructed in the Poyang Lake wetland study area in an embodiment of the present invention. The pink area represents the boundary of the study area, and the black dots represent the locations of sampling points.

[0033] Figures 5-8 This is a spatial distribution map of soil organic carbon at different depths in Poyang Lake wetland, predicted using the 3D-INLA-SPDE model, the existing 3D RF model, the 3D EBK model, and the 3D IDW model in this embodiment of the invention.

[0034] Figures 9-12 The diagram shows the upper and lower bounds and confidence intervals of the uncertainty distribution of soil organic carbon content predicted by the 3D-INLA-SPDE model at a 95% confidence level in this embodiment of the invention.

[0035] Figure 13 This is a comparison and verification chart showing the prediction accuracy of the model in this embodiment of the invention with that of the traditional method.

[0036] Figure 14 This is a ranking chart of importance based on SHAP values ​​in an embodiment of the present invention.

[0037] Figures 15-18 This is a partial dependency graph based on SHAP values ​​in an embodiment of the present invention. Detailed Implementation

[0038] The present invention will be described in detail below with reference to specific embodiments and accompanying drawings. This embodiment takes Poyang Lake Wetland (28°11'~29°51' N, 115°31'~117°06' E), China's largest freshwater wetland, as the study area. Figure 1 (a), (b), and (c) in the figures illustrate specific embodiments of the present invention, using three-dimensional spatial prediction of soil organic carbon (SOC) as an example. Figure 2 However, the scope of protection of this invention is not limited to this example.

[0039] (1) Soil sample data collection and environmental covariate preparation. The data preparation in this invention example includes the following specific steps: 1.1 Soil sample collection and SOC content determination: 1.1.1 Sampling time and scope: From December 2023 to January 2024, 127 soil profile samples (0-100 cm depth) were deployed and collected in the Poyang Lake wetland system.

[0040] 1.1.2 Sampling Method: A Rhino S1 soil drill was used for sampling. Each profile was divided into 10 layers: 0–10 cm, 10–20 cm, 20–30 cm, 30–40 cm, 40–50 cm, 50–60 cm, 60–70 cm, 70–80 cm, 80–90 cm, and 90–100 cm. The precise latitude and longitude coordinates of each sampling point were recorded using a portable GPS device.

[0041] 1.2 Laboratory Analysis: Soil samples were air-dried, ground, and sieved through a 2 mm sieve. The soil organic matter (SOM) content was determined using the potassium dichromate oxidation-external heating method. Subsequently, the SOM content was converted to SOC content (unit: g / kg) using the Bemmelen conversion factor (0.58). This yielded a dataset of 127 soil profiles (a total of 1270 stratified soil profile samples) containing three-dimensional spatial coordinates (longitude, latitude, and depth) and SOC content.

[0042] 1.3 Environmental covariate data collection and preprocessing: 1.3.1 Data Sources: Four major categories of environmental covariates were collected: soil properties, topography, climate, and biology (see Table 1 for details), including: Soil properties, such as sand content, silt content, clay content, pH, total nitrogen (TN), total phosphorus (TP), total potassium (TK), available nitrogen (AN), available phosphorus (AP), available potassium (AK), and cation exchange capacity (CEC), with a spatial resolution of 90 meters or 1 kilometer, are obtained from the Resource and Environmental Science Data Center of the Chinese Academy of Sciences (http: / / data.tpdc.ac.cn).

[0043] Topographic factors include elevation (Ele), slope, aspect, topographic moisture index (TWI), topographic location index (TPI), valley depth (VD), and multi-resolution valley floor flatness (MrVBF). These factors are primarily calculated using 90-meter resolution SRTM DEM data (https: / / srtm.csi.cgiar.org / ) in SAGA GIS or ArcGIS software.

[0044] Climate factors, such as mean annual temperature (MAT), mean annual precipitation (MAP), soil temperature (ST), and soil moisture content (SMC), are sourced from NASA (https: / / disc.gsfc.nasa.gov / ) and the Resource and Environmental Science Data Center of the Chinese Academy of Sciences (https: / / www.resdc.cn / ), with a resolution of 0.1° or 1 km.

[0045] Biological factors, such as Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), and Net Primary Productivity (NPP), are derived from satellite imagery products such as MODIS, with a resolution of 250 meters or 500 meters.

[0046] 1.3.2 Data Preprocessing: All environmental covariate data were preprocessed uniformly in ArcGIS 10.8 software. The steps included: converting the projection to the WGS 1984 UTM coordinate system; resampling all covariates to a spatial resolution of 90 meters using the nearest neighbor method (discrete data) or bilinear interpolation method (continuous data); and finally using the mask extraction tool to unify the spatial extent of all data to the boundary of Poyang Lake wetland.

[0047] (2) Feature variable selection and dataset partitioning optimize model input and improve prediction efficiency. The present invention implements the following steps: 2.1 Recursive Feature Elimination Algorithm for Variable Selection: 2.1.1 Implementation Method: Use the "..." function in R language (v4.3.3) caret "In the bag" rfe The function performs recursive feature elimination.

[0048] 2.1.2 Parameter settings: A 10-fold cross-validation strategy is adopted, with RMSE as the evaluation metric.

[0049] 2.1.3 Selection Process: The algorithm starts with all environmental covariates and iteratively performs the following operations: train the model - evaluate variable importance - remove the least important variables - retrain the model with the remaining set of variables. This process is repeated until the number of remaining variables is 5.

[0050] 2.1.4 Determining the optimal variable set: such as Figure 3 As shown, the model's cross-validation RMSE was minimized when the number of retained variables was 30 (2.3875 g / kg). Therefore, these 30 covariates were ultimately selected as the optimal feature set for subsequent modeling (the specific variables are as follows: Sand, AK, Silt, pH, TK, AN, SP, SG, AP, CEC, Clay, TP, TN, TE, MAP, SAT, SR, MAT, SMC, Slope, PC, MrVBF, TWI, VD, Aspect, LSF, TPI, NDVI, EVI, NPP).

[0051] 2.2 Dataset Partitioning: Data from 1270 soil stratified profile samples were randomly divided into a training set (1020 stratified profile samples for model training) and a validation set (250 stratified profile samples for model accuracy validation) at an 8:2 ratio using the Latin hypercube sampling method.

[0052] (3) 3D-INLA-SPDE model construction, training and prediction. This step is the core of the method of this invention, and the specific implementation is as follows: 3.1 3D Mesh Construction: 3.1.1 Method: Using R language... INLA Based on the three-dimensional spatial coordinates of the study area boundary and 1270 soil profile stratified samples, a constrained Delaunay triangulation was performed to generate an unstructured three-dimensional mesh.

[0053] 3.1.2 Key Parameter Settings: After optimization, the key mesh parameters set in this embodiment are as follows: offset (Grid boundary extension distance): 6 km, used to minimize boundary effects. cutoff (Minimum allowable distance between vertices): 150 m, to prevent the generation of triangles that are too small. maxedge (Maximum side length of the triangle): 6 km, used to control the mesh density in the outer extension region. The generated mesh effectively covers the entire study area (see...). Figure 4 And it maintains high accuracy in areas with dense sampling points.

[0054] 3.2 Definition and Training of 3D-INLA-SPDE Model: 3.2.1 Model Form: Establish a three-dimensional Bayesian model in the following form: (2) in, It refers to the observed soil organic carbon content; It is a matrix consisting of 30 preferred environmental covariates; It is the fixed effects coefficient; It is a Gaussian random field that characterizes the correlation in three-dimensional space; It is a Gaussian measurement error.

[0055] Specific embodiments of the present invention provide a method for constructing a three-dimensional Gaussian random field, including: Methods for constructing three-dimensional Gaussian random fields include: A Gaussian random field is constructed in the horizontal dimension using the stochastic partial differential equation method. By defining an exponential model or a Gaussian model based on the constructed Gaussian random field in the horizontal dimension, the spatial correlation in the horizontal direction is obtained. The spatial correlation refers to the characteristic that geographically close things or phenomena tend to have similar or related observations. Anisotropic covariance functions are used to separate and model the spatial correlation in the horizontal direction from the vertical correlation (calculated through the covariance function in the vertical direction) to reflect the differences in soil properties in different soil layers or depth ranges. By establishing the fluctuation range constraint or correlation decay rule in the vertical direction through the cross covariance model, the indirect influence of deep soil on shallow parameters is captured. The final generated model can simultaneously characterize the horizontal spatial variation and vertical distribution characteristics.

[0056] The anisotropic covariance function provided in this specific embodiment of the invention is implemented by applying a linear transformation matrix A to spatial coordinates. Its general form is: for a standard (isotropic) covariance function C(||h||) (where ||h|| is the Euclidean distance of vector h), its anisotropic version is: C_aniso(h) = C( ||A * h|| ), where the transformation matrix A is crucial. It maps the anisotropic distances in physical space to an isotropic "standard space" for calculating correlations through rotation and scaling operations.

[0057] Two main types include: 1. Geometric anisotropy, defined as follows: the correlation decays to different extents in different directions, but the variance at which it reaches the sill value (correlation decay is 0) is the same. Manifestation: Spatially, its isotropic profile is elliptical. Correlation decays slowly (larger range) along the major axis and rapidly (smaller range) along the minor axis. The role of matrix A: A is typically a rotation and scaling matrix. It first rotates the coordinates to the major axes of the ellipse, then scales them along different axes by different proportions (i.e., the reciprocal of the range parameter), and finally calculates the scaled Euclidean distance.

[0058] 2. Banded anisotropy: Definition: Different directions not only have different ranges, but also different sill values ​​(variance). Representation: It can be understood as the superposition of multiple structures with different directions, resulting in more complex isotropic contours. Modeling: It is usually constructed through a linear combination of multiple basic structures. Applications in specific models: The anisotropic covariance function can be embedded in almost all covariance-based Gaussian random field models, such as: 1. Kriging: In geostatistics, the anisotropic variogram (the counterpart of the covariance function) is directly used. 2. Gaussian process regression: In machine learning, it is achieved by modifying the distance metric of the kernel function (covariance function). 3. INLA-SPDE model: Anisotropy can be achieved by defining the coefficients of the differential operator in SPDE as varying with direction. During discretization, this is equivalent to using anisotropic basis functions on a triangular mesh or defining direction-dependent weights, so that the accuracy matrix of the generated GMRF reflects the direction-dependent structure.

[0059] The cross-covariance model provided in the specific embodiment of the present invention is based on the existing model, and simultaneously fits the spatial correlation of soil properties in the vertical and horizontal directions through the covariance function, that is, the spatial correlation has different decay rates or structures in different directions.

[0060] 3.2.2 Discretization using the SPDE method: The continuous Gaussian random field is discretized using the stochastic partial differential equation (SPDE) method. Discretize the constructed triangular mesh and transform it into a computationally efficient Gaussian Markov random field.

[0061] The SPDE (Stochastic Partial Differential Equation) method provided in this invention can be successfully discretized into a triangular mesh because its mathematical foundation is deeply compatible with the geometric properties of the triangular mesh. The solution of SPDE (such as a Gaussian field with Matérn covariance) is essentially a Gaussian Markov Random Field (GMRF), whose core characteristic is the local Markov property—that is, the value of any point depends only on the values ​​of its neighboring points and is independent of points further away. This locality corresponds precisely to the local connectivity structure of the triangular mesh: each node in the triangular mesh is directly related to only a few neighboring nodes (connected by edges), which perfectly matches the local dependency assumption of the Gaussian Markov Random Field.

[0062] From a numerical perspective, finite element analysis provides a robust mathematical tool for discretizing SPDEs. By dividing the continuous spatial domain into triangular elements and defining basis functions on each element, SPDEs can be transformed into a grid-node-based discrete system. This transformation not only preserves the spatial correlation of the original continuous field but also generates a sparse precision matrix, significantly improving computational efficiency. The flexibility of triangular meshes allows them to adapt to various complex spatial regions and boundary shapes, whether regular geographical areas or irregular three-dimensional volumes; accurate discretization can be achieved by adjusting the mesh density.

[0063] In practical applications, such as the three-dimensional spatial prediction of soil organic carbon in Poyang Lake wetland in this study, the combination of the SPDE method and triangular meshes demonstrated its advantages. By constructing an appropriate three-dimensional triangular mesh and setting reasonable mesh parameters (such as cutoff distance and maximum side length), the SPDE method can transform the continuous soil organic carbon field into a discrete representation based on mesh nodes. This discretization not only accurately captures the spatial variation characteristics of soil properties in the horizontal and vertical directions, but also significantly reduces computational complexity through sparse matrix operations, making efficient processing of large-scale three-dimensional spatial data possible.

[0064] 3.2.3 INLA Algorithm Inference: Bayesian inference is performed using the INLA algorithm. INLA calculates the model parameters (fixed effects) through analytical approximation and numerical integration. Random fields The posterior distribution of hyperparameters (such as variance, spatial correlation range, etc.) avoids the problems of large computational cost and difficult convergence in traditional MCMC methods.

[0065] 3.2.4 Software and Computation: In the R v4.3.3 environment, call... INLA The package performs model fitting. The model runs on a computer configured with an Intel i7-8750H CPU (quad-core @ 2.20 GHz), and a complete model training and inference cycle takes approximately 50 minutes.

[0066] 3.3 Three-dimensional spatial prediction and uncertainty quantification: 3.3.1 SOC Content Prediction: Using the trained model, the SOC content in the three-dimensional space (90-meter horizontal resolution, 10-centimeter vertical interval) of the entire study area is predicted, generating a series of three-dimensional spatial distribution maps from the surface layer (0~10 cm) to the bottom layer (90~100 cm) (see...). Figures 5-8 The 3D-INLA-SPDE model comprehensively captures the vertical decay and horizontal heterogeneity of SOC in Poyang Lake wetlands. Vertical SOC decay refers to the systematic exponential or logarithmic decrease in SOC content with increasing soil depth. This is determined by the regularity of organic matter input, decomposition, and migration processes along the soil profile. In the model, this is typically treated as a fixed effect. Horizontal SOC heterogeneity refers to the patchy, continuously varying variability of SOC content in horizontal space at the same or similar depths. This is driven by spatial heterogeneity factors such as topography, vegetation, microclimate, soil texture, and human activities. In the model, this is typically captured by spatial random effects. SOC levels decreased significantly in the 0–10 cm soil layer (0.34–8.18 g / kg) and decreased to 0.19–4.12 g / kg below 30 cm. High organic carbon zones (>6.0 g / kg) were persistently accumulated in the western, central, and eastern parts of the study area. However, their size decreased with increasing soil depth. In contrast, the low-carbon organic carbon zone (< 4.0 g / kg) expands in the periphery and southern regions. This spatial distribution pattern is consistent with the predictions of the 3D-EBK, 3D-IDW, and 3D-RF models.

[0067] 3.3.2 Uncertainty Mapping: Based on the Bayesian posterior prediction distribution, the prediction interval for each location in the study area at a 95% confidence level is calculated to quantify the uncertainty of the three-dimensional spatial prediction results. For example... Figure 9-11As shown, prediction uncertainty exhibits a regular spatial differentiation pattern within the study area. Spatially, high uncertainty areas are mainly distributed at the edges of the study area. This is primarily attributed to two factors: first, sample points at the regional edges are usually sparse, lacking effective constraints from neighboring observation data during spatial interpolation or extrapolation; second, the spatial covariance structure upon which the model relies has inherent limitations at the computational boundaries, failing to fully utilize symmetric neighborhood information, thus increasing prediction variance. Furthermore, prediction uncertainty shows significant vertical differentiation with soil depth. In the topsoil layer (0–30 cm), uncertainty is relatively high. This is partly due to the inherently high and highly variable SOC content in this layer, and also because the topsoil is more directly affected by surface dynamic environmental factors and human activities, resulting in more complex spatial processes. In contrast, while deeper soil layers (>30 cm) face greater challenges due to data sparsity, their lower organic carbon content and less direct surface disturbance result in relatively lower uncertainty at this depth.

[0068] In summary, the uncertainties in this embodiment mainly stem from two aspects: first, data limitations, particularly the difficulty in deep soil sampling leading to sparse data that cannot fully represent complex deep soil conditions; second, the necessary simplification of the model structure, as the soil system involves complex interactions among multiple factors such as climate, vegetation, topography, soil type, and human activities, and the model cannot fully incorporate all processes, which also limits the prediction accuracy to some extent. This technical solution, through the explicit quantification of three-dimensional uncertainty, provides crucial information for assessing the reliability of prediction results, identifying key data collection areas, and supporting risk-based management decisions.

[0069] (4) Model accuracy verification and comparative analysis: To evaluate the performance of the method of the present invention, the following verifications were performed: 4.1 Accuracy Validation Metrics: An evaluation was conducted using a retained validation set (250 stratified soil profile samples). Two metrics were selected: the coefficient of determination (R²). 2 ) and root mean square error (RMSE).

[0070] 4.2 Validation Results: Model accuracy validation results ( Figure 13 Figures (a)-(d) show the differences in performance from four subgraphs (3D-IDW, 3D-EBK, 3D-RF, 3D-INLA-SPDE). The 3D-INLA-SPDE model ( Figure 13 The middle (d) prediction accuracy is the highest, R 2 The R² value of 0.80 indicates a strong consistency between the predicted and measured values. Furthermore, the model has the lowest RMSE (1.04 g / kg), indicating the smallest prediction error. Conversely, the R² values ​​of the 3D-IDW, 3D-EBK, and 3D-RF models are significantly lower. 2The lower values ​​(0.37, 0.45, 0.40) and higher values ​​(1.85, 1.74, 1.88) indicate poor predictive ability. Overall, the 3D-INLA-SPDE model has high interpolation accuracy and reliability, making it the best choice for three-dimensional soil organic carbon prediction mapping tasks.

[0071] (5) Analysis of environmental driving factors and estimation of soil carbon storage 5.1 Driver factor analysis based on SHAP values: 5.1.1 Implementation Method: Using the R language's " kernelshap "Based on a trained 3D-INLA-SPDE model and training dataset, this study calculates the SHAP value of each environmental covariate, and variable importance ranking and partial dependency graph analysis elucidate the influence mechanisms of soil properties, climate, topography, and biological factors on organic carbon dynamics."

[0072] 5.1.2 Global Importance: The results show that Sand is the most important factor for soil organic carbon uptake, followed by AK, Silt, pH, etc. Soil properties contribute the most to the mapping of soil organic carbon (53.7%). Figure 14 The most significant contribution was (a) of the population, followed by topography (23.0%), climate (12.1%), and biota (11.2%). This was based on the global contribution analysis and importance ranking chart of SHAP. Figure 14 In (b) of the model, the importance of different environmental covariate categories to the three-dimensional soil organic carbon changes in the 3D-INLA-SPDE model is ranked as follows: soil properties (0.594) > topography (0.255) > climate (0.134) > biota (0.124). Specifically, the top ten variables with the highest absolute values ​​of SHAP are: Sand (0.085) > AK (0.066) > Slit (0.064) > pH (0.059) > TK (0.057) > AN (0.057) > NDVI (0.046) > Slope (0.046) > EVI (0.046) > SP (0.042).

[0073] 5.1.3 Direction and extent of influence: Draw a partial dependency graph (see...) Figures 15-18 As shown in the figure, the influence of individual variables on SOC prediction under different values ​​is analyzed to determine the direction and extent of their impact. For example, the influence of soil property variables on SOC accumulation is illustrated: the partial dependence plot shows that when Sand exceeds 25%, its impact on soil organic carbon changes from positive to negative. The increase in Sand initially disrupts the binding of organic matter with the soil, weakening its positive impact on organic carbon, but later forms a stable microenvironment, promoting SOC accumulation. AK content is significantly positively correlated with SHAP value (…). p< 0.001). AK indirectly enhances the input and stabilization of organic matter by promoting plant root growth and microbial activity. Soil pH bias-dependent plots show that the positive effect of pH on soil organic carbon is greatest between 5.4 and 5.6. Beyond this range, the microbial community structure becomes unbalanced, and the decomposition rate increases.

[0074] Topographic factors: The partial dependence plot of TWI shows that in low-lying areas with a TWI of 110, TWI has a positive impact on SOC. This is because long-term flooding inhibits the activity of aerobic microorganisms, reducing the decomposition of organic carbon. The partial dependence plot of MrVBF shows that when MrVBF exceeds 8.0, the area is typically a depression or valley floor with very flat terrain and a large catchment area. Such topographic conditions usually lead to long-term or seasonal water accumulation in the soil, forming an anaerobic environment. This slows down the decomposition rate of organic matter. At the same time, this stable topographic unit is also more likely to capture and deposit organic matter, thus jointly promoting the accumulation of soil organic carbon.

[0075] Climate factors: The partial dependence plot of MAP shows that as MAP increases, its impact on soil organic carbon (SOC) increases positively. This is mainly due to increased precipitation improving plant productivity and increasing litter return. The partial dependence plot of MAT shows that as MAT increases, the impact of SHAP on SOC gradually becomes negative, but the overall impact remains relatively small. This is because the study area is relatively small, and the overall variation range of MAT is less than 1°C. However, as temperature rises, microbial activity increases, thereby accelerating the decomposition of organic matter and increasing SOC loss.

[0076] Biological factors: The partial dependence plot of NDVI and EVI shows that in Poyang Lake wetland, the positive effects of NDVI and EVI on organic carbon gradually weaken with the increase of vegetation cover. When NDVI>0.7 and EVI>0.5, the positive effects of NDVI and EVI on organic carbon are slightly weakened.

[0077] 5.2 Soil Organic Carbon Storage (SOCS) Estimation: Based on the obtained three-dimensional SOC content distribution map, combined with the estimated values ​​from BD and gravel content data, the soil organic carbon density (SOCD, kg / m³) of each grid cell was calculated. 2 ).

[0078] The soil organic carbon concentration (SOCS) of the 0–100 cm soil layer in Poyang Lake wetland was obtained by integrating (summing) the SOCD of all grid cells in the three-dimensional space of the entire study area. The estimation results of this embodiment are shown in Table 2. Although the estimated values ​​of SOCS differ slightly (92.24 Tg C for 3D-INLA-SPDE, while 107.80–114.83 Tg C for other models), these differences reflect different methods for handling spatial autocorrelation and data distribution. Despite the differences in SOCS quantification, all models consistently emphasize the advantage of surface organic carbon accumulation and the fragmented distribution of deep high-organic carbon patches. The 3D-INLA-SPDE method demonstrates significant robustness in capturing vertical trends and spatial clustering, providing important insights for wetland carbon management and soil conservation strategies.

[0079] In summary, this invention provides a complete, detailed, and highly operable three-dimensional Bayesian spatial mapping and driven analysis scheme for soil properties. Those skilled in the art can implement this invention by appropriately adjusting the parameters based on the above description and the specific research area and data conditions. The above embodiments are merely illustrative of the technical solutions of this invention and are not intended to limit the scope of protection of this invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this invention should be included within the scope of protection of this invention.

[0080] Table 1. Soil properties and environmental covariate data Table 2 Soil Carbon Storage (SOCS) Results

Claims

1. A method for three-dimensional mapping of soil properties based on integrated nested Bayesian and Laplace approximation, characterized in that, include: Multiple soil profile samples from the study area are obtained, along with the three-dimensional spatial coordinates, soil property measurements, and environmental covariates for each soil profile sample, thereby constructing a sample set. The sample set is then divided into a training sample set and a validation sample set. The three-dimensional spatial coordinates include longitude, latitude, and depth. The 3D-INLA-SPDE model is constructed as follows: based on the three-dimensional spatial coordinates of the study area boundary and multiple soil profile samples, the constrained Delaunay triangulation method is used to obtain a three-dimensional unstructured mesh covering the entire study area; the SPDE method is used to discretize the three-dimensional Gaussian random field on the three-dimensional unstructured network to obtain a Gaussian Markov random field; based on the Gaussian Markov random field, the training samples are used as input, and Bayesian inference is performed based on the three-dimensional Bayesian model using the INLA algorithm to calculate the posterior distribution parameters of the model parameters. The 3D-INLA-SPDE model is trained using a training sample set. Based on the environmental covariates and three-dimensional spatial coordinates of the area under study, the model is trained to obtain the predicted values ​​of soil properties in three-dimensional space of the area under study, thereby realizing three-dimensional mapping of soil properties.

2. The method for three-dimensional mapping of soil properties based on integrated nested Bayesian and Laplace approximation as described in claim 1, characterized in that, In a three-dimensional unstructured mesh covering the entire study volume using the constrained Delaunay triangulation method, the fineness of the three-dimensional unstructured mesh is controlled by key parameters, including offset parameters, truncation parameters, and longest edge length parameters. The offset parameters are used to define the mesh boundary extension distance, the truncation parameters are used to set the minimum allowable distance between vertices, and the longest edge parameter is used to set the maximum allowable boundary of the triangle.

3. The method for three-dimensional mapping of soil properties based on integrated nested Bayesian and Laplace approximation as described in claim 1, characterized in that, The three-dimensional Bayesian model is: ,in, These are soil property observation values; This is an environmental covariate matrix constructed from environmental covariates; This is the fixed effects coefficient; To characterize a three-dimensional Gaussian random field with spatial dependence, For Gaussian measurement error, the posterior distribution parameters include The variance of S and the spatial correlation range.

4. The method for three-dimensional mapping of soil properties based on integrated nested Bayesian and Laplace approximation as described in claim 1, characterized in that, After obtaining the environmental covariates, a recursive feature elimination method is used to select features for the environmental covariates.

5. The method for three-dimensional mapping of soil properties based on integrated nested Bayesian and Laplace approximation according to claim 4, characterized in that, A recursive feature elimination method is used for feature selection of environmental covariates, including: A k-fold cross-validation strategy is adopted, with RMSE and other indicators as evaluation criteria. The optimal subset of feature variables with the smallest prediction error for the target soil properties is iteratively selected from the obtained set of environmental covariates. Then, the soil property measurement values ​​and the corresponding values ​​of the optimal feature variable subset at the corresponding spatial locations are associated with the nodes of the constructed three-dimensional unstructured grid to form the structured dataset required for model training, i.e., the sample set.

6. The method for three-dimensional mapping of soil properties based on integrated nested Bayesian and Laplace approximation as described in claim 1, characterized in that, The SPDE method is used to discretize the three-dimensional Gaussian random field onto the three-dimensional unstructured network to obtain a Gaussian Markov random field, including: Based on the parameters of a three-dimensional unstructured mesh (i.e., the mesh cells do not have a fixed arrangement rule or topological structure), a SPDE model is constructed by setting the Matérn covariance function, and the SPDE model is solved by the finite element method to obtain a discretized Gaussian Markov random field.

7. The method for three-dimensional mapping of soil properties based on integrated nested Bayesian and Laplace approximation as described in claim 1, characterized in that, Methods for constructing three-dimensional Gaussian random fields include: A Gaussian random field is constructed in the horizontal dimension using the stochastic partial differential equation method. The spatial correlation in the horizontal direction is obtained by defining an exponential model or a Gaussian model based on the constructed Gaussian random field in the horizontal dimension. The horizontal correlation and vertical correlation are separated and modeled using an anisotropic covariance function. The fluctuation range constraint or correlation decay rule in the vertical direction is established by using a cross covariance model to capture the indirect influence of deep soil on shallow parameters. The final model can simultaneously characterize the horizontal spatial variation and vertical distribution characteristics.

8. The method for three-dimensional mapping of soil properties based on integrated nested Bayesian and Laplace approximation according to claim 1, characterized in that, When the soil property is soil organic carbon content, soil bulk density and gravel content data are calculated based on the predicted soil property values ​​in three-dimensional space and empirical formulas. The soil organic carbon density and soil organic carbon storage in the study area are then calculated using these formulas. The empirical formula for soil bulk density (BD) is as follows: （1） The soil organic carbon density (SOCD) of the i-th soil sample is: （2） The soil organic carbon (SOCS) storage of the i-th soil sample is: （3） In the formula, Let be the sampling depth of the i-th soil sample. This represents the side length of a cell in a three-dimensional unstructured mesh. Let be the gravel content of the i-th soil sample.

9. A three-dimensional soil property mapping device based on integrated nested Bayesian and Laplace approximation, characterized in that, include: One or more processors; Memory, used to store one or more programs; When the one or more programs are executed by the one or more processors, the one or more processors implement the three-dimensional mapping method for soil properties based on integrated nested Bayesian and Laplace approximation as described in any one of claims 1-9.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by the processor, the program implements the three-dimensional mapping method for soil properties based on integrated nested Bayesian and Laplace approximation as described in any one of claims 1-9.

Images