Machine learning-based large-scale prediction method and system for cadmium-arsenic in rice field soil

Abstract

This invention discloses a large-scale cadmium-arsenic prediction method and system for paddy soil based on machine learning, belonging to the field of soil environmental monitoring. The method includes: acquiring an initial dataset of the large-scale prediction area; mapping the geographically weighted non-stationary spatial distance matrix by introducing a kernel function based on spatial distance decay to calculate and generate a local spatial autocorrelation weight vector; constructing a geographically weighted ensemble tree model; obtaining the gridded spatial coordinates of the area to be predicted and the corresponding target soil physicochemical properties, inputting them into the geographically weighted ensemble tree model, and outputting a predicted cadmium-arsenic concentration distribution map. This invention effectively eliminates the distortion interference of spatial outliers on local modeling weight allocation and gradient estimation.

Description

Technical Field

[0001] This invention relates to the field of soil environmental monitoring, and in particular to a large-scale method and system for predicting cadmium-arsenic levels in paddy field soil based on machine learning. Background Technology

[0002] Heavy metal pollution in soil is a major problem threatening food security and the ecological environment. Cadmium (Cd) and arsenic (As), two typical harmful heavy metals found in paddy field soils, require precise spatial distribution prediction for effective risk management, safe utilization zoning, and precise remediation decisions. With the ongoing detailed investigation of soil pollution and the classification management of agricultural land, there is an urgent need to efficiently and accurately obtain spatial distribution information of cadmium and arsenic concentrations in paddy field soils over a large spatial scale. Traditional methods for spatial prediction of heavy metal pollution in soil mainly include geostatistical interpolation methods (such as Kriging interpolation) and classical machine learning regression methods (such as random forests, support vector machines, and gradient boosting decision trees). While geostatistical methods inherently possess spatial autocorrelation modeling capabilities, their reliance on the assumption of data stationarity severely limits their applicability under large-scale, highly spatially heterogeneous conditions. Classical machine learning methods, though possessing powerful nonlinear fitting capabilities, employ a globally uniform parameter estimation strategy during model construction, treating all spatially located samples as independent observations with the same statistical distribution. This ignores a fundamental fact widely verified in environmental geochemistry: the driving mechanism of soil physicochemical properties on heavy metal concentrations is not uniform across geographic space, but exhibits significant non-stationarity with spatial location. For example, in areas surrounding industrial and mining enterprises, soil pH may be the dominant driving factor influencing cadmium migration activity; while in natural background areas far from anthropogenic pollution sources, the geochemical composition of parent material may be the key factor determining cadmium background levels. This spatial non-stationarity is particularly prominent in large-scale prediction scenarios because the study area often spans multiple landform types, climatic conditions, and land use patterns, resulting in environmental gradient variations far greater than in local small-scale studies.

[0003] Currently, there is a lack of a soil heavy metal prediction method that can effectively characterize the spatial nonstationarity driving relationship of soil physicochemical parameters on cadmium-arsenic concentration distribution on a large-scale spatial scale, while taking into account the robust handling of spatial outliers, collaborative modeling of cadmium-arsenic geochemical coupling relationship, and the need for large-scale, efficient parallel computing. Summary of the Invention

[0004] One of the objectives of this invention is to provide a large-scale cadmium-arsenic prediction method and system for paddy soil based on machine learning, in order to solve the problem that existing machine learning models, which use a globally unified parameter estimation strategy, cannot capture the driving relationship between soil physicochemical properties and the spatial nonstationarity of cadmium-arsenic concentration.

[0005] This invention is achieved through the following technical solution: a large-scale cadmium-arsenic prediction method for paddy field soil based on machine learning, comprising the following steps: S1, obtaining an initial dataset for the large-scale prediction area, wherein the initial dataset includes the spatial coordinates of multiple known sampling points, a soil physicochemical property parameter matrix, and corresponding cadmium-arsenic concentration label values; S2, constructing a geographically weighted non-stationary spatial distance matrix based on the spatial coordinates, and mapping the geographically weighted non-stationary spatial distance matrix by introducing a kernel function based on spatial distance decay, thereby dynamically calculating and generating local spatial distances for any target prediction point within the large-scale prediction area. S3. Based on the initial dataset and the local spatial autocorrelation weight vector, construct a geographic weighted ensemble tree model; wherein, in the process of constructing the geographic weighted ensemble tree model, the local spatial autocorrelation weight vector and the derivative of the loss function are multiplicatively combined to calculate the geographic weighted gradient statistic, and the spatial weighted structure gain is maximized based on the geographic weighted gradient statistic as the node splitting criterion; S4. Obtain the gridded spatial coordinates of the area to be predicted and the corresponding target soil physicochemical property parameters, input them into the geographic weighted ensemble tree model, and output a cadmium-arsenic concentration distribution prediction map.

[0006] Furthermore, after step S1 and before step S2, a data cleaning step is also included: for the parameter values ​​of each known sampling point in the soil physicochemical property parameter matrix or the cadmium-arsenic concentration label values, the local Moran index of each known sampling point is calculated based on the spatial adjacency weight matrix. The local Moran index is used to measure the local spatial association pattern between each known sampling point and its spatial neighborhood. The local Moran index of each known sampling point is compared with the dynamic spatial saliency distribution table to identify spatial outliers that exhibit high spatial isolation characteristics. When constructing the geographic weighted non-stationary spatial distance matrix in step S2, the spatial distance between the spatial outliers and other spatial locations is forcibly set to infinity.

[0007] Furthermore, the calculation of the local Moran index includes: for any known sampling point, calculating the deviation between the parameter value of the known sampling point and the global mean; based on the weight elements corresponding to the known sampling point in the spatial adjacency weight matrix, performing weighted aggregation on the deviations between the parameter values ​​of all spatial neighbors of the known sampling point and the global mean; multiplying the deviation of the known sampling point by the weighted aggregation result, and normalizing it with the global variance to obtain the local Moran index of the known sampling point.

[0008] Furthermore, the identification of spatial outliers includes: when the local Moran index of a known sampling point is negative, and the absolute value of the negative value exceeds the critical value of the corresponding significance level in the dynamic spatial significance distribution table, the known sampling point is determined to be a spatial outlier exhibiting a high-low pattern or a low-high pattern.

[0009] Furthermore, the kernel function based on spatial distance attenuation is an adaptive bisquare spatial attenuation kernel function, which is constructed based on an adaptive bandwidth variable. The adaptive bandwidth variable is a distance parameter determined individually for each target prediction point, used to control the effective range of spatial weight attenuation. The adaptive bisquare spatial attenuation kernel function has a tight support characteristic. When the distance between the target prediction point and the known sampling point reaches or exceeds the adaptive bandwidth variable, the corresponding weight element in the local spatial autocorrelation weight vector is zero.

[0010] Furthermore, the calculation steps for each weight element in the local spatial autocorrelation weight vector include: calculating the Euclidean spatial distance between the target prediction point and each known sampling point; calculating the ratio of the Euclidean spatial distance to the adaptive bandwidth variable; when the ratio is less than one, using one minus the square of the ratio as the base, squaring the base to obtain the corresponding weight element; when the ratio is greater than or equal to one, the corresponding weight element is zero.

[0011] Furthermore, the adaptive bandwidth variable is determined through the following dynamic adjustment steps: the adaptive bandwidth variable is initialized to a basic distance threshold; the number of valid known sampling points whose distance to the target prediction point is less than the current adaptive bandwidth variable is counted; if the number of valid known sampling points is less than a preset minimum spatial degree of freedom threshold, the current adaptive bandwidth variable is amplified and recounted according to a preset step size coefficient, and the process is iterated until the number of valid known sampling points reaches or exceeds the preset minimum spatial degree of freedom threshold.

[0012] Furthermore, after step S2 and before step S3, the method further includes: a targeting constraint step based on the synergistic and antagonistic relationship between cadmium and arsenic, wherein for each target prediction point, a weighted Pearson correlation coefficient between the cadmium concentration label value and the arsenic concentration label value is calculated based on the local spatial autocorrelation weight vector; a targeting penalty term is calculated based on the weighted Pearson correlation coefficient; and when calculating the spatial weighted structural gain in step S3, the targeting penalty term is subtracted from the spatial weighted structural gain.

[0013] Furthermore, the calculation of the weighted Pearson correlation coefficient includes:

[0014] For each known sampling point within the local spatial range of the target prediction point, using the weight elements in the local spatial autocorrelation weight vector as weights, calculate the first weighted arithmetic mean of the cadmium concentration label value and the second weighted arithmetic mean of the arsenic concentration label value, respectively; calculate the difference between the cadmium concentration label value and the first weighted arithmetic mean for each known sampling point, and the difference between the arsenic concentration label value and the second weighted arithmetic mean for each known sampling point; sum the products of the two differences using the weight elements to obtain the weighted covariance; sum the squares of the two differences using the weight elements and take the square root to obtain the weighted standard deviation of cadmium concentration and the weighted standard deviation of arsenic concentration; divide the weighted covariance by the product of the weighted standard deviation of cadmium concentration and the weighted standard deviation of arsenic concentration to obtain the weighted Pearson correlation coefficient.

[0015] Furthermore, the calculation of the targeted penalty term includes: calculating the absolute value of the weighted Pearson correlation coefficient; performing an exponential function operation with the negative value of the absolute value as the exponent to obtain an exponential mapping value; and multiplying the exponential mapping value by the targeted penalty coefficient to obtain the targeted penalty term.

[0016] Further, the derivative of the loss function includes the first derivative and the second derivative. The steps for calculating the geographically weighted gradient statistics include: for each known sampling point within the current node, calculating the first and second derivatives of the global loss function relative to the current predicted value of the model; performing the multiplicative combination operation between each weight element in the local spatial autocorrelation weight vector and the first derivative of the corresponding known sampling point to obtain the weighted first derivative of each known sampling point; performing the multiplicative combination operation between each weight element in the local spatial autocorrelation weight vector and the second derivative of the corresponding known sampling point to obtain the weighted second derivative of each known sampling point; according to the candidate splitting scheme, accumulating the weighted first derivatives corresponding to the samples to be divided into the left subtree to obtain the weighted first gradient sum of the left subtree, and accumulating the weighted second derivatives corresponding to the samples to be divided into the left subtree to obtain the weighted second gradient sum of the left subtree; similarly calculating the weighted first gradient sum and the weighted second gradient sum of the right subtree.

[0017] Further, the calculation of the spatially weighted structural gain includes: dividing the square of the weighted first-order gradient sum of the left subtree by the sum of the weighted second-order gradient sum of the left subtree and the regularization coefficient to obtain the left subtree gain component; dividing the square of the weighted first-order gradient sum of the right subtree by the sum of the weighted second-order gradient sum of the right subtree and the regularization coefficient to obtain the right subtree gain component; dividing the square of the sum of the weighted first-order gradient sum of the left subtree and the sum of the weighted first-order gradient sum of the right subtree by the sum of the weighted second-order gradient sum of the left subtree, the weighted second-order gradient sum of the right subtree, and the regularization coefficient to obtain the gain component before splitting; adding the gain component of the left subtree and the gain component of the right subtree, subtracting the gain component before splitting, multiplying by half, and subtracting the leaf node splitting penalty constant to obtain the spatially weighted structural gain.

[0018] Further, step S3 may also include the step of determining the weights of the leaf nodes: for a leaf node that has completed splitting, multiply each weight element in the local spatial autocorrelation weight vector by the first derivative of each sample falling into the leaf node and sum them up to obtain the weighted first-order gradient sum of the leaf node; multiply each weight element in the local spatial autocorrelation weight vector by the second derivative of each sample falling into the leaf node and sum them up to obtain the weighted second-order gradient sum of the leaf node; divide the negative value of the weighted first-order gradient sum of the leaf node by the sum of the weighted second-order gradient sum of the leaf node and the regularization coefficient to obtain the leaf node weight of the leaf node.

[0019] Furthermore, a distributed computing mechanism is employed in step S3: using a spatial quadtree data structure, the large-scale prediction area is divided into multiple local sub-computation domains, with an overlap width set between adjacent local sub-computation domains; data fragments corresponding to each local sub-computation domain are allocated to independent computing nodes, and local geographic weighted ensemble tree models are constructed in parallel within each independent computing node. Further, for target prediction points located within the overlap width range, a spatial inverse distance weighted mixing function is used to fuse the prediction results output by adjacent local sub-computation domains: the distance from the target prediction point to the geometric center of each local sub-computation domain covering the target prediction point is calculated; using the negative powers of each distance as weights, a weighted average is applied to the predicted values ​​output by the geographic weighted ensemble tree models corresponding to each local sub-computation domain for the target prediction point, yielding the fused prediction value for the target prediction point.

[0020] Further, the step S4 of inputting the data into the geographic weighted ensemble tree model for prediction includes: for each target prediction point corresponding to the gridded spatial coordinates in the area to be predicted, calculating the local spatial autocorrelation weight vector of the target prediction point based on the distance relationship between the target prediction point and each known sampling point; using the local spatial autocorrelation weight vector as a spatial modulation parameter, inputting the target soil physicochemical property parameters corresponding to the target prediction point into the geographic weighted ensemble tree model; multiplying the weights of the leaf nodes output by each decision tree in the geographic weighted ensemble tree model by the learning rate and then summing them to obtain the predicted cadmium-arsenic concentration value of the target prediction point; and assembling the predicted cadmium-arsenic concentration values ​​of all target prediction points into a rasterized array according to their spatial location to generate the predicted cadmium-arsenic concentration distribution map.

[0021] Another aspect of the present invention provides a large-scale cadmium-arsenic prediction system for paddy soil based on machine learning, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the large-scale cadmium-arsenic prediction method for paddy soil based on machine learning as described above.

[0022] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0023] 1. This invention prioritizes fitting the local data distribution patterns in the neighborhood surrounding the target prediction point when evaluating candidate splitting schemes, thereby adaptively characterizing the spatial nonstationarity driving relationship between soil physicochemical parameters and cadmium-arsenic concentration distribution over a large spatial scale. This overcomes the inherent defect of traditional global unified parameter estimation strategies that cannot capture spatial heterogeneity, and significantly improves the prediction accuracy of cadmium-arsenic concentration under large-scale spatial nonstationarity conditions.

[0024] 2. This invention can identify spatial outliers that are not significantly abnormal in global statistical distribution but exhibit inconsistent high-low or low-high patterns with their spatial neighborhood attribute values. By forcibly setting the distance mapping value of spatial outliers to infinity, soft deletion is achieved. That is, without physically removing the original data, the spatial weight of spatial outliers after kernel function mapping is deterministically set to zero. Thus, while preserving data integrity and traceability, it effectively eliminates the distortion interference of spatial outliers on local modeling weight allocation and gradient estimation, avoids non-physical jumps in prediction results in transitional zones with drastic spatial variations, and improves the spatial robustness of the model.

[0025] 3. This invention employs an adaptive bandwidth dynamic search mechanism to dynamically adjust the effective range of the spatial weight decay kernel function based on the local sampling density of each target prediction point. This ensures that regardless of whether the target prediction point is located in a densely sampled or sparsely sampled region, a sufficient number of effective known sampling points are always included in its local spatial gravitational field to guarantee the statistical degrees of freedom for local model parameter estimation. This solves the bandwidth selection dilemma caused by the uneven spatial distribution of sampling density in large-scale prediction regions, avoids the local model underdeterminacy problem caused by insufficient information in sparsely sampled regions, and ensures the robustness and statistical reliability of the prediction results throughout the entire study area. Attached Figure Description

[0026] The accompanying drawings, which are included to provide a further understanding of embodiments of the invention and form part of this application, do not constitute a limitation thereof. In the drawings:

[0027] Figure 1 This is a flowchart of the method provided in Embodiment 1 of the present invention.

[0028] Figure 2 This is a comparison chart of the spatial anomaly detection effects provided in Embodiment 1 of the present invention.

[0029] Figure 3 This is a four-quadrant scatter plot for local Moran's index spatial anomaly detection provided in Embodiment 1 of the present invention.

[0030] Figure 4 This is a comparison chart of spatial weight decay provided in Embodiment 1 of the present invention.

[0031] Figure 5 This is a comparison chart of the coverage effects of adaptive bandwidth dynamic search and fixed bandwidth search in dense areas, as provided in Embodiment 1 of the present invention.

[0032] Figure 6 This is a comparison chart of the coverage effects of adaptive bandwidth dynamic search and fixed bandwidth search in sparse regions provided in Embodiment 1 of the present invention.

[0033] Figure 7 The response curve of the cadmium-arsenic local correlation intensity provided in Embodiment 1 of the present invention.

[0034] Figure 8 This is a schematic diagram comparing the selection of splitting thresholds provided in Embodiment 1 of the present invention.

[0035] Figure 9 This is a comparison chart of the computational complexity of the double logarithmic axis provided in Embodiment 1 of the present invention. Detailed Implementation

[0036] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, 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. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0037] Example 1

[0038] This embodiment discloses a large-scale cadmium-arsenic prediction method for paddy field soil based on machine learning. Figure 1 The overall method flowchart of this embodiment is shown. As can be seen from the figure, this embodiment includes the following steps:

[0039] S1: Obtain the initial dataset for the large-scale prediction area. The initial dataset includes the spatial coordinates of multiple known sampling points, the soil physicochemical property parameter matrix, and the corresponding known cadmium-arsenic concentration label values.

[0040] Among them, the large-scale prediction area refers to the target research area that covers a large range in the geographical space dimension. This area can be one or more provincial-level concentrated areas of paddy fields, or it can be the coverage area of ​​a national soil pollution survey that spans multiple geographical and climatic zones. For example, the large-scale prediction area can be multiple major rice-producing areas in the middle and lower reaches of the Yangtze River Plain, the paddy field belts around industrial and mining enterprises in the Pearl River Delta region, or the scope of a national soil environmental quality survey covering multiple provinces.

[0041] The initial dataset refers to a set of structured raw data collected, organized, and summarized from multiple known sampling points within a large-scale prediction area before formally constructing the geographic weighted ensemble tree model. This dataset constitutes the basic data source for all subsequent spatial modeling and predictive analysis.

[0042] Spatial coordinates refer to a set of positioning parameters used to uniquely determine the absolute location of each known sampling point in geographic space; spatial coordinates include, but are not limited to, longitude, latitude, and elevation information. In the case of a two-dimensional plane, spatial coordinates are a binary pair containing longitude and latitude. , Longitude In three dimensions, spatial coordinates can be further incorporated with the elevation of the sampling points to form a triple. , Altitude.

[0043] A soil physicochemical property parameter matrix is ​​a multidimensional set of parameters organized in matrix form to quantitatively describe the soil environmental attributes at the location of each known sampling point. Each row of the matrix corresponds to a known sampling point, and each column corresponds to a soil physicochemical property index. For example, soil physicochemical property parameters include, but are not limited to, soil pH, soil organic matter content (SOM), cation exchange capacity (CEC), soil texture components (percentage of sand, silt, and clay), redox potential (Eh), iron and manganese oxide content, soil moisture content, and the content of nutrients such as total nitrogen, total phosphorus, and total potassium.

[0044] The known cadmium-arsenic concentration label values ​​refer to the concentration values ​​of cadmium (Cd) and arsenic (As) in soil obtained by actual measurement at each known sampling point through field sampling and laboratory chemical analysis. These concentration label values ​​serve as the supervised learning target in the model training process, guiding the geographically weighted ensemble tree model to learn the spatial non-stationary mapping relationship between soil physicochemical properties and heavy metal concentrations.

[0045] It is understandable that the spatial coordinates, soil physicochemical properties, and cadmium-arsenic concentration labels of each known sampling point in the initial dataset are one-to-one correspondences. That is, each sampling point simultaneously possesses definite geographical location information, a complete description of environmental attributes, and a true value of heavy metal concentration verified by the laboratory.

[0046] In this embodiment, the initial dataset may be provided by the National Soil Environmental Quality Monitoring Network, regional soil pollution status detailed investigation projects, or field sampling surveys organized by research institutions.

[0047] After obtaining the initial dataset, a non-black-box data cleaning step can be included to eliminate spatial outliers. The necessity of this data cleaning step stems from the fact that various data quality issues are inevitably introduced during the soil field sampling process: some sampling points may exhibit extreme values ​​that are inconsistent with their surrounding spatial environment due to human error, laboratory analysis bias, or the instantaneous impact of local point source contamination on the sample.

[0048] Traditional data cleaning paradigms typically employ threshold truncation methods based on global statistical characteristics, such as the 3σ principle, which removes samples that deviate from the global mean by more than three standard deviations to identify and eliminate outliers. However, these methods have a fundamental logical blind spot: they only focus on anomalies in the numerical domain, completely ignoring anomalies in the spatial domain.

[0049] Specifically, the concentration of a soil heavy metal sampling point may fall entirely within the reasonable range of the global statistical distribution. However, if the point exhibits concentration characteristics that are diametrically opposed to its immediate neighbors—for example, if the point itself has an extremely high concentration while all its neighbors have extremely low concentrations (i.e., a high-low anomaly pattern, or HL pattern), or if the point itself has an extremely low concentration while all its neighbors have extremely high concentrations (i.e., a low-high anomaly pattern, or LH pattern)—then the point constitutes a typical spatial outlier. If such spatial outliers are included in the subsequent spatial weighted modeling process without discrimination, it will severely distort and mislead the model's weight allocation and gradient estimation in the neighborhood of that point, thereby causing violent oscillations and instability in the local model parameters. Ultimately, this will lead to non-physical jumps in the prediction results in the transitional zone of drastic spatial variation that do not conform to geochemical laws.

[0050] Therefore, to rigorously define and identify such spatial anomalies mathematically, the data cleaning step in this embodiment introduces the classic Local Moran's Index (LMI) from the field of spatial autocorrelation analysis as a diagnostic tool. The LMI is a statistical measure used in spatial autocorrelation analysis to measure the local spatial correlation pattern between a sampling point and its spatial neighborhood. Its core idea is that it not only measures the deviation of a sampling point's own attribute value from the global mean, but also, by introducing a spatial adjacency weight matrix, weights and aggregates the deviations of the point's attribute values ​​from those of its spatial neighbors, thereby characterizing the local spatial correlation pattern between the point and its spatial neighborhood. When a sampling point exhibits a significant high-low or low-high anomaly pattern, its LMI will obtain a significant negative value, which statistically provides a quantifiable basis for spatial anomaly judgment with a clear confidence level.

[0051] In this embodiment, the data cleaning step may specifically include:

[0052] 1) Calculate the local Moran index of soil physicochemical properties at each known sampling point.

[0053] For example, in this embodiment, for sampling points Local Moran index The following formula is used to calculate:

[0054]

[0055] in, Indicates the first The local Moran index calculated from each soil sampling point reveals the correlation pattern and strength between the point and its spatial neighborhood, based on the sign of the index and its absolute value. Indicates at the sampling point The specific value of a certain soil physicochemical parameter or heavy metal content measured at the site is the core object of the local Moran index calculation; This indicates that the soil parameter applies to all areas within the study region. The global arithmetic mean of the sampling points serves as a global reference benchmark, used to measure the direction and magnitude of each sampling point's deviation from the overall average level. This indicates the total number of soil sampling points within the study area; The spatial adjacency weight matrix represents the connection of sampling points. With sampling points The element value, this weight is usually defined as the reciprocal of the Euclidean distance between the two points, or based on the K-nearest neighbor criterion, when the sampling points Sampling point When the K nearest neighbors are By taking a fixed positive value and otherwise taking zero, the spatial adjacency weight matrix mathematically formalizes the core geographical concept of spatial proximity.

[0056] It is understandable that when the local Moran's index is positive and the absolute value is large, it indicates that the sampling point and its spatial neighbors exhibit a high-high or low-low spatial clustering pattern, i.e., a positive spatial correlation. When the local Moran's index is negative and the absolute value is large, it indicates that the sampling point and its spatial neighbors exhibit a high-low or low-high spatial dispersion pattern, i.e., a negative spatial correlation. In this case, the sampling point is likely to constitute a spatial outlier.

[0057] 2) After calculating the local Moran index for each sampling point, a conditional decision flow is further constructed. The local Moran index of each known sampling point is compared with the dynamic spatial significance distribution table. If a known sampling point exhibits high spatial isolation characteristics of high-low (HL) or low-high (LH) with confidence exceeding the limit, the distance mapped to the point in the spatial distance matrix is ​​forcibly adjusted to infinity before generating the local spatial autocorrelation weight vector, thereby truncating its abnormal spatial interference on the modeling of surrounding local nonstationarity.

[0058] For example, in this embodiment, the updated spatial distance mapping function As shown in the following formula:

[0059]

[0060] in, This indicates the update from the sampling point after spatial anomaly diagnosis. A generalized distance mapping function from the origin to any other spatial location; the output of this function will be directly used as the input of the subsequent geographic weighted kernel function. This represents the confidence threshold used to determine whether a spatially anomalous pattern is significant. This threshold is typically selected based on the statistical distribution of the local Moran's index under the null hypothesis (i.e., complete spatial randomness), corresponding to a specific significance level (e.g., ...). The critical value of ). Indicates sampling point The original Euclidean or geodesic distance between other spatial locations; symbol This represents positive infinity. When a sampling point is identified as a spatial outlier, its distance from all other locations is forcibly set to infinity. Figure 2 This figure shows a comparison of the spatial anomaly detection results of the data cleaning steps in this embodiment with those of the traditional 3σ method. Figure 2 The 3σ method misses outliers that are within a reasonable range but spatially abnormal; while the method in this embodiment accurately detects all three spatial outliers. Figure 3 The four-quadrant scatter plot of the local Moran's index spatial anomaly detection in this embodiment is shown. Figure 3 The X-axis represents the standardized z-score of each sample point, and the Y-axis represents the standardized spatial lag value (i.e., the weighted average of the spatial neighbor attribute values). The graph is divided into four quadrants centered on the origin: Quadrant 1 (HH: High-High Clustering), Quadrant 2 (LH: Low-High Anomaly, Spatial Outliers), Quadrant 3 (LL: Low-Low Clustering), and Quadrant 4 (HL: High-Low Anomalies, Spatial Outliers). Spatial outliers exceeding the confidence threshold in the LH and HL quadrants are marked with different colors.

[0061] Understandably, by forcibly setting the distance mapping value of spatial outliers to infinity, when this infinite value is substituted into the bisquared spatial decay kernel function (to be introduced in subsequent steps), the spatial weights output by the kernel function will deterministically and unambiguously decay to zero due to the domain limitation of the kernel function. This mechanism achieves soft deletion of spatial outliers, while the original data of the point remains in the dataset. However, when building a local model for any target prediction point, the contribution weight of the outlier is zero, thus effectively excluding it from the local modeling process. Compared to traditional hard deletion methods (i.e., directly physically removing outlier records from the dataset), this soft deletion mechanism offers greater flexibility and traceability, and ensures that the model's weight allocation exhibits a smooth and gradual transition in transitional zones with drastic spatial gradient variations, avoiding the risks of local weight overfitting and non-physical jumps in prediction results caused by outliers.

[0062] The dynamic spatial significance distribution table is a lookup table constructed based on the theoretical statistical distribution of the local Moran index under the null hypothesis of complete spatial randomness. Each entry in the distribution table records the critical value at a given significance level, which is used to determine whether the local Moran index of each sampling point reaches the level of statistical significance anomaly.

[0063] S2: Based on spatial coordinates, construct a geographically weighted non-stationary spatial distance matrix to cope with large-scale environmental gradients. By introducing a kernel function based on spatial distance decay to map the spatial distance matrix, dynamically calculate and generate local spatial autocorrelation weight vectors for any target prediction point in the prediction area.

[0064] The geographically weighted non-stationary spatial distance matrix is ​​a two-dimensional numerical structure organized in matrix form, recording the geographic spatial distance relationships between any target prediction point and all known sampling points. Each element of this matrix represents the Euclidean or geodesic distance between a pair of spatial locations, serving as the basic data input for subsequent spatial weight allocation. This matrix is ​​called a non-stationary spatial distance matrix because its purpose is not merely to record a globally uniform distance relationship, but rather to serve as a basis for subsequent mapping through kernel functions to local weight vectors with spatial heterogeneity. That is, for different target prediction points, the weights assigned to the same known sampling point are different and dynamically change with spatial location.

[0065] The kernel function based on spatial distance decay is a mathematical function that maps the continuous physical quantity of spatial distance to normalized spatial weights. The core characteristic of this function is that its output value (i.e., spatial weights) decreases monotonically as the input value (i.e., spatial distance) increases, and strictly decays to zero after the distance exceeds a certain critical threshold.

[0066] The local spatial autocorrelation weight vector is a one-dimensional vector generated by a kernel function mapping for a specific target prediction point within the prediction region. It contains the spatial weight relationships between the target point and all known sampling points. Each element of this weight vector takes a value in a closed interval between zero and one, where the weight elements corresponding to known sampling points closer to the target prediction point approach one, and the weight elements corresponding to known sampling points farther from the target prediction point approach zero.

[0067] In this embodiment, the specific steps for generating the local spatial autocorrelation weight vector may include:

[0068] S2.1: For any target prediction point Calculate its relationship with all known sampling points. European spatial distance between .

[0069] Among them, target prediction points This refers to any spatial location within a large-scale prediction area where cadmium-arsenic concentration prediction is needed. This location can be the center coordinates of a gridded raster cell or any user-specified geographic coordinate point. (Euclidean spatial distance) Refers to the target prediction point With known sampling points In geographic coordinate space, this straight-line distance measure reflects the absolute separation between two spatial locations. In a two-dimensional plane, this distance is calculated using the Euclidean distance formula based on the longitude and latitude coordinates of the two points; in a three-dimensional plane, elevation differences can be further incorporated.

[0070] S2.2: Set an adaptive bandwidth variable An adaptive bisquared spatial decay kernel function is constructed based on the adaptive bandwidth variable.

[0071] Among them, adaptive bandwidth variable This refers to the target prediction point A separately determined distance parameter that controls the effective range of the spatial weight decay function. This variable is called adaptive because its value is not a globally uniform fixed constant, but is dynamically adjusted according to the local sampling density of each target prediction point. In densely sampled regions, the bandwidth is smaller to preserve fine spatial resolution; in sparsely sampled regions, the bandwidth is larger to ensure sufficient information sources participate in local modeling.

[0072] The adaptive bi-square spatial decay kernel function is a nonlinear function that maps the ratio of spatial distance to adaptive bandwidth to spatial weights. The mathematical form of this function is bi-square. Its key characteristic is that it has a bounded compact support set. When the distance between two points exceeds the bandwidth, the kernel function output is strictly zero (rather than an infinitesimal quantity approaching zero). This means that the spatial weight matrix is ​​naturally sparse, which is beneficial for reducing the storage overhead and computational complexity of large-scale computations.

[0073] For example, in this embodiment, when At that time, the sample For target prediction points Local spatial autocorrelation weights The following formula is used to calculate:

[0074]

[0075] when At that time, the local spatial autocorrelation weight .

[0076] in, Indicates the target prediction point In a local spatial gravitational field constructed around a central point, the known sampling points The assigned spatial autocorrelation weights, whose values ​​range from a closed interval. , where 1 represents the maximum influence (obtained when the sampling point coincides with the target point), and 0 represents zero influence (obtained when the sampling point is outside the bandwidth range); This represents the normalized ratio of the distance between two points relative to the adaptive bandwidth. This ratio maps the distances to different target points and different bandwidths to a unified value. On comparable scales, the outer squaring operation provides the kernel function with second-order smooth boundary conditions, ensuring that the derivative of the weights at the bandwidth boundary is zero, thereby achieving a smooth and gradual transition in space and avoiding the introduction of artificial step artifacts in the prediction results due to abrupt changes in weights. Figure 4 This diagram shows a comparison of the spatial weight decay of the double-squared kernel function with other kernel functions in this embodiment. Figure 4 As can be seen from this, the compact support characteristic of the bisquare kernel function in this embodiment makes the weight matrix naturally sparse, and its storage and computational costs are far superior to those of the Gaussian kernel in large-scale computation.

[0077] S2.3: Predict the target point All corresponding Combined into a local spatial autocorrelation weight vector .

[0078] Among them, the local spatial autocorrelation weight vector This refers to the target prediction point With all A one-dimensional vector is formed by arranging the spatial weight values ​​calculated between the known sampling points in the order of their indices. This vector fully encodes the target prediction point. The spatial neighborhood structure and influence distribution pattern within the entire study area serve as the foundational weight inputs for core operations such as gradient weighting and feature splitting evaluation during the subsequent construction of the geographic weighted ensemble tree model.

[0079] Understandably, due to the compact support property of the bisquare kernel function, for most known sampling points far from the target prediction point, the corresponding weight elements... Strictly zero, therefore the local space autocorrelation weight vector Essentially, it is a highly sparse vector. This sparsity feature has significant engineering value in large-scale computing scenarios, as it can significantly reduce storage overhead and computational load.

[0080] In this embodiment, for unfavorable edge technology scenarios with sparse distribution of large-scale sampling points, the adaptive bandwidth variable in step S2.2... The following dynamic search algorithm can be used for adjustment:

[0081] Initialize candidate bandwidth The basic distance threshold is the starting value of the adaptive bandwidth search process. This starting value is usually set according to the average or median sampling interval of known sampling points in the study area, and serves as the initial seed of the search algorithm.

[0082] Statistical spatial distance Number of valid known sampling points within the range Among them, a valid known sampling point refers to a known sampling point within the spatial distance defined by the current candidate bandwidth that has not been marked as a spatial outlier by the aforementioned spatial anomaly cleaning step and has complete and compliant soil physicochemical parameter data and cadmium-arsenic concentration label values. This refers to the number of sampling points that meet the above conditions.

[0083] Execute conditional logical branch judgment: If Less than the preset minimum spatial degrees of freedom threshold Then, according to the step size coefficient, the augmentation is performed. And re-statistics This continues until the number of valid samples within the calculation range of the decay kernel function meets the requirements of adaptive modeling, thereby avoiding spatial non-stationary calculations from falling into singular matrices or overfitting of weights. A preset minimum spatial degrees of freedom threshold is used. This refers to a pre-defined positive integer parameter that ensures that, within the local spatial gravitational field at any target prediction point, at least [the gravitational field is included]. This ensures that there are several valid known sampling points, thereby guaranteeing sufficient statistical degrees of freedom for subsequent local model parameter estimation. In this embodiment, The typical value range is 30 to 50, which is usually determined based on the number of parameters in the local model with a certain margin.

[0084] For example, in this embodiment, the above-described adaptive bandwidth dynamic search process can be formally represented as the following iterative recursive formula:

[0085]

[0086] in, Indicates the target prediction point In the The adaptive bandwidth candidate value for the next iteration, the iteration process starts from the above basic distance threshold; Indicates the first In this embodiment, the update bandwidth value during the next iteration uses a geometric increment factor of 1.2. The selection of this factor balances search efficiency and accuracy—too small a factor will lead to slow convergence of the search process, while too large a factor may lead to excessive bandwidth expansion and loss of spatial locality. and They represent the target prediction points respectively. With known sampling points The spatial coordinate vector; This represents the Euclidean distance between two points; This indicates an indicator function that outputs 1 when the conditional expression within the parentheses is true, and 0 otherwise; the iteration process first causes the number of effective sampling points in the space gravitational field to reach or exceed [a certain threshold]. The adaptive bandwidth is determined at the specified time. .

[0087] Understandably, the adaptive bandwidth dynamic search algorithm is designed to address the bandwidth selection dilemma caused by the uneven spatial distribution of sampling density within a large-scale prediction region. Using a globally uniform fixed bandwidth might cover too many sampling points in densely sampled areas, leading to an overly smooth local model, while covering too few sampling points in sparsely sampled areas might cause the local model to fall into an underdetermined state due to insufficient information. By allowing each target prediction point to dynamically adjust its dedicated bandwidth based on the local sampling density of its location, the algorithm ensures that regardless of whether the target point is located in a densely or sparsely sampled region, its local spatial gravitational field always incorporates a sufficient number of effective information sources, thereby guaranteeing the robustness and statistical reliability of the local model parameter estimation.

[0088] Figure 5 This diagram shows a comparison of the coverage performance of adaptive bandwidth dynamic search and fixed bandwidth search in dense areas in this embodiment. Figure 6 This diagram shows a comparison of the coverage performance of adaptive bandwidth dynamic search and fixed bandwidth search in sparse regions; from Figure 5 and Figure 6 It can be seen that the fixed bandwidth covers too many points in the dense area (overly smooth), and only covers 1-2 points in the sparse area (under-fixed); the adaptive bandwidth automatically expands to a sufficient number through a geometric increment of 1.2 times.

[0089] In this embodiment, before step S3, step S2 further includes a targeted constraint process based on the synergistic and antagonistic relationship between cadmium and arsenic. The purpose of this targeted constraint process is that, in actual soil environmental systems, cadmium (Cd) and arsenic (As), as two typical soil heavy metal pollutants, do not migrate and transform independently in the soil. Instead, they are jointly regulated by multiple environmental factors such as soil redox potential (Eh), pH, organic matter content, and iron and manganese oxides, exhibiting complex synergistic or antagonistic relationships. Specifically, in strongly reducing paddy soil environments (such as long-term flooded paddy fields), the soil redox potential decreases, and the bioavailability of arsenic is significantly enhanced. Meanwhile, under reducing conditions, cadmium tends to combine with sulfides to form insoluble cadmium sulfide precipitates, potentially exhibiting an antagonistic relationship where one increases while the other decreases. However, in areas with compound pollution affected by mining activities or industrial emissions, cadmium and arsenic may originate from the same pollution source simultaneously, exhibiting a synergistic relationship. Traditional machine learning models model cadmium and arsenic as two completely independent prediction targets, methodologically severing the intrinsic connection between the two elements in their geochemical mechanisms.

[0090] In this embodiment, the targeting constraint process may specifically include:

[0091] Construct a double-penalty covariance matrix and calculate the autocorrelation weight vectors of cadmium and arsenic label values ​​in the local space, respectively. The weighted Pearson correlation coefficient under constraints. The weighted Pearson correlation coefficient refers to the correlation coefficient at the target prediction point. The spatial weights calculated using the double-square kernel function within the local spatial gravitational field. The statistic, denoted by , is a weighted correlation measure between the concentrations of cadmium and arsenic, and its value ranges from a closed interval. .

[0092] For example, in this embodiment, for the target prediction point Locally weighted Pearson correlation coefficient The following formula is used to calculate:

[0093]

[0094] in, Indicates at the target prediction point Within the local spatial gravitational field, the weighted Pearson correlation coefficient between the concentrations of the two heavy metal elements cadmium and arsenic is given, where +1 indicates a perfectly positive correlation (cooperative relationship), -1 indicates a perfectly negative correlation (antagonistic relationship), and 0 indicates no linear correlation. This represents the spatial weighting coefficient calculated by the aforementioned bisquare kernel function, which ensures that when calculating local correlation, the sampling point closer to the target point contributes more. and These represent the known sampling points. The measured concentrations of cadmium and arsenic at the site; and These represent the target points respectively. The weighted arithmetic mean of cadmium and arsenic concentrations within a local spatial gravitational field, where the weighted mean is... The weights are calculated to ensure that the mean estimate itself also has spatial locality.

[0095] After calculating the weighted Pearson correlation coefficient, the weighted Pearson correlation coefficient is dynamically injected into the loss function in step S3 as a structured penalty term, forcing the geographic weighted integrated tree model to jointly evaluate the synergistic nonstationary response of the threshold to changes in cadmium and arsenic concentrations when extracting the local soil physicochemical property split threshold.

[0096] For example, in this embodiment, the targeted penalty term The following formula is used to calculate:

[0097]

[0098] in, Indicates the target prediction point The targeted penalty term, which will be subtracted from the structural gain value in the subsequent step S3, thereby constraining the model's splitting decision; The target penalty coefficient is represented by this hyperparameter, which controls the relative importance of geochemical coupling constraints in the overall objective function. Its value needs to be adjusted according to the prior knowledge of the synergistic relationship between cadmium and arsenic in the specific application scenario. This represents the absolute value of the locally weighted correlation coefficient, and its range is... ; It is a monotonically decreasing exponential mapping function—the function outputs its minimum value when the absolute value of the correlation coefficient is 1 (strong correlation). The penalty is smaller because the model only needs to capture a common driving mechanism to explain the spatial distribution of the two elements at the same time; when the absolute value of the correlation coefficient is 0 (no correlation), the function outputs a maximum value of 1, which has a larger penalty, forcing the model to be more cautious in feature selection in order to avoid making a one-sided decision that only benefits one element while ignoring the other. Figure 7 The response curve of the targeting penalty term as a function of the cadmium-arsenic local correlation intensity in this embodiment is shown. Figure 7 This demonstrates how the monotonically decreasing nature of the penalty function in this embodiment forces the model to be more cautious in selecting splitting features in weakly correlated regions.

[0099] S3: Process the initial dataset and local spatial autocorrelation weight vectors, construct and iteratively generate a geographically weighted ensemble tree model to adaptively analyze the spatial nonstationarity driving relationship between soil physicochemical properties and cadmium-arsenic distribution at a large scale.

[0100] Among them, the geographically weighted ensemble tree model refers to an ensemble learning model with spatial awareness formed by introducing local weights generated by the spatial distance decay kernel function into the Taylor expansion term of the loss function on the basis of the classic gradient boosting decision tree framework, and fully integrating the spatial location of the target prediction point as an implicit parameter into the entire construction process of the decision tree.

[0101] Spatially nonstationary driving relationships refer to the fact that the driving mechanisms of soil physicochemical parameters on cadmium-arsenic concentrations in large-scale geographic spaces are not completely identical or constant at all locations within the study area, but rather exhibit significant heterogeneity with spatial location. For example, in areas surrounding industrial and mining enterprises, changes in pH value may be the dominant factor influencing cadmium migration; while in natural background areas far from pollution sources, the geochemical composition of the parent material may be the key driving force determining the background cadmium content.

[0102] In this embodiment, the construction and iterative generation of the geographic weighted ensemble tree model specifically includes: during the node splitting process of each decision tree, extracting the sample data falling into the current node, multiplying the local spatial autocorrelation weight vector with the derivative of the loss function of each sample to calculate the geographic weighted gradient statistic, and performing deterministic feature threshold traversal based on the geographic weighted gradient statistic to maximize the spatial weighted structure gain as the splitting criterion to generate the tree topology.

[0103] The derivative of the loss function refers to the partial derivative of the global loss function with respect to the current predicted value of the model, including the first and second derivatives. The first derivative reflects the direction and magnitude of the residual between the model's predicted value and the true value, while the second derivative reflects the curvature information of the loss function at the current predicted value.

[0104] Multiplicative associative operation refers to performing element-wise scalar multiplication of the weight value corresponding to each sampling point in the local spatial autocorrelation weight vector with the derivative value of the loss function of that sampling point, thereby injecting spatial proximity information into the gradient statistics in the form of a multiplicative modulation factor. Geographically weighted gradient statistics refer to the aggregated gradient value after spatial weight modulation, including the weighted first-order gradient sum and the weighted second-order gradient sum. Spatial weighted structural gain refers to a numerical indicator calculated based on the geographically weighted gradient statistics, used to evaluate the merits of each candidate node splitting scheme. Tree topology refers to the complete hierarchical connection relationship of the decision tree from the root node to all leaf nodes, including the splitting features and splitting thresholds selected by each internal node, and the predicted values ​​output by each leaf node.

[0105] In this embodiment, the specific steps for performing deterministic feature threshold traversal based on geographic weighted gradient statistics may include:

[0106] S3.1: For each known sampling point within the current node Calculate the first derivative of the global loss function with respect to the current model predictions. With the second derivative .

[0107] Wherein, the first derivative Indicates the global loss function for the first... The first partial derivative of each sample's predicted value is the negative residual under the squared loss function. This value reflects the direction and magnitude of the deviation between the current model's prediction and the true label value; the second derivative... Indicates the global loss function for the first... The second-order partial derivative of each sample prediction is always a constant 1 under the squared loss function. This value is related to the local curvature of the loss function and plays an adaptive role in adjusting the learning rate during the optimization process.

[0108] S3.2: Combine the first and second derivatives with the local spatial autocorrelation weights. Coupled, calculate the spatially weighted first-order gradient of the left subtree. and second-order gradient and Similarly, calculate the spatially weighted first-order gradient of the right subtree. With second-order gradient and .

[0109] in, This represents the weighted sum of first-order gradients for all samples assigned to the left subtree of the current candidate splitting scheme for the target prediction point. This value reflects the spatial weighting potential of the left subtree samples for further reducing the loss function. This represents the weighted sum of second-order gradients for all samples assigned to the left subtree with respect to the target prediction point; This represents the set of sample indices that have been assigned to the left child node based on the current candidate's splitting features and splitting threshold. The spatial weights are calculated using the aforementioned bisquare kernel function. The introduction of these weights means that the gradient values ​​corresponding to the residuals of samples closer to the target prediction point are amplified to a greater extent during aggregation, while the gradient contributions of samples farther away are compressed accordingly or even approach zero.

[0110] Understandably, this weighted gradient aggregation mechanism differs fundamentally from the classic XGBoost approach of simply summing the gradients of all samples. In the classic formula, the contribution of each sample is treated equally; however, in this model, the gradient contribution of each sample is modulated by its spatial weight. Because spatial weights decay non-linearly with distance, the decision tree, in each evaluation of candidate splitting schemes, is mathematically compelled to prioritize fitting the local data distribution patterns within the neighborhood of the target point, rather than pursuing a global compromise solution that is roughly good for all spatial locations but actually not good enough anywhere.

[0111] S3.3: For a given feature splitting threshold, calculate the spatially weighted structure gain:

[0112]

[0113] in, This represents the spatial weighted structure gain that the current candidate splitting scheme can bring. The larger the value, the more valuable the splitting scheme is for fitting samples in the neighborhood of the target point. , and , Let represent the weighted first-order and second-order gradient sums of the left and right subtrees, respectively; This represents the L2 regularization penalty coefficient, which is applied over the second-order gradient and used to control the magnitude of the leaf node weights to prevent the model from overfitting—in a physical sense, Its function is similar to applying a frictional force on the optimization path, preventing the model from generating numerically unstable extreme weight values ​​in order to pursue the ultimate fit of local residuals. This represents the leaf node splitting penalty constant, which is deducted from the gain as a fixed cost of splitting. It controls the tree depth and the number of leaf nodes, achieving model structural simplicity—it only occurs when the gain exceeds a certain threshold. The split will only be executed when the time is right; otherwise, the current node will be retained as a leaf node.

[0114] It is understandable that, in embodiments that include the targeting constraint process in step S2, the spatially weighted structural gain needs to be further subtracted from the targeting penalty term. The actual formula for calculating the structure gain is:

[0115]

[0116] This formula directly incorporates the cadmium-arsenic geochemical coupling constraint into the evaluation criteria for splitting decisions, ensuring that the feature importance ranking learned by the model is consistent with prior knowledge in the field of environmental geochemistry.

[0117] Figure 8This diagram illustrates a comparison of the splitting threshold selection between global structural gain and spatially weighted structural gain in this embodiment. Figure 8 The globally optimal split point (pH≈6.5) is a compromise solution; this scheme selects pH≈5.8 in the mining area and pH≈7.2 in the background area, with higher peak values ​​and different locations; this shows how the spatially weighted gradient statistics make the splitting strategy of the decision tree vary depending on the location of the target point.

[0118] S3.4: Traverse all features and potential partition points of the soil physicochemical property parameter matrix, and select those that make the following... Deterministic node splitting is performed to maximize the feature and threshold values. Deterministic node splitting refers to the process of exhaustively searching all possible feature-threshold combinations, given spatial weight configurations and gradient statistics, to find the unique optimal splitting scheme that maximizes the global structural gain, and then executing the split accordingly. This process is completely deterministic and does not involve any random sampling or approximate estimation.

[0119] In this embodiment, when traversing potential partition points in step S3.4, a fixed-size bitmap binning acceleration operation is performed on the continuous parameters.

[0120] The technical motivation for this bitmap binning acceleration operation is that, in large-scale soil data scenarios, if gradient aggregation calculation is performed again on all samples in the current node for each candidate splitting threshold, the computational complexity will increase sharply with the increase of sample size and feature dimension, becoming the main bottleneck restricting the efficiency of model training.

[0121] For example, in this embodiment, the bitmap binning acceleration operation specifically includes:

[0122] Mapping continuously distributed soil physicochemical properties to Discretized binary fixed-length histogram bins.

[0123] in, This represents the total number of bins, which in this embodiment is typically a fixed constant of 256 or 512. This mapping process is performed as a one-time preprocessing step before model training begins, dividing the value range of each continuous environmental feature into equal-width or equal-frequency groups. The system divides samples into non-overlapping bins and replaces the continuous values ​​of each sample on each feature with the index number of its respective bin.

[0124] When calculating the spatial weighted first-order gradient sum and second-order gradient sum of the left and right subtrees, the prefix sum query is directly performed on the weighted gradient bitmap corresponding to the fixed-length histogram bins, which eliminates redundant memory access caused by sorting massive floating-point numbers while maintaining the accuracy of spatial heterogeneity parsing.

[0125] For example, in this embodiment, the weighted first-order gradient based on bitmap binning prefix and technology is used. As shown in the following formula:

[0126]

[0127] in, For target prediction points All eigenvalues ​​fall into the first The sum of the weighted first-order gradients of the samples in each bin, which is pre-computed once and stored as a lookup table before model training; This indicates that when the splitting threshold is set to the first... When the bin boundary is reached, the weighted first-order gradient summation of all samples falling into the left subtree; This represents the bin index corresponding to the candidate split threshold, obtained by traversing... From 1 to the total number of boxes This allows us to enumerate all possible splitting schemes.

[0128] Understandably, the introduction of bitmap binning prefixes and techniques reduces the computational complexity of traditional sample-by-sample sorting from... Significantly reduced to ,in The number of bins is a fixed constant. This optimization completely eliminates the CPU cache breakdown and memory bandwidth bottlenecks caused by a large number of floating-point multiplication operations when constructing high-dimensional weight matrices from the underlying algorithm architecture, laying a solid computational foundation for efficient training of the model on a scale of millions of sampling points. Figure 9 This diagram shows a double logarithmic axis comparison of the computational complexity of the bitmap binning prefix sum in this embodiment and the traditional sorting method. Figure 9 It can be seen that the time consumption of the traditional sorting method increases superlinearly with N, while the time consumption of the bitmap binning method in this embodiment is a constant horizontal line.

[0129] In this embodiment, step S3 further includes determining the output weight of the leaf node according to the following rules:

[0130] For a leaf node that has completed splitting, calculate and assign weight values ​​to the leaf node based on the spatially weighted first-order and second-order gradients of all samples within that node. .

[0131] For example, in this embodiment, the geographically weighted prediction weight of the leaf node The following formula is used to calculate:

[0132]

[0133] in, This represents the increment of the final geographically weighted prediction value output by the current leaf node when the target prediction point is routed to it. This value will be multiplied by the learning rate and then added to the prediction result of the target point. and Let represent the weighted first and second gradients of all samples falling into the leaf node with respect to the target prediction point, respectively; This is the L2 regularization coefficient, and its function is the same as described above.

[0134] In this embodiment, to address the memory overflow and high-concurrency computing blocking bottlenecks caused by large-span, large-scale spaces, step S3 employs the following customized memory scheduling mechanism:

[0135] S3.a: Using a spatial quadtree data structure, the global coordinate boundary of a large-scale prediction area is divided into multiple local sub-computation domains with a set overlap width. The spatial quadtree data structure is a classic spatial indexing structure that recursively divides the two-dimensional geographic space into four equal parts, decomposing the entire study area into several hierarchically nested spatial units.

[0136] Local sub-computation domains refer to independent computational partitions that are obtained after spatial quadtree decomposition, each covering a portion of the geographic space of the study area. Each sub-computation domain contains known sampling point data and raster pixels to be predicted within that partition.

[0137] The overlap width refers to the spatial width of the pre-defined overlap buffer zone between two adjacent local sub-computation domains at the spatial boundary. The purpose of setting up this overlap area is to provide data redundancy for subsequent boundary fusion processing.

[0138] For example, in this embodiment, the spatial quadtree decomposition and sub-computation domain partitioning can be represented by the following formula:

[0139]

[0140] in, This represents the spatial domain spanned by the entire study area in a two-dimensional geographic coordinate system. This represents the result obtained after quadtree decomposition. Local sub-computation domains The union of these sub-computational domains collectively covers the entire study area; This indicates that there is a non-empty intersection region between any two adjacent sub-computation domains; This represents the spatial width of the buffer overlap band. The setting of this width needs to be balanced between computational redundancy and boundary smoothness, and is usually taken as 1 to 2 times the adaptive bandwidth.

[0141] S3.b: The soil physicochemical property parameter matrix corresponding to the local sub-computation domain is sharded and loaded into the L3 cache of the independent computing nodes. The geographically weighted ensemble tree model construction of the specific local sub-computation domain is then executed in parallel within each independent computing node. Here, an independent computing node refers to a single working process or physical computing unit in a distributed computing cluster. Each computing node has independent processor, memory, and cache resources.

[0142] L3 cache refers to the third-level cache in a compute node processor. Its capacity is typically tens of megabytes. It is used to provide high-speed data transfer between the processor core and main memory. Loading the parameter matrix of a local sub-computation domain into the L3 cache can significantly reduce data access latency and improve computational throughput.

[0143] S3.c: For coordinate points within the overlap width, a spatial inverse distance weighted fusion function is used to deterministically fuse the predicted leaf node results of two adjacent local sub-computation domains to eliminate boundary abrupt artifacts caused by block computation. The spatial inverse distance weighted fusion function is a nonlinear weighted averaging method based on the distance from the target point to the geometric center of each sub-computation domain. This method follows the principle that sub-computation domain models with closer distances receive higher fusion weights.

[0144] Boundary abrupt change artifacts refer to the discontinuous numerical jumps that may occur at the boundaries of predicted values ​​when spatially stitching together predictions made independently by each sub-computational domain. These jumps are caused by differences in the local model parameters relied upon by adjacent sub-computational domains and appear as clearly visible seams in the final output spatial distribution map. This artifact is a purely artificial trace introduced by the computational partitioning strategy and has no relation to the actual spatial distribution patterns of heavy metals in the soil.

[0145] For example, in this embodiment, the target point is located within the overlapping area. Fusion prediction value The following formula is used to calculate:

[0146]

[0147] in, This represents the target point located within the overlapping buffer zone after spatial inverse distance weighted fusion. The final deterministic prediction value obtained; This indicates that its jurisdiction covers the target point. The set of indices of all sub-computation domains. For a point located at the boundary of two adjacent sub-computation domains, the set usually contains two elements, while for a point located at the corner of four sub-computation domains, the set may contain four elements. Indicates the sub-computational domain The predicted value output by the built-in local geographic weighted decision tree model for the target point; Subcomputation domain The geometric center coordinates of the point can be understood as the optimal applicable location of the local model carried by the sub-computation domain; This represents the Euclidean distance between the target point and the geometric center of the sub-computation domain; Represents the negative of distance The exponentiation achieves a non-linear mapping from distance to weight; The power exponent representing distance decay is typically set to 2 in this embodiment, consistent with the classic inverse distance weight interpolation method, which enables a smooth weight transition that conforms to the spatial autocorrelation law.

[0148] Understandably, the boundary fusion mechanism based on inverse distance weights described above does not simply select the output from a single sub-computation domain model for the final predicted value of any target point located within the overlapping buffer zone. Instead, it performs a spatially location-sensitive weighted average of the outputs from all relevant sub-computation domain models. When the target point is closer to the geometric center of a sub-computation domain, the model's fusion weight is greater; as the target point moves further away from the center of a sub-computation domain and approaches its boundary, the model's weight smoothly decreases. Through this spatially location-based weight gradient mechanism, the predicted value exhibits continuous and gradual spatial transition characteristics at the boundaries of sub-computation domains, thereby completely eliminating the jigsaw artifacts that may be introduced by block computation and ensuring the physical continuity and visual integrity of the large-scale prediction map.

[0149] S4: Obtain the gridded spatial coordinates of the area to be predicted and the corresponding target soil physicochemical property parameters, input them into the completed geographically weighted integrated tree model, and output and generate a large-scale cadmium-arsenic concentration distribution prediction map of the target paddy field soil.

[0150] The region to be predicted refers to the target geographic area from which the spatial distribution prediction results of cadmium-arsenic concentrations need to be generated. This area can be the same as the collection area of ​​the initial dataset, or it can be a subset or an extended area. Gridded spatial coordinates refer to the set of geographic coordinates corresponding to the center point of each raster cell after the region to be predicted is regularly rasterized according to a preset spatial resolution (e.g., 100m × 100m, 250m × 250m, or 1000m × 1000m). This gridding process discretizes the continuous geographic space into a finite set of coordinate points that can be enumerated one by one, allowing the model to output a predicted value for each raster cell. Target soil physicochemical property parameters refer to the set of parameters that correspond one-to-one with the gridded spatial coordinates and describe the soil environmental attribute characteristics at each raster cell in the region to be predicted. This parameter set has the same feature dimensions and variable types as the soil physicochemical property parameter matrix in the initial dataset.

[0151] In this embodiment, the target soil physicochemical properties parameters can be obtained in the following ways: based on existing national soil attribute databases, remote sensing inversion products (such as remote sensing estimation maps of soil organic carbon content, surface temperature inversion products, etc.), topographic parameters derived from digital elevation models (DEMs) (such as slope, aspect, topographic humidity index, etc.), and multi-source auxiliary data such as land use / land cover classification maps, complete soil physicochemical property parameter estimates are obtained at each raster cell in the area to be predicted through spatial interpolation or data fusion techniques.

[0152] The target large-scale cadmium-arsenic concentration distribution prediction map for paddy fields refers to a two-dimensional map product generated by spatially visualizing and rendering the point-by-point prediction results of a geographically weighted integrated tree model at all gridded spatial coordinates in the form of a raster layer. This prediction map uses different color levels or grayscale values ​​to represent the predicted cadmium or arsenic concentration values ​​at different raster pixels, intuitively presenting the distribution pattern and concentration gradient of heavy metals in large-scale geographic space.

[0153] In this embodiment, when inputting into the completed geographically weighted ensemble tree model for prediction, it specifically includes: for each gridded spatial coordinate point (i.e., the target prediction point) in the region to be predicted. First, based on the distance relationship between its spatial coordinates and all known sampling points, the local spatial autocorrelation weight vector of the target prediction point is calculated according to the method in step S2. Subsequently, using this weight vector As a spatial modulation parameter, the target prediction point The corresponding soil physicochemical properties are input into the geographic weighted ensemble tree model. Each decision tree within the model performs node-by-node feature discrimination and routing based on its tree topology determined during the training phase using spatial weighted structure gain. Ultimately, the target prediction point is routed to the corresponding leaf node, and the geographic weighted prediction weight value of that leaf node is obtained. The prediction weight values ​​of all leaf nodes output by all decision trees are then multiplied sequentially by the learning rate and summed to obtain the target prediction point. The final predicted values ​​of cadmium or arsenic concentrations are obtained; the predicted values ​​of all gridded spatial coordinate points are summarized and rasterized according to spatial location to generate a large-scale predicted map of cadmium-arsenic concentration distribution in paddy field soil.

[0154] It is understandable that, since the geographically weighted ensemble tree model has already deeply embedded spatial location information into the splitting strategy and leaf node weights of each decision tree through the spatial weight mechanism during the training phase, different prediction results may be obtained for target prediction points in different spatial locations even if the same soil physicochemical property parameters are input. This is the core manifestation of the model's ability to characterize the driving relationship of spatial nonstationarity.

[0155] Example 2

[0156] This embodiment discloses a large-scale cadmium-arsenic prediction system for paddy field soil based on machine learning. Specifically, this system can be integrated into an electronic device, such as a terminal or server.

[0157] The terminal can be a mobile phone, tablet computer, smart Bluetooth device, laptop computer, or personal computer; the server can be a single server or a server cluster composed of multiple servers.

[0158] In this embodiment, the machine learning-based large-scale paddy soil cadmium-arsenic prediction system can also be integrated into multiple electronic devices. For example, the machine learning-based large-scale paddy soil cadmium-arsenic prediction system can be integrated into multiple servers, and multiple servers can implement the method disclosed in Embodiment 1 of this application.

[0159] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A large-scale cadmium-arsenic prediction method for paddy field soil based on machine learning, characterized in that, The prediction method includes: S1. Obtain the initial dataset for the large-scale prediction area. The initial dataset includes the spatial coordinates of multiple known sampling points, the soil physicochemical property parameter matrix, and the corresponding cadmium-arsenic concentration label values. S2. Construct a geographic weighted non-stationary spatial distance matrix based on the spatial coordinates. Map the geographic weighted non-stationary spatial distance matrix by introducing a kernel function based on spatial distance decay. Dynamically calculate and generate a local spatial autocorrelation weight vector for any target prediction point in the large-scale prediction area. S3. Based on the initial dataset and the local spatial autocorrelation weight vector, construct a geographically weighted ensemble tree model; wherein... In the process of constructing the geographic weighted ensemble tree model, the local spatial autocorrelation weight vector and the derivative of the loss function are multiplicatively combined to calculate the geographic weighted gradient statistic, and the spatial weighted structure gain is maximized based on the geographic weighted gradient statistic as the node splitting criterion. S4. Obtain the gridded spatial coordinates of the area to be predicted and the corresponding target soil physicochemical property parameters, input them into the geographic weighted integrated tree model, and output the cadmium-arsenic concentration distribution prediction map.

2. The large-scale cadmium-arsenic prediction method for paddy field soil based on machine learning according to claim 1, characterized in that, The kernel function based on spatial distance attenuation is an adaptive bisquare spatial attenuation kernel function, which is constructed based on an adaptive bandwidth variable, wherein: The adaptive bandwidth variable is a distance parameter determined individually for each target prediction point, used to control the effective range of spatial weight decay. The adaptive bisquared spatial decay kernel function has a tight support characteristic. When the distance between the target prediction point and the known sampling point reaches or exceeds the adaptive bandwidth variable, the corresponding weight element in the local spatial autocorrelation weight vector is zero.

3. The large-scale cadmium-arsenic prediction method for paddy field soil based on machine learning according to claim 2, characterized in that, The adaptive bandwidth variable is determined through the following dynamic adjustment steps: The adaptive bandwidth variable is initialized to a base distance threshold; The distance between the statistical and target prediction points is less than the number of valid known sampling points for the current adaptive bandwidth variable; If the number of valid known sampling points is less than the preset minimum spatial degrees of freedom threshold, the current adaptive bandwidth variable is amplified and re-counted according to the preset step size coefficient, and the process is iterated until the number of valid known sampling points reaches or exceeds the preset minimum spatial degrees of freedom threshold.

4. The large-scale cadmium-arsenic prediction method for paddy field soil based on machine learning according to claim 1, characterized in that, After step S2 and before step S3, the method further includes a targeted constraint step based on the cadmium-arsenic synergistic and antagonistic relationship. For each target prediction point, the weighted Pearson correlation coefficient between the cadmium concentration label value and the arsenic concentration label value is calculated based on the local spatial autocorrelation weight vector. Calculate the targeted penalty term based on the weighted Pearson correlation coefficient; When calculating the spatial weighted structure gain in step S3, the targeting penalty term is subtracted from the spatial weighted structure gain.

5. The large-scale cadmium-arsenic prediction method for paddy field soil based on machine learning according to claim 4, characterized in that, The calculation of the weighted Pearson correlation coefficient includes: For each known sampling point within the local spatial range of the target prediction point, the first weighted arithmetic mean of the cadmium concentration label value and the second weighted arithmetic mean of the arsenic concentration label value are calculated using the weight elements in the local spatial autocorrelation weight vector as weights. Calculate the difference between the cadmium concentration label value at each known sampling point and the first weighted arithmetic mean, and the difference between the arsenic concentration label value at each known sampling point and the second weighted arithmetic mean; The weighted covariance is obtained by weighting the product of the two differences with the weighted elements. The two differences are weighted and summed using the weighted elements respectively, and the square root is taken to obtain the weighted standard deviation of cadmium concentration and the weighted standard deviation of arsenic concentration. The weighted Pearson correlation coefficient is obtained by dividing the weighted covariance by the product of the weighted standard deviation of cadmium concentration and the weighted standard deviation of arsenic concentration.

6. The large-scale cadmium-arsenic prediction method for paddy field soil based on machine learning according to claim 1, characterized in that, The derivative of the loss function includes the first derivative and the second derivative, and the steps for calculating the geographic weighted gradient statistic include: For each known sampling point within the current node, calculate the first and second derivatives of the global loss function with respect to the current model prediction value; The weight elements in the local spatial autocorrelation weight vector are multiplied and combined with the first derivative of the corresponding known sampling point to obtain the weighted first derivative of each known sampling point. The weight elements in the local spatial autocorrelation weight vector are multiplicatively combined with the second derivatives of the corresponding known sampling points to obtain the weighted second derivatives of each known sampling point. According to the candidate splitting scheme, the weighted first derivatives corresponding to the samples to be assigned to the left subtree are accumulated to obtain the weighted first gradient sum of the left subtree, and the weighted second derivatives corresponding to the samples to be assigned to the left subtree are accumulated to obtain the weighted second gradient sum of the left subtree. Similarly, calculate the weighted first-order gradient sum and the weighted second-order gradient sum of the right subtree.

7. The large-scale cadmium-arsenic prediction method for paddy field soil based on machine learning according to claim 1, characterized in that, In step S3, when traversing the candidate splitting thresholds, bitmap binning is used to accelerate the operation: Map the continuous parameters in the soil physicochemical property parameter matrix to a preset number of fixed-length histogram bins; For each target prediction point, the weighted gradient accumulation value corresponding to the sample in each fixed-length histogram bin is pre-calculated and stored; When calculating the geographic weighted gradient statistics corresponding to the candidate splitting scheme, the values ​​are obtained by performing a prefix sum query on the weighted gradient accumulation values ​​of the fixed-length histogram bins.

8. The large-scale cadmium-arsenic prediction method for paddy field soil based on machine learning according to claim 7, characterized in that, A distributed computing mechanism is used in step S3: Using a spatial quadtree data structure, the large-scale prediction region is divided into multiple local sub-computation domains, with an overlap width set between adjacent local sub-computation domains; The data corresponding to each of the local sub-computation domains is distributed to independent computing nodes, and the local geographic weighted integrated tree model is constructed in parallel within each independent computing node. Furthermore, for the target prediction point located within the overlap width range, a spatial inverse distance weighted fusion function is used to fuse the prediction results output by adjacent local sub-computation domains: the distance from the target prediction point to the geometric center of each local sub-computation domain covering the target prediction point is calculated; using the negative power of each distance as the weight, the predicted values ​​output by the geographic weighted ensemble tree model corresponding to each local sub-computation domain for the target prediction point are weighted and averaged to obtain the fused prediction value of the target prediction point.

9. The large-scale cadmium-arsenic prediction method for paddy field soil based on machine learning according to claim 1, characterized in that, The steps in step S4 for inputting data into the geographic weighted ensemble tree model for prediction include: For each target prediction point corresponding to the gridded spatial coordinates in the region to be predicted, the local spatial autocorrelation weight vector of the target prediction point is calculated based on the distance relationship between the target prediction point and each known sampling point. Using the local spatial autocorrelation weight vector as a spatial modulation parameter, the target soil physicochemical property parameters corresponding to the target prediction point are input into the geographic weighted integrated tree model; The weights of the leaf nodes output by each decision tree in the geographic weighted ensemble tree model are multiplied by the learning rate and then summed to obtain the predicted cadmium-arsenic concentration value of the target prediction point. The predicted cadmium-arsenic concentration values ​​of all target prediction points are rasterized and assembled according to their spatial location to generate the cadmium-arsenic concentration distribution prediction map.

10. A large-scale cadmium-arsenic prediction system for paddy field soil based on machine learning, characterized in that, The prediction system includes: processor; The memory stores a computer program that, when executed by a processor, implements the machine learning-based large-scale cadmium-arsenic prediction method for paddy soil as described in any one of claims 1 to 9.

Images