Soil thickness type prediction method based on feature ensemble learning

Abstract

The application discloses a soil thickness type prediction method based on feature ensemble learning and relates to the technical field of machine learning and digital soil mapping. The method comprises the following steps: S1, based on the measured thickness data of a target region, calibrating a parameterized geomorphology-soil process model, generating simulated thickness data, and fusing the simulated thickness data with the measured thickness data to generate a virtual-real fusion training set; and S2, extracting basic environmental variables of the target region and constructing proxy features of quantitative key soil-forming processes based on a pedogenesis function. The application performs constraint causal analysis for feature screening, embeds domain knowledge in the form of hard constraints into the causal discovery algorithm, eliminates false correlation features, and obtains an integrated feature set with clear causality and strong interpretability, thereby improving the physical reliability and extrapolation robustness of the model prediction results.

Description

Technical Field

[0001] This invention relates to the field of interdisciplinary technology of machine learning and digital soil mapping, specifically a method for predicting soil thickness type based on feature ensemble learning. Background Technology

[0002] Soil thickness is a key attribute determining soil moisture retention, nutrient reserves, vegetation productivity, and suitability for engineering construction. Obtaining high-precision, large-scale spatial distribution information on soil thickness types is crucial for agricultural production, ecological restoration, soil and water conservation, and land use planning. Currently, the "soil-landscape model," based on statistical learning using environmental variables and measured soil thickness, is the mainstream technical approach. However, it suffers from limitations such as limited data sources, low reliability of feature correlations, and static, rigid models.

[0003] For example, the invention patent with patent publication number CN114970934A discloses a soil thickness type prediction method based on feature integration learning. It effectively utilizes environmental variables that affect soil thickness changes to construct a prediction model. By screening the optimal set of environmental variables, it constructs a continuous soil thickness integrated prediction model and a soil depth interval integrated prediction model, and finally obtains the spatial distribution of soil thickness intervals covering the target area.

[0004] However, the above and similar technical solutions still have the following shortcomings: feature selection relies on statistical correlation, which may result in the selected optimal environmental variables containing spurious correlation features, leading to poor reliability and weak physical interpretability when the model extrapolates and predicts; the prediction model is a static mapping and does not encode key physical processes such as soil erosion-deposition, making it difficult for the model to simulate the dynamic evolution of soil thickness under the dynamic influence of multiple factors such as climate and human activities; model ensemble remains at the decision-level weighted average, making it difficult to capture spatial physical correlations and continuous temporal evolution patterns. Summary of the Invention

[0005] The purpose of this invention is to provide a soil thickness type prediction method based on feature ensemble learning to solve the problems mentioned in the background art.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a soil thickness type prediction method based on feature ensemble learning, comprising:

[0007] S1. Based on the measured thickness data of the target area, calibrate the parameterized geomorphic-soil process model, generate simulated thickness data, and fuse it with the measured thickness data to generate a virtual-real fusion training set;

[0008] S2. Extract the basic environmental variables of the target area and construct surrogate features of key soil-forming processes based on soil genesis functions; use constrained causal analysis to select a subset of features that are causally related to soil thickness from the basic environmental variables and surrogate features to form an integrated feature set.

[0009] S3. Using spatial grids as nodes, construct a spatiotemporal graph skeleton with physically meaningful edge weights based on hydrological flow direction; input the integrated feature set of multiple time phases as node attributes into the spatiotemporal fusion network guided by physical mechanisms, and train the network using the virtual-real fusion training set as supervision labels; the network achieves prediction through iterative spatial aggregation and temporal evolution, and is optimized using a gradient-aware adaptive multi-task loss function.

[0010] S4. The integrated feature set at the current moment is used as a node attribute and loaded into the spatiotemporal graph skeleton to form the current spatiotemporal graph. It is then input into the trained spatiotemporal fusion network and, combined with the preset soil thickness classification standard, a soil thickness type prediction result is generated. The time-series integrated feature set under the future scenario is used as a node attribute sequence to construct a future spatiotemporal graph sequence. This sequence is then input into the spatiotemporal fusion network for multi-step rolling prediction to generate a prediction result of the dynamic evolution trend of soil thickness type and a soil degradation risk warning information.

[0011] Furthermore, S1 specifically includes:

[0012] The soil thickness data of the measured profile points were obtained and cleaned using an anomaly detection method combined with spatial autocorrelation.

[0013] Based on the dominant soil-forming process in the target area, a parameterized geomorphic-soil process model is selected, and prior probability distributions are set for its key parameters.

[0014] The target area is divided into geomorphic units. For each geomorphic unit, the posterior probability distribution of the model parameters is solved by variational inference algorithm using measured thickness data within the unit, thus completing the probabilistic calibration of the model.

[0015] Based on the posterior distribution of the calibrated parameters, simulated thickness data and its uncertainty are generated; a confidence evaluation neural network is constructed and trained, and the confidence weight of the simulated point is output according to the local environmental characteristics, simulation uncertainty and nearby measured residuals.

[0016] The measured thickness data and the weighted simulated data are combined to form a virtual-real fusion training set.

[0017] Furthermore, S2 specifically includes:

[0018] The basic environmental variables of the target area are extracted to form a basic environmental feature layer; based on the preset soil genesis function, the proxy features that quantify the intensity of key soil-forming processes are calculated from the basic environmental variables to form a process proxy feature layer.

[0019] At the measured soil thickness points, data of the corresponding basic feature layer and process proxy feature layer are collected; a domain knowledge constraint matrix is ​​constructed, and a causal relationship network between variables is generated by a constraint causal discovery algorithm based on the matrix; when searching for the causal relationship network structure, the constraint causal discovery algorithm treats the constraints in the matrix as inviolable hard conditions.

[0020] From the causal relationship network, features that have direct or indirect causal paths with soil thickness are identified and formed into feature subsets; based on the feature subsets, feature data are extracted and integrated on the whole-area grid to form an integrated feature set.

[0021] Furthermore, the basic environmental variables include climate variables, parent material type variables, and slope variables; the proxy features include the erosion and transport capacity index; the constraints of the domain knowledge constraint matrix in S22 include prohibiting the existence of direct causal edges between climate variables and parent material type variables, and forcing the slope variable to influence the soil thickness variable through the erosion and transport capacity index.

[0022] Furthermore, in step S3, a spatiotemporal graph skeleton with physically meaningful edge weights is constructed based on the hydrological flow direction. Specifically, this includes: calculating the flow direction of each grid node based on the digital elevation model; for any grid node, establishing a directed edge from itself to its downstream adjacent grid node; the weight of the directed edge is calculated based on the elevation difference and horizontal distance between upstream and downstream nodes.

[0023] Furthermore, the spatiotemporal fusion network guided by the physical mechanism in step S3 specifically includes:

[0024] The spatial aggregation layer employs a physically constrained graph attention mechanism, where the calculation of attention weights is based on the correlation of node feature learning and the physical weights pre-calculated based on hydrological connectivity.

[0025] The time evolution layer employs a gated recurrent unit network to independently process the state sequence of each node and output the predicted state for future time steps.

[0026] The iterative spatiotemporal fusion module inputs the output of the time evolution layer back into the spatial aggregation layer to perform at least one round of time prediction-spatial correction loop to generate the final prediction result.

[0027] Furthermore, the construction and optimization method of the gradient-aware adaptive multi-task loss function in S3 includes: the loss function is composed of a weighted sum of a state prediction loss term and a change direction prediction loss term; the weight coefficients of the loss function are dynamically adjusted inversely during training based on the gradient magnitude information generated by each loss term on the spatiotemporal fusion network parameters.

[0028] Furthermore, the generation method for the dynamic evolution trend results of soil thickness type and early warning information of soil degradation risk in S4 includes:

[0029] Acquire temporal environment data under future scenarios, process it in step S2 to obtain temporal integrated feature sets for each future time step, and construct an initial future spatiotemporal map;

[0030] The future spatiotemporal map of the current time step is input into the spatiotemporal fusion network to predict the soil thickness value of the next time step; the predicted soil thickness value is used as a feedback feature to update the attributes of the corresponding node in the next time step.

[0031] By combining updated node attributes and future scenario data, a new future spatiotemporal map is constructed as the input for the next round of prediction; the above iterative steps are repeated until the simulation of the preset duration is completed.

[0032] The soil thickness type change sequence of each spatial grid within the statistical extrapolation period is used to generate a soil thickness type transition matrix and a dynamic evolution diagram of soil thickness type at key time nodes.

[0033] Based on the simulation results, the annual average rate of change of soil thickness for each grid is calculated; areas with a rate of change exceeding a preset degradation threshold are identified as risk degradation zones, and soil degradation risk warning information is generated.

[0034] Furthermore, the method also includes model decision analysis, which uses Shapley additive interpretation based on the trained spatiotemporal fusion network to analyze the contribution of each feature in the integrated feature set to the prediction result.

[0035] Compared with the prior art, the beneficial effects of the present invention are:

[0036] The soil thickness type prediction method based on feature ensemble learning constructs a virtual-real fusion training set, which integrates simulated data generated by the physical model with measured data through credibility weighting, significantly expanding the high-quality training samples and improving the generalization ability of the prediction model in data-sparse areas and complex terrain areas. The Shapley additive interpretation method is used to analyze the model decision, quantifying the contribution of each environmental feature to the prediction results, enhancing the interpretability of the prediction results and the credibility of decision support.

[0037] Meanwhile, by performing constrained causal analysis for feature selection, domain knowledge is embedded into the causal discovery algorithm in the form of hard constraints, eliminating false association features and obtaining an integrated feature set with clear causality and strong interpretability, thus improving the physical reliability and extrapolation robustness of the model's prediction results. By constructing a spatiotemporal fusion network guided by physical mechanisms, using hydrological flow direction to build a graph skeleton, and fusing graph attention and gated recurrent units for iterative prediction, the physical laws of surface material transport are encoded into the network learning, realizing the transformation from static statistical mapping to dynamic process simulation. By adopting a gradient-aware adaptive multi-task loss function to optimize the prediction model, the weights are dynamically adjusted according to the gradient of the sub-task loss, automatically balancing the optimization process of state prediction and trend prediction of change direction, resulting in more stable and accurate dynamic prediction results, providing reliable support for early warning of soil degradation risks. Attached Figure Description

[0038] Figure 1 This is a schematic diagram of the soil thickness type prediction method of the present invention;

[0039] Figure 2 This is a schematic diagram of the virtual-real fusion training set generation method of the present invention;

[0040] Figure 3 This is a schematic diagram of the binding causal discovery and feature selection method of the present invention;

[0041] Figure 4 This is a schematic diagram of the spatiotemporal fusion network architecture and training method guided by the physical mechanism of the present invention. Detailed Implementation

[0042] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0043] like Figure 1 As shown, this invention provides a technical solution: a soil thickness type prediction method based on feature ensemble learning, comprising:

[0044] S1. Based on the measured thickness data of the target area, calibrate the parametric geomorphic-soil process model, generate simulated thickness data, and fuse it with the measured thickness data to generate a virtual-real fusion training set.

[0045] like Figure 2 As shown, the present invention provides a method for generating a virtual-real fusion training set;

[0046] Specifically:

[0047] The soil thickness data of the measured profile points were obtained and cleaned using an anomaly detection method that combines spatial autocorrelation.

[0048] It is important to note that the dataset of measured points containing geographic coordinates and soil thickness values ​​should be obtained from soil surveys or related research projects in the target area. To overcome the neglect of spatial autocorrelation by conventional statistical cleaning methods, this step introduces a graph-based intelligent cleaning method. Specifically, a spatial adjacency graph is constructed based on the geographic proximity of all measured points (e.g., using the K-nearest neighbor method). Then, an unsupervised learning process is performed on this adjacency graph structure using a graph autoencoder; this encoder consists of a two-layer graph convolutional network used to learn low-dimensional embedding representations of nodes, while the decoder reconstructs the original adjacency relationships from the embeddings. After training, the adjacency relationship reconstruction error for each node is calculated. Nodes with reconstruction errors significantly higher than the average level indicate anomalies in the spatial association pattern between their thickness values ​​and surrounding neighbors. Combining basic geographical knowledge (e.g., checking whether the point is located on an extremely steep slope but reports a thick soil layer), these nodes are identified as spatial context anomalies and removed, ultimately yielding a reliable dataset of measured soil thickness.

[0049] Based on the dominant soil-forming processes in the target area, a parameterized geomorphic-soil process model is selected, and prior probability distributions are set for its key parameters.

[0050] It is important to note that a parametric geomorphic-soil process model (such as the SCAPE model) should be selected based on the dominant soil-forming process in the target area (e.g., slope erosion-deposition process). To quantify the uncertainty of the model, key parameters in the model (e.g., soil productivity, material diffusion coefficient) are represented as probability distributions. In this embodiment, they are defined as log-normal distributions to ensure that the parameters are positive. The hyperparameters of the distribution (e.g., mean μ and standard deviation σ) are set based on relevant geological, climatological, and soil science literature as prior knowledge of the parameters. For example, in temperate hilly areas, the prior mean of log(soil productivity) can be set to log(0.05), and the standard deviation can be set to 0.5.

[0051] The target area is divided into geomorphic units. For each geomorphic unit, the posterior probability distribution of the model parameters is solved using the measured thickness data within the unit through variational inference algorithm, thus completing the probabilistic calibration of the model.

[0052] It is important to note that terrain classification algorithms (such as DEM-based clustering algorithms) are used to divide the target area into multiple relatively homogeneous geomorphic units. For each geomorphic unit, a variational inference method is employed for probabilistic calibration, i.e., a variational distribution (usually assumed to be Gaussian) parameterized by a lightweight neural network is constructed to approximate the true posterior distribution. The input of the neural network is the aggregated terrain features of the geomorphic unit, and the output is the parameters of the variational distribution, which are optimized by maximizing the lower bound of evidence. The variational network parameters of all units are optimized using gradient descent, ultimately resulting in a series of inference networks that can output parameter probability distributions based on the unit's terrain features, thus completing the localization of model parameters and the quantification of uncertainty.

[0053] Based on the posterior distribution of the calibrated parameters, simulated thickness data and its uncertainty are generated; a confidence evaluation neural network is constructed and trained, and the confidence weight of the simulated point is output according to the local environmental characteristics, simulation uncertainty and nearby measured residuals.

[0054] It is important to note that, using the calibrated model, multiple sets of parameters are sampled from the posterior distribution of each landform unit, and the landform-soil process model is run to generate multiple sets of simulated thickness data for the entire region. The simulation results are aggregated, and the mean and variance of the simulated thickness on each spatial grid are calculated. For each simulated data point, a confidence weight is dynamically assigned through a confidence assessment neural network. The input features of this network are a comprehensive vector, including: environmental features of the point (such as slope, lithology), physical model uncertainty features (variance), local fitting features (residual statistics of simulated mean and neighboring measured values), and data density features. The network output is a weight value in the interval [0,1]. The training of the confidence assessment neural network adopts a meta-learning framework: in the inner loop, the simulated data and a portion of the measured data (meta-training set) are fused using the weights generated by the current network, and a surrogate prediction model (such as MLP) is quickly trained; in the outer loop, the prediction error of the surrogate prediction model on another portion of the measured data (meta-validation set) is used as the meta-loss, and the parameters of the confidence assessment neural network are updated through gradient backpropagation.

[0055] The measured thickness data and the weighted simulated data are combined to form a virtual-real fusion training set.

[0056] It is important to note that the cleaned measured thickness dataset (each point is assigned a fixed high weight of 1.0) is merged with the simulated thickness dataset with dynamic weights assigned by the aforementioned meta-learning credibility evaluation neural network to form a virtual-real fusion training set for subsequent spatiotemporal fusion network training.

[0057] S2. Extract the basic environmental variables of the target area and construct proxy features for quantifying key soil-forming processes based on soil genesis functions; use constrained causal analysis to select a subset of features that are causally related to soil thickness from the basic environmental variables and proxy features to form an integrated feature set.

[0058] like Figure 3 As shown, the present invention provides a constrained causal discovery and feature selection method;

[0059] Specifically:

[0060] The basic environmental variables of the target area are extracted to form a basic environmental feature layer; based on the preset soil genesis function, the proxy features that quantify the intensity of key soil formation processes are calculated from the basic environmental variables to form a process proxy feature layer.

[0061] It is important to note that, firstly, basic environmental variables are extracted using public databases and local geographic information systems on a unified geographic grid of the target area. These basic environmental variables include, but are not limited to: climate variables, such as annual mean temperature (MAT) (°C) and annual precipitation (MAP) (mm); parent material type, represented by unique thermal encoding based on digital geological maps or converted into continuous variables of potential element content; and topographic variables, such as slope, planar curvature, and topographic humidity index. Then, based on soil genesis principles, a series of functions are pre-defined to map these basic environmental variables into proxy features characterizing the intensity of key soil-forming processes. For example, the weathering potential index, used to comprehensively characterize the influence of temperature and moisture on the weathering of primary minerals, is calculated using the following formula:

[0062]

[0063] Where k is an adjustment constant, such as 10; the erosion transport capacity index is used to quantify the potential of slope runoff to strip and transport soil particles, and the calculation formula is:

[0064]

[0065] The deposition probability index, used to indicate the depositional environment, can be directly proxies using the topographic humidity index. All calculations are performed independently on each spatial grid, generating raster data with the same resolution as the underlying environmental feature layer.

[0066] At the measured soil thickness points, data of the corresponding basic feature layer and process proxy feature layer are collected; a domain knowledge constraint matrix is ​​constructed, and a causal relationship network between variables is generated by a constraint causal discovery algorithm based on the matrix; when searching for the causal relationship network structure, the constraint causal discovery algorithm treats the constraints in the matrix as inviolable hard conditions.

[0067] It is important to note that, at the spatial location corresponding to the reliable measured thickness dataset, the values ​​of all basic environmental variables and process proxy features are collected to form a multivariate observation data matrix, the dimension of which is the number of samples multiplied by the number of features.

[0068] Construct a square constraint matrix corresponding to the number of features. The rows and columns of this matrix correspond one-to-one with all candidate features in a defined order. The value of each element in the matrix represents the type of prior constraint on the causal relationship between a pair of features: when the element value is -1, the corresponding row feature is prohibited from being the direct cause of the column feature; when the element value is 1, the corresponding row feature is required to be the direct cause of the column feature; when the element value is 0, no prior constraint is imposed on the causal relationship between the two features, and discovery is entirely data-driven.

[0069] Based on the fundamental principles of soil genesis, this matrix includes at least the following core constraints: Elements corresponding to the parent material type variable for climate variables are set to -1 to encode the domain knowledge that climate cannot alter the geological parent material; elements corresponding to the erosion transport capacity index for slope variables and elements corresponding to the erosion transport capacity index for soil thickness are both set to 1 to force the encoding that slope affects soil thickness by influencing erosion capacity. A causal network structure learning algorithm based on scoring function optimization is employed, such as a greedy search algorithm using the Bayesian information criterion as the scoring standard. During algorithm execution, the domain knowledge constraint matrix must be integrated into its search logic as an inviolable hard condition. Specifically, during algorithm initialization, corresponding directed edges must be added to the initial network structure at all constraint matrix elements that are 1. During the iterative search process, any modifications to the network structure must first be verified to ensure they do not conflict with the constraints. Any modification that adds an edge at a constraint of -1 or deletes an edge at a constraint of 1 will be rejected. Modifications are only accepted if they involve feature pairs with zero constraints or meet mandatory constraints, and are ultimately decided based on their improvement to the overall network score.

[0070] In the obtained causal relationship network, soil thickness is used as the target node, and a depth-first search algorithm is used to traverse all directed paths pointing to this node. All feature variables that affect soil thickness through direct or indirect paths are identified, forming a candidate feature subset. The data of these features across all grids in the entire region are extracted and integrated to form an integrated feature set, which serves as the input to the subsequent spatiotemporal fusion network.

[0071] S3. Using spatial grids as nodes, construct a spatiotemporal graph skeleton with physically meaningful edge weights based on hydrological flow direction; input the integrated feature set of multiple time phases as node attributes into the spatiotemporal fusion network guided by physical mechanisms, and train the network using the virtual-real fusion training set as supervision labels; the network achieves prediction through iterative spatial aggregation and temporal evolution, and is optimized using a gradient-aware adaptive multi-task loss function.

[0072] It is important to note that the spatial grid is a regular square grid system covering the target area, and its resolution is determined based on the resolution of the input digital elevation model data and application requirements (e.g., 30m × 30m). All environmental feature data are preprocessed to the exact same spatial reference and resolution as this grid system. Based on the digital elevation model aligned with the spatial grid, a multi-flow direction algorithm is used to calculate the flow direction of each grid node. For each grid node, a directed edge is established from that node to its directly downstream adjacent grid node according to its flow direction. Grids outside the target area or parts within the area that do not need to be considered (such as large bodies of water or urban built-up areas) can be defined as masks, and these nodes are excluded when constructing the spatiotemporal graph skeleton. A physical weight W is assigned to each directed edge (a, b). ab This weight is used to quantify the potential material and information transmission intensity from upstream nodes to downstream nodes. The weight is calculated based on hydrodynamic principles, and the formula is as follows:

[0073]

[0074] Where, p ab It is the probability allocation value from upstream node a to downstream node b in the multi-flow algorithm. The elevation difference between upstream node a and downstream node b. Let be the horizontal distance from upstream node a to downstream node b, and let β be an adjustable scale factor (e.g., β=1.0). This formula makes the connection weight between units with large elevation differences and short distances higher, and vice versa, thus accurately reflecting the physical law of the change of surface material transport efficiency with hydraulic gradient.

[0075] Thus, a spatiotemporal graph skeleton with physically meaningful weights is constructed, where nodes represent spatial locations, directed edges represent the main paths of matter transport, and the weights on the edges represent the transmission efficiency of the paths.

[0076] The input to the spatiotemporal fusion network is the node feature vector of each node at each time step, which originates from the ensemble feature set generated in step S2. Specifically, the ensemble feature set is a multi-dimensional data cube with dimensions of (number of time steps, number of grid rows, number of grid columns, and number of features). For a given time step and spatial location, the values ​​of all features at that location are extracted from the ensemble feature set to form the feature vector of that node at that time step.

[0077] like Figure 4 As shown, this invention provides a spatiotemporal fusion network architecture and training method guided by physical mechanisms;

[0078] Specifically:

[0079] The spatial aggregation layer employs a physically constrained graph attention mechanism, where the calculation of attention weights is based on the correlation of node feature learning and the physical weights pre-calculated based on hydrological connectivity.

[0080] It is important to note that for the target node, its effective neighbor set is strictly defined as all upstream nodes that point to it via directed edges, forcing information to be transmitted only upstream to downstream in the hydrological process. The attention weight calculation process for the target node to aggregate information from upstream nodes includes: linearly transforming and concatenating the feature vectors of the target node and upstream nodes respectively, and calculating the original relevance score based on data features through a neural network; then, the calculated feature relevance score is compared with the pre-calculated physical connection weight W from the upstream node to the target node. ba Multiplication is performed. The strength of physical connectivity is used as the modulation signal; even if two nodes are highly correlated, their attention will be forcibly suppressed if they lack a physical connection path (weights approaching 0). The modulated scores of all upstream nodes are normalized to obtain the final attention weights. The updated features of the target unit are the weighted sum of the features of all its upstream nodes.

[0081] The time evolution layer employs a gated recurrent unit network to independently process the state sequence of each node and output the predicted state for future time steps.

[0082] It is important to note that after the above aggregation, each spatial grid node yields a feature that incorporates information from upstream nodes. The feature sequence of each node across consecutive historical time steps is independently input into a gated recurrent unit (GRU) network. This GRU network learns and memorizes the time-varying patterns of each node's state and outputs a preliminary prediction of the state for future moments. This layer processes all nodes in parallel, efficiently capturing temporal dependencies.

[0083] The iterative spatiotemporal fusion module inputs the output of the time evolution layer back into the spatial aggregation layer to perform at least one round of time prediction-spatial correction loop to generate the final prediction result.

[0084] It is important to note that the preliminary predicted states of all nodes output by the temporal evolution layer at future moments are used as new node features and input again into the spatial aggregation layer. In the second round of aggregation, the nodes, based on these preliminary predicted states, aggregate spatial information from upstream nodes again according to the physically constrained attention mechanism, achieving a spatial correction. Ultimately, the output state at future moments is the result of at least one round of spatiotemporal iterative fusion.

[0085] The gradient-aware adaptive multi-task loss function is defined as follows:

[0086] S 总 =λ1×S 状态 +λ2×S 变化方向

[0087] Among them, S 状态 For state prediction loss, a loss function robust to outliers, such as the Huber loss function, can be used to supervise the spatiotemporal fusion network's accuracy in predicting the thickness at the end of the input sequence. Its expression is:

[0088]

[0089] in, Let this be the actual thickness value of sample i at the current time T. To simulate the thickness value, the parameter δ is the threshold for transforming the Huber loss from squared loss to linear loss, which can be based on the training set. Set the standard deviation, for example, to 1.5 times the standard deviation. 变化方向 To predict the loss in the direction of change, an interval loss function, such as the interval ranking loss function, can be used to guide the spatiotemporal fusion network to learn from the current time T to future times. The thickness variation trend is expressed as follows:

[0090]

[0091] in, To label the direction of real change, the true thickness value of the same spatial grid node at a known future time is calculated in the multi-temporal virtual-real fusion training set. Compared with the actual thickness value at the current moment The difference is obtained by taking its sign function; the interval value is used to control the judgment boundary for meaningful changes, which can be determined based on the statistics (such as the median) of the absolute value sequence of thickness changes of all samples in the training set. For sample i at future time The simulated soil thickness value. The multi-temporal virtual-real fusion training set is obtained by using historical climate, remote sensing, and other data to drive the geomorphological-soil process model in step S1, generating simulated thickness data for several past time steps, and fusing it with the measured thickness data for the corresponding period.

[0092] λ1 and λ2 are the weight coefficients of the two loss terms, both initially set to 1. During training, the weights are dynamically adjusted based on the instantaneous optimization difficulty of the loss for each task. Specifically, this includes: calculating the gradients of the state prediction loss and the change direction prediction loss with respect to all parameters of the spatiotemporal fusion network before each parameter update; calculating the L2 norm (i.e., the length of the gradient vectors) of these two gradient vectors; and dynamically updating the weights using an inverse proportional balancing rule.

[0093] λ 1新 =λ1×State prediction loss gradient norm / (Change direction prediction loss gradient norm + minimum value)

[0094] λ 2新 =λ² × gradient norm of change direction prediction loss / (gradient norm of state prediction loss + minimum value)

[0095] The minimum value is a very small positive constant introduced to prevent the denominator from being zero. Its value should be much smaller than the typical value of the gradient norm during training, such as 1e-8.

[0096] S4. The integrated feature set at the current moment is used as a node attribute and loaded into the spatiotemporal graph skeleton to form the current spatiotemporal graph. It is then input into the trained spatiotemporal fusion network and, combined with the preset soil thickness classification standard, a soil thickness type prediction result is generated. The time-series integrated feature set under the future scenario is used as a node attribute sequence to construct a future spatiotemporal graph sequence. This sequence is then input into the spatiotemporal fusion network for multi-step rolling prediction to generate a prediction result of the dynamic evolution trend of soil thickness type and a soil degradation risk warning information.

[0097] It is important to note that the method for generating the soil thickness type prediction results is as follows: Environmental data (such as climate or remote sensing imagery) at the current moment is processed through step S2 to obtain an integrated feature set for the current moment. This integrated feature set is used as node attributes and loaded into the spatiotemporal map skeleton constructed in the preceding steps to form the current spatiotemporal map. This map is then input into a pre-trained spatiotemporal fusion network. After forward propagation, the network outputs the predicted soil thickness value for each network node. A preset soil thickness classification standard is used, such as extremely thin layer (0-15cm), thin layer (15-30cm), medium layer (30-60cm), thick layer (60-100cm), and extremely thick layer (>100cm). Based on the predicted thickness value, each grid is assigned to the corresponding thickness type, generating a spatial distribution map of soil thickness types for the entire region at the current moment.

[0098] The specific method for generating the dynamic evolution trend results of soil thickness type and the early warning information of soil degradation risk is as follows:

[0099] Acquire future scenario data, such as precipitation, temperature series, and land use change maps. For each future time step (e.g., every 5 years), recalculate the integrated feature set in S2 using this data, especially climate-related and vegetation cover-related features, to form a future time-series integrated feature set sequence. Load the integrated feature set of the first time step (e.g., 2025) as node attributes into the spatiotemporal map skeleton to form the initial future spatiotemporal map. Then, load the current time step t... 当前 The future spacetime graph is input into the spacetime fusion network, and the next time step t is output. 当前 +1 Predicted soil thickness value H t当前+1 The predicted soil thickness value is used as a new important node feature to update the time step t. 当前 +1 to the node attribute. The updated t 当前 The attributes of the node at time +1 are loaded into the spatiotemporal graph skeleton to construct a new future spatiotemporal graph, which serves as the input for the next round of prediction. This process is repeated until the prediction for the preset time period (e.g., 2100) is completed.

[0100] After completing the rolling simulation, the predicted thickness sequence for each grid at different future time steps is obtained. This sequence is then analyzed as follows: First, for each grid, the path of its thickness type change over time is analyzed. At the regional scale, the area or proportion of grids where soil thickness type changes from A to B during the simulation period is statistically analyzed, forming a transition matrix. This matrix quantitatively describes the main direction and intensity of soil degradation or improvement. Then, several key future time nodes (such as 2030, 2050, and 2100) are selected, and the predicted thickness values ​​for each node are mapped according to the same classification standard, generating a series of spatial distribution maps of soil thickness types. By comparing these maps, the dynamic evolution trend of soil thickness types is revealed.

[0101] Based on the simulation results, soil degradation risks are identified and early warnings are issued. For each grid, the annual average rate of change of soil thickness is calculated using the predicted thickness values ​​at the beginning and end of the simulation. The formula is: Annual average rate of change of soil thickness = (Ending thickness value - Initial thickness value) / Number of simulation years.

[0102] Negative values ​​indicate degradation, while positive values ​​indicate accumulation. A degradation risk threshold is preset, such as an average annual thinning rate exceeding 0.5 cm / year. All grids with an average annual thinning rate exceeding this threshold are identified as high-risk degradation zones. By combining the spatial distribution, contiguousness, and topographical location of these zones (e.g., whether they are located on steep slope farmland or in important ecological function zones), a spatial distribution map of soil degradation risk warning is generated.

[0103] To gain a deeper understanding of the decision-making basis of the spatiotemporal fusion network and quantify the contribution of each environmental factor to the spatial differentiation of soil thickness, model decision analysis is performed after the spatiotemporal fusion network has been trained or applied. This embodiment performs model decision analysis after the spatiotemporal fusion network has been applied.

[0104] It is important to note that this step is based on the Shapley additive interpretation framework and adapted to the architectural characteristics of the spatiotemporal fusion network. Specifically, it includes: selecting the spatiotemporal graph input to be interpreted, such as a graph representing the current situation or a graph under a future scenario, and obtaining its node feature matrix and spatiotemporal graph skeleton. Baseline values ​​of node features are sampled and calculated from the training set. By switching specific features of the target node between their true values ​​and baseline values, and considering their influence on downstream nodes through graph edges, a series of intervention samples are constructed. Each intervention sample is input into the trained spatiotemporal fusion network to obtain the predicted value changes of the target node. Using the Shapley value calculation formula, the contribution of each feature to the predicted value of each node is estimated through sampling approximation. The calculated contributions of all nodes are mapped back to geographic space to generate a spatial distribution map of the contribution of each feature. By analyzing these maps, the key features that dominate the formation of high soil thickness areas and their direction of action (positive / negative contribution) are identified, thereby demonstrating the dominant controlling factors and geographic variation patterns of spatial differentiation of soil thickness.

[0105] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended embodiments and their equivalents.

Claims

1. A method for predicting soil thickness type based on feature set ensemble learning, characterized by, include: S1. Based on the measured thickness data of the target area, calibrate the parameterized geomorphic-soil process model, generate simulated thickness data, and fuse it with the measured thickness data to generate a virtual-real fusion training set; S2. Extract the basic environmental variables of the target area and construct surrogate features of key soil-forming processes based on soil genesis functions; use constrained causal analysis to select a subset of features that are causally related to soil thickness from the basic environmental variables and surrogate features to form an integrated feature set. S3. Using spatial grids as nodes, construct a spatiotemporal graph skeleton with physically meaningful edge weights based on hydrological flow direction; input the integrated feature set of multiple time phases as node attributes into the spatiotemporal fusion network guided by physical mechanisms, and train the network using the virtual-real fusion training set as supervision labels; the network achieves prediction through iterative spatial aggregation and temporal evolution, and is optimized using a gradient-aware adaptive multi-task loss function. S4. The integrated feature set at the current moment is used as a node attribute and loaded into the spatiotemporal graph skeleton to form the current spatiotemporal graph. The graph is then input into the trained spatiotemporal fusion network and, combined with the preset soil thickness classification standard, a soil thickness type prediction result is generated. The time-series integrated feature set under future scenarios is used as a node attribute sequence to construct a future spatiotemporal graph sequence. This sequence is then input into the spatiotemporal fusion network for multi-step rolling prediction to generate a prediction result of the dynamic evolution trend of soil thickness type and a soil degradation risk warning information. The spatiotemporal fusion network guided by the physical mechanism in S3 specifically includes: The spatial aggregation layer employs a physically constrained graph attention mechanism, where the calculation of attention weights is based on the correlation of node feature learning and the physical weights pre-calculated based on hydrological connectivity. The time evolution layer employs a gated recurrent unit network to independently process the state sequence of each node and output the predicted state for future time steps. The iterative spatiotemporal fusion module inputs the output of the time evolution layer back into the spatial aggregation layer, performs at least one round of time prediction-spatial correction loop, and generates the final prediction result. The construction and optimization method of the gradient-aware adaptive multi-task loss function in S3 includes: the loss function is composed of a weighted sum of a state prediction loss term and a change direction prediction loss term; the weight coefficients of the loss function are dynamically adjusted inversely during training based on the gradient magnitude information generated by each loss term on the spatiotemporal fusion network parameters. 2.The feature set ensemble learning based soil thickness type prediction method according to claim 1, characterized in that: S1 specifically includes: S11. Obtain the measured soil thickness data of the profile and clean it using an anomaly detection method combined with spatial autocorrelation. S12. Based on the dominant soil-forming process in the target area, select a parameterized geomorphological-soil process model and set prior probability distributions for its key parameters. S13. Divide the target area into geomorphic units. For each geomorphic unit, use the measured thickness data within the unit to solve the posterior probability distribution of the model parameters through variational inference algorithm to complete the probabilistic calibration of the model. S14. Based on the posterior distribution of the calibrated parameters, generate simulated thickness data and its uncertainty; construct and train a confidence evaluation neural network, and output the confidence weight of the simulated point according to the local environmental characteristics, simulation uncertainty and nearby measured residuals. S15. Combine the measured thickness data with the weighted simulated data to form a virtual-real fusion training set. 3.The feature set ensemble learning based soil thickness type prediction method according to claim 1, wherein: S2 specifically includes: S21. Extract the basic environmental variables of the target area to form a basic environmental feature layer; based on the preset soil genesis function, calculate the proxy features that quantify the intensity of key soil-forming processes from the basic environmental variables to form a process proxy feature layer. The basic environmental variables include climate variables, parent material type variables, and slope variables; the proxy characteristics include the erosion transport capacity index. S22. At the measured soil thickness points, collect data of the corresponding basic feature layer and process proxy feature layer; construct a domain knowledge constraint matrix, and generate a causal relationship network between variables through a constraint causal discovery algorithm based on the matrix; when searching for the causal relationship network structure, the constraint causal discovery algorithm treats the constraints in the matrix as inviolable hard conditions. The specific methods for constructing the domain knowledge constraint matrix include: S221. Based on the soil thickness modeling requirements, determine all candidate basic environmental variables and process proxy variables, and construct a feature set; S222. Construct a square constraint matrix corresponding to the number of features in the feature set; each element in the matrix represents a preset prior constraint type for the causal relationship between feature pairs. S223. Based on the basic principles of soil genesis, set core constraints in the matrix; the core constraints include at least: setting the corresponding element of the climate variable pointing to the parent material type variable to -1, indicating that climate cannot directly change the geological parent material; setting the corresponding element of the slope variable pointing to the erosion transport capacity index to 1, indicating that slope directly affects erosion transport capacity; setting the corresponding element of the erosion transport capacity index pointing to the soil thickness to 1, indicating that erosion transport capacity directly affects soil thickness. S224. After assigning constraint values ​​to all feature pairs, output the domain knowledge constraint matrix. S23. From the causal relationship network, identify features that have direct or indirect causal paths with soil thickness and form a feature subset; based on the feature subset, extract and integrate feature data on the whole area grid to form an integrated feature set. 4.The feature set ensemble learning based soil thickness type prediction method according to claim 1, wherein: The S3 section describes the construction of a spatiotemporal graph skeleton with physically meaningful edge weights based on hydrological flow direction. Specifically, this includes: calculating the flow direction of each grid node based on a digital elevation model; establishing a directed edge from the node itself to its downstream adjacent grid node for any given grid node; and calculating the weight of the directed edge based on the elevation difference and horizontal distance between upstream and downstream nodes. 5.The feature set ensemble learning based soil thickness type prediction method according to claim 1, wherein: The prediction results of the dynamic evolution trend of soil thickness type and the early warning information of soil degradation risk in S4 are generated by the following methods: S41. Obtain the temporal environment data under the future scenario. Step S2 processes the data to obtain the temporal integrated feature set of each future time step and constructs the initial future spatiotemporal map. S42. Input the future spatiotemporal map of the current time step into the spatiotemporal fusion network to predict the soil thickness value of the next time step; use the predicted soil thickness value as a feedback feature to update the attributes of the corresponding node in the next time step. S43. Combine the updated node attributes and future scenario data to construct a new future spatiotemporal map as the input for the next round of prediction; repeat steps S42 and S43 until the prediction of the preset duration is completed. S44. The soil thickness type change sequence of each spatial grid within the statistical extrapolation period is used to generate a soil thickness type transition matrix and a dynamic evolution diagram of soil thickness type at key time nodes. S45. Based on the simulation results, calculate the annual average rate of change of soil thickness for each grid; identify areas where the rate of change exceeds the preset degradation threshold as risk degradation zones, and generate soil degradation risk warning information. 6.The feature set ensemble learning based soil thickness type prediction method according to claim 1, wherein: The method also includes model decision analysis, which uses the Shapley additive interpretation method based on the trained spatiotemporal fusion network to analyze the contribution of each feature in the integrated feature set to the prediction result.

Images