Dynamic monitoring method and system for soil trace elements for land remediation

Abstract

The present invention relates to the field of land consolidation technology, and discloses a method and system for dynamic monitoring of soil trace elements for land consolidation. The method comprises obtaining a soil element dataset; performing autocorrelation analysis on the soil element dataset, and performing spatial interpolation based on the obtained autocorrelation characteristics to obtain an element spatial distribution map; performing time series analysis on the soil element dataset to obtain temporal variation characteristics; performing spatiotemporal coupling based on the element spatial distribution map and the temporal variation characteristics to obtain spatiotemporal prediction model parameters; performing Bayesian optimization on the spatiotemporal prediction model parameters to obtain optimized prediction model parameters; and performing distribution prediction based on the optimized prediction model parameters, the element spatial distribution map, and the temporal variation characteristics to obtain trace element prediction results, and performing visualization. The method has the following effects: the method can improve the accuracy of dynamic monitoring of soil trace elements.

Description

Technical Field

[0001] The present invention relates to the technical field of land consolidation, and in particular to a soil trace element dynamic monitoring method and system for land consolidation. Background Art

[0002] In the fields of modern agriculture and environmental protection, land consolidation projects are crucial for improving soil quality and promoting sustainable agricultural development. With the rapid development of industrialization and urbanization, the release of large amounts of industrial waste and chemicals has caused serious soil pollution, with an increasing problem of trace element imbalance. These trace elements (such as zinc, iron, and manganese) are essential for plant growth, but either excess or deficiency can affect crop yield and quality. Therefore, establishing an effective method and system for dynamic soil trace element monitoring is crucial for scientifically guiding land consolidation efforts, optimizing fertilization programs, and protecting the ecological environment.

[0003] One existing technique for monitoring soil trace elements uses laboratory analysis. The specific steps include: first, collecting soil samples according to a specific grid pattern in the target area; second, sending the collected samples to a specialized laboratory for processing and analysis, where trace element content is determined using methods such as atomic absorption spectroscopy (AAS) and inductively coupled plasma mass spectrometry (ICP-MS); finally, mapping the soil trace element distribution based on the measurement results, and formulating appropriate soil improvement measures accordingly. This method provides accurate data support, helping to gain a deeper understanding of soil nutrient status and its changing trends.

[0004] However, the entire process from sampling to the release of final results takes at least several weeks. During this period, the soil environment will change significantly, resulting in the data obtained not accurately reflecting the true state of the current soil. Each monitoring can only capture information at a specific time point, making it difficult to track short-term fluctuations in trace elements due to seasonal changes or farming activities. Assuming a typical planting cycle of 6 months, during which time the fertilization strategy needs to be adjusted multiple times to adapt to the needs of different growth stages, traditional monitoring methods cannot respond to these needs in a timely manner. In summary, due to spatial heterogeneity and temporal dynamics, the existing methods for monitoring soil trace elements have low accuracy. Summary of the Invention

[0005] The present invention provides a soil trace element dynamic monitoring method and system for land consolidation, so as to improve the accuracy of soil trace element dynamic monitoring.

[0006] In a first aspect, in order to solve the above technical problems, the present invention provides a method for dynamic monitoring of soil trace elements for land remediation, comprising:

[0007] Obtain soil element dataset;

[0008] Performing autocorrelation analysis on the soil element data set, and performing spatial interpolation based on the obtained autocorrelation characteristics to obtain an element spatial distribution map;

[0009] Performing time series analysis on the soil element dataset to obtain temporal variation characteristics;

[0010] Performing spatiotemporal coupling according to the element spatial distribution map and the time variation characteristics to obtain spatiotemporal prediction model parameters;

[0011] Performing Bayesian optimization on the spatiotemporal prediction model parameters to obtain optimized prediction model parameters;

[0012] Distribution prediction is performed based on the optimized prediction model parameters, the element spatial distribution diagram and the time variation characteristics to obtain trace element prediction results and visualize them.

[0013] In an optional embodiment, performing autocorrelation analysis on the soil element dataset and performing spatial interpolation based on the obtained autocorrelation characteristics to obtain an element spatial distribution map includes:

[0014] Calculating the Euclidean distance between each sampling point according to the soil element dataset and constructing a distance matrix;

[0015] constructing a semivariogram based on the distance matrix;

[0016] Calculating a function value according to a preset direction angle based on the semivariogram and drawing a rose diagram;

[0017] Calculating anisotropy ratio according to the rose diagram;

[0018] When the anisotropy ratio is greater than or equal to a preset anisotropy threshold, anisotropic Kriging interpolation is performed to obtain interpolation data;

[0019] When the anisotropy ratio is less than the anisotropy threshold, ordinary Kriging interpolation is performed to obtain interpolation data;

[0020] The spatial grid is divided according to the interpolation data, and a spatial distribution diagram of the elements is drawn.

[0021] In an optional embodiment, performing time series analysis on the soil element dataset to obtain time variation characteristics includes:

[0022] sorting the soil element dataset in ascending order according to timestamps to obtain a time series;

[0023] Performing first-order difference on the time series to obtain a first-order stationary series;

[0024] Perform a stationary test on the first-order stationary sequence to obtain a stationary index;

[0025] When the stability index is lower than a preset stability threshold, a stationary time series is obtained;

[0026] When the stability index is higher than the stability threshold, performing secondary difference processing and then continuing the stability test until the stability index is lower than or equal to the stability threshold;

[0027] Perform autocorrelation calculation according to the stationary time series to obtain an autocorrelation function value;

[0028] Performing partial correlation calculation based on the stationary time series to obtain a partial correlation function value;

[0029] Determining the autoregressive term order and the sliding average term order according to the autocorrelation function value and the partial correlation function value;

[0030] Performing maximum likelihood estimation based on the autoregressive term order and the sliding average term order to obtain a time variation feature;

[0031] The time variation characteristics include autoregressive coefficients, sliding mean coefficients and seasonal parameters.

[0032] In an optional embodiment, performing spatiotemporal coupling according to the element spatial distribution diagram and the time variation characteristics to obtain spatiotemporal prediction model parameters includes:

[0033] extracting a spatial covariant component according to the element spatial distribution map;

[0034] constructing a time covariant component according to the time variation characteristics;

[0035] Multiplying the spatial covariant component and the temporal covariant component to obtain a spatiotemporal covariant function;

[0036] Performing least squares optimization based on the soil element data set and the spatiotemporal covariance function to obtain an optimized covariance function;

[0037] The optimized covariance function is substituted into the Kriging equations to solve the Kriging coefficients as the parameters of the spatiotemporal prediction model.

[0038] In an optional embodiment, performing Bayesian optimization on the spatiotemporal prediction model parameters to obtain optimized prediction model parameters includes:

[0039] Obtain historical monitoring data and generate prior distribution;

[0040] Define mean square error as the objective function and define the hyperparameter search space;

[0041] activating the spatiotemporal prediction model according to the spatiotemporal prediction model parameters;

[0042] generating candidate parameters according to the hyperparameter search space and the prior distribution;

[0043] Substituting the candidate parameters into the spatiotemporal prediction model to obtain an objective function value;

[0044] Incorporating the objective function value and the candidate parameters into the prior distribution and generating a next set of candidate parameters;

[0045] Repeat the iteration until the preset maximum number of iterations is reached or the change in the objective function value is less than the preset threshold, and the optimized prediction model parameters are obtained.

[0046] In an optional embodiment, performing distribution prediction based on the optimized prediction model parameters, the element spatial distribution map, and the time variation characteristics to obtain trace element prediction results and visualize them includes:

[0047] activating the optimized prediction model according to the optimized prediction model parameters;

[0048] Predicting the uncollected area within the target area according to the optimized prediction model to obtain an unknown prediction value;

[0049] Performing Kriging correction based on the unknown predicted value in combination with the element spatial distribution map and the time variation characteristics to obtain a trace element prediction result;

[0050] A dynamic evolution graph is generated according to the trace element prediction results according to the timestamps and visualized.

[0051] In an optional embodiment, before obtaining the soil element dataset, the method further includes:

[0052] Obtain soil texture distribution map;

[0053] Dividing the regions according to the soil texture distribution map and calculating the coefficient of variation within different regions;

[0054] When the coefficient of variation is greater than or equal to a preset coefficient of variation threshold, sampling is performed using a high-density sampling scheme;

[0055] When the coefficient of variation is less than the coefficient of variation threshold, a low-density sampling scheme is adopted for sampling.

[0056] In a second aspect, the present invention provides a soil trace element dynamic monitoring system for land remediation, comprising:

[0057] Data acquisition module, used to obtain soil element data set;

[0058] A spatial distribution module is used to perform autocorrelation analysis on the soil element data set and perform spatial interpolation based on the obtained autocorrelation characteristics to obtain an element spatial distribution map;

[0059] A time series analysis module, configured to perform time series analysis based on the soil element dataset to obtain time variation characteristics;

[0060] A spatiotemporal coupling module, configured to perform spatiotemporal coupling according to the element spatial distribution diagram and the time variation characteristics to obtain spatiotemporal prediction model parameters;

[0061] A model optimization module, configured to perform Bayesian optimization based on the spatiotemporal prediction model parameters to obtain optimized prediction model parameters;

[0062] The result output module is used to perform distribution prediction based on the optimized prediction model parameters, the element spatial distribution map and the time variation characteristics, obtain the trace element prediction results, and visualize them.

[0063] In a third aspect, the present invention also provides an electronic device comprising a processor, a memory, and a computer program stored in the memory and configured to be executed by the processor, wherein when the processor executes the computer program, the method for dynamic monitoring of soil trace elements for land reclamation as described above is implemented.

[0064] In a fourth aspect, the present invention also provides a computer-readable storage medium, which includes a stored computer program, wherein when the computer program is running, the device where the computer-readable storage medium is located is controlled to execute any one of the above-mentioned methods for dynamic monitoring of soil trace elements for land reclamation.

[0065] Compared with the prior art, the present invention has the following beneficial effects:

[0066] (1) The process of acquiring the soil element dataset ensures the integrity of the basic information collected on the distribution of soil trace elements. High-precision sampling and detection methods are used to obtain multi-dimensional soil element data, providing a high-quality data foundation for subsequent spatial analysis and temporal evolution modeling. This step not only improves the accuracy and representativeness of the original data but also lays a solid foundation for the stability of the entire prediction process.

[0067] (2) Autocorrelation analysis was performed on the soil element dataset, and spatial interpolation was performed based on the obtained autocorrelation characteristics to obtain an element spatial distribution map. By quantitatively evaluating the spatial correlation between different sampling points, a continuous element spatial distribution map was generated using the Kriging method. This method effectively compensates for the information loss caused by sparse sampling points and improves the resolution and accuracy of spatial prediction.

[0068] (3) Time series analysis was performed on the soil element dataset to obtain temporal variation characteristics. Maximum likelihood estimation was used to extract the temporal evolution patterns of each element, identifying key features such as autoregressive coefficients, sliding mean coefficients, and seasonal parameters. This step enhanced the system's ability to understand the dynamic changes of soil elements and provided temporal dimension support for building more accurate spatiotemporal models.

[0069] (4) Based on the spatial distribution map of the elements and the temporal variation characteristics, spatiotemporal coupling is performed to obtain spatiotemporal prediction model parameters. The spatial distribution characteristics are combined with the temporal evolution law to construct a spatiotemporal consistent prediction model framework, and core parameters are extracted from it. This process achieves the effective fusion of spatial and temporal information, significantly improving the model's adaptability to complex environmental changes and predictive robustness.

[0070] (5) Performing Bayesian optimization on the spatiotemporal prediction model parameters to obtain optimized prediction model parameters. Model hyperparameters are automatically tuned based on the Bayesian optimization algorithm to find the optimal parameter combination under limited computing resources. This approach not only improves model training efficiency but also enhances the stability and generalization ability of prediction results, enabling the model to better adapt to diverse soil environmental conditions.

[0071] (6) Distribution prediction is performed based on the optimized prediction model parameters, the element spatial distribution map, and the temporal variation characteristics to obtain trace element prediction results, which are then visualized. By combining the optimized model parameters with the multidimensional feature input, an accurate prediction of the future soil trace element distribution is achieved, which is then visualized in the form of a map. This step not only improves the interpretability and intuitiveness of the prediction results, but also facilitates agricultural managers to make scientific decisions quickly, thereby improving overall accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] Figure 1 This is a flow chart of a method for dynamic monitoring of soil trace elements for land consolidation provided by the first embodiment of the present invention;

[0073] Figure 2 This is a structural diagram of a soil trace element dynamic monitoring system for land consolidation provided by the second embodiment of the present invention. DETAILED DESCRIPTION

[0074] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0075] Reference Figure 1 The first embodiment of the present invention provides a method for dynamic monitoring of soil trace elements for land consolidation, comprising the following steps:

[0076] S11, obtain soil element dataset;

[0077] S12, performing autocorrelation analysis on the soil element dataset, and performing spatial interpolation based on the obtained autocorrelation characteristics to obtain an element spatial distribution map;

[0078] S13, performing time series analysis based on the soil element dataset to obtain time variation characteristics;

[0079] S14, performing spatiotemporal coupling according to the element spatial distribution map and the time variation characteristics to obtain spatiotemporal prediction model parameters;

[0080] S15, performing Bayesian optimization based on the spatiotemporal prediction model parameters to obtain optimized prediction model parameters;

[0081] S16, performing distribution prediction based on the optimized prediction model parameters, the element spatial distribution diagram and the time variation characteristics, obtaining the trace element prediction results, and visualizing them.

[0082] In step S11 , a soil element dataset is acquired.

[0083] In one embodiment, before obtaining the soil element dataset, the method further includes:

[0084] Obtain soil texture distribution map;

[0085] Dividing the regions according to the soil texture distribution map and calculating the coefficient of variation within different regions;

[0086] When the coefficient of variation is greater than or equal to a preset coefficient of variation threshold, sampling is performed using a high-density sampling scheme;

[0087] When the coefficient of variation is less than the coefficient of variation threshold, a low-density sampling scheme is adopted for sampling.

[0088] It's worth noting that the soil texture distribution map is a thematic map that reflects the spatial distribution of different particle size ratios within the soil. Its core function is to provide a scientific basis for land consolidation and agricultural management. The map uses color or symbol demarcation to visually display the distribution characteristics of texture types such as sand, loam, and clay, helping to identify differences in soil physical properties.

[0089] It is worth noting that after dividing the regions into 1km grids based on the soil texture distribution map, the coefficient of variation (CV) for each region must first be calculated. The CV is calculated by first collecting soil texture data (e.g., the proportion of sand, silt, and clay) from all sampling points within the region, finding the average of these data, then calculating the difference between each data point and the average (i.e., the standard deviation). Finally, the CV is divided by the average to obtain the CV. This coefficient reflects the spatial variability of soil texture within the region—larger values ​​indicate more significant texture differences. The preset CV thresholds are: if the CV is ≥ 0.5 (e.g., for high variability of certain heavy metals), it indicates strong spatial heterogeneity in texture within the region and requires a high-density sampling strategy, such as denser sampling points (e.g., reducing the grid spacing to less than 50 meters) to capture detailed local variations. If the CV is < 0.5 (e.g., for more stable indicators such as pH), a low-density sampling strategy is adopted, such as increasing the grid spacing to more than 100 meters and reducing the number of sampling points to improve efficiency. For example, in a study of purple soil hilly areas, the coefficient of variation of available phosphorus was as high as 106%, requiring high-density sampling; while the coefficient of variation of pH was only 12.9%, so a low-density scheme was suitable.

[0090] It is worth noting that the process of acquiring a soil element dataset requires a combination of field sampling and sensor monitoring techniques. First, before sampling, appropriate sensors, including soil moisture, conductivity, or pH sensors, must be selected based on the target elements, such as heavy metals, pH, or organic matter. These sensors must be installed according to specifications at representative locations within the soil, ensuring they are stable and protected from external interference, such as foot traffic or rain. After installation, sensors must be regularly calibrated, verifying measurement accuracy with standard solutions or soil samples with known parameters to eliminate errors. During the data collection phase, raw sensor data, such as soil moisture content and conductivity, is collected in real time using data loggers or wireless transmission equipment. Laboratory analysis is then used to supplement the precise determination of key elements. Data storage utilizes a structured format. Short-term data is recorded in CSV or Excel spreadsheets for preliminary analysis, while long-term or large-scale data is stored in an SQL database for efficient querying and dynamic updating. The entire process adheres to standardized sampling methods, recording metadata such as the latitude, longitude, depth, and sampling time of sampling points to ensure data traceability and temporal consistency.

[0091] In step S12, an autocorrelation analysis is performed on the soil element data set, and spatial interpolation is performed based on the obtained autocorrelation characteristics to obtain an element spatial distribution map.

[0092] In one embodiment, the Euclidean distance between each sampling point is calculated based on the soil element dataset, and a distance matrix is ​​constructed;

[0093] constructing a semivariogram based on the distance matrix;

[0094] Calculating a function value according to a preset direction angle based on the semivariogram and drawing a rose diagram;

[0095] Calculating anisotropy ratio according to the rose diagram;

[0096] When the anisotropy ratio is greater than or equal to a preset anisotropy threshold, anisotropic Kriging interpolation is performed to obtain interpolation data;

[0097] When the anisotropy ratio is less than the anisotropy threshold, ordinary Kriging interpolation is performed to obtain interpolation data;

[0098] The spatial grid is divided according to the interpolation data, and a spatial distribution diagram of the elements is drawn.

[0099] It's worth noting that the distance matrix is ​​constructed by calculating the Euclidean distance between all sampling points, forming a symmetric matrix in which each element represents the straight-line distance between two points. For example, if there are 10 sampling points, the matrix contains 10 rows and 10 columns, where the value in row i, column j represents the distance between points i and j. (Due to the symmetry of distance, the lower left triangle and upper right triangle of the matrix are identical.) To construct the semivariogram, all sampling points are first paired according to the distance matrix. The distance between each pair of points and the corresponding soil element content difference (for example, the difference in heavy metal concentration between one point and another) are calculated. These pairs are then grouped into distance ranges (e.g., 0-10 meters, 10-20 meters, etc.), and the squared difference values ​​for all pairs within each group are averaged. This average value is then divided by two to obtain the semivariogram for the corresponding distance range. For example, if a group contains five pairs of points and the squared differences are 4, 9, 16, 25, and 36, respectively, then the mean is (4 + 9 + 16 + 25 + 36) / 5 = 18, and the semivariance is 18 / 2 = 9. Finally, the semivariance values ​​corresponding to each distance range are plotted as a curve, namely the semivariogram, which is used to characterize the spatial autocorrelation of soil elements.

[0100] It's worth noting that in semivariogram analysis, the semivariogram values ​​are first calculated for each direction according to the preset directional angles (0°, 45°, 90°, and 135°). Specifically, all pairs of sampling points are grouped according to the corresponding directional angle range. The distance between each pair and the squared difference in elemental content are then averaged to obtain the semivariograms for each direction. Subsequently, the semivariogram values ​​for these directions are plotted in polar coordinates as a rose diagram, visually demonstrating the differences in spatial variation across different directions. Each petal in the rose diagram represents the strength of the semivariogram in a particular direction, with longer petals indicating weaker spatial correlation in that direction. The rose diagram can be used to further calculate the anisotropy ratio: the ratio of the range (i.e., the distance where the semivariogram reaches a plateau) in the primary direction (the direction of maximum semivariance) to the range in the secondary direction (the perpendicular direction). For example, if the range in the primary direction is 100 meters and the range in the secondary direction is 50 meters, the ratio is 2. When the ratio is greater than or equal to the preset anisotropy threshold (2), it indicates that the data is significantly anisotropic and anisotropic kriging interpolation is required to generate more accurate interpolation results by adjusting the weight parameters in different directions. If the ratio is less than the anisotropy threshold (such as 1.5), it indicates that the spatial variation is close to isotropy and ordinary kriging interpolation can be used directly.

[0101] It is worth noting that when soil element data exhibit significant spatial variation in different directions (e.g., when variation in one direction is much greater than in the vertical direction), anisotropic kriging interpolation is required. This involves geometrically transforming the original spatial coordinates: compressing or stretching the coordinate axes along the principal directions to make the variation characteristics in different directions more consistent, thus converting the model to an isotropic one. On this basis, the computational logic of ordinary kriging interpolation is used, but the weights are calculated based on the adjusted coordinate distances. Finally, the interpolated results are projected back to the original coordinate system to generate a spatial distribution map.

[0102] It's worth noting that ordinary kriging interpolation assumes that the spatial variation of soil elements is uniform in all directions (i.e., isotropic). The process first calculates the semivariogram based on the distance matrix of all sampling points. Point pairs are grouped by distance ranges, and the squared attribute differences between pairs within each group are averaged and divided by two to obtain the semivariogram. Subsequently, a semivariogram model (such as a spherical model) is fitted to determine the range and nugget value of spatial correlation. Next, the semivariogram is calculated for the target point relative to all known points. A linear system of weighting coefficients is constructed, and the optimal weights are determined by minimizing the variance of the estimation error. Weights are assigned based on the principle of "greatest contribution from neighboring points," i.e., closer points receive higher weights. Finally, the attribute values ​​of the known points are weighted and summed to obtain the estimated value for the target point. This method eliminates directional differences and directly interpolates based on spatial distance and semivariogram relationships. It is suitable for datasets with nearly isotropic variation.

[0103] It's worth noting that ordinary kriging interpolation weights are assigned based on the semivariogram model and the spatial relationship between the known points and the target point. Specifically, the optimal weights are determined by constructing a system of linear equations. First, the semivariogram is calculated for each pair of points. Then, the semivariograms of all known points are combined with those of the target point to form a matrix. Constraints are added (e.g., the sum of the weights must be 1), forming a linear system consisting of multiple equations. Solving this system yields a weight for each known point: points closer to the target point have smaller semivariograms and higher weights. Conversely, points farther away have larger semivariograms and lower weights. For example, if three known points are located 10, 30, and 50 meters from the target point, respectively, and their corresponding semivariograms are 2, 4, and 6, respectively, the weights are proportionally assigned to 0.5, 0.3, and 0.2. Finally, the attribute value of the target point is calculated through a weighted sum (e.g., 0.5 × A + 0.3 × B + 0.2 × C). This process ensures that the weights reflect spatial proximity (closer neighbors contribute more) while also ensuring unbiased estimation.

[0104] In step S13, a time series analysis is performed based on the soil element dataset to obtain time variation characteristics.

[0105] In one embodiment, the soil element dataset is sorted in ascending order according to timestamps to obtain a time series;

[0106] Performing first-order difference on the time series to obtain a first-order stationary series;

[0107] Perform a stationary test on the first-order stationary sequence to obtain a stationary index;

[0108] When the stability index is lower than a preset stability threshold, a stationary time series is obtained;

[0109] When the stability index is higher than the stability threshold, performing secondary difference processing and then continuing the stability test until the stability index is lower than or equal to the stability threshold;

[0110] Perform autocorrelation calculation according to the stationary time series to obtain an autocorrelation function value;

[0111] Performing partial correlation calculation based on the stationary time series to obtain a partial correlation function value;

[0112] Determining the autoregressive term order and the sliding average term order according to the autocorrelation function value and the partial correlation function value;

[0113] Performing maximum likelihood estimation based on the autoregressive term order and the sliding average term order to obtain a time variation feature;

[0114] The time variation characteristics include autoregressive coefficients, sliding mean coefficients and seasonal parameters.

[0115] It's worth noting that after sorting the soil element dataset in ascending timestamp order, the observations at two adjacent time points are first subtracted from each other to produce a first-order difference sequence (for example, if the original data is [10, 12, 15, 14], the first-order difference is [2, 3, -1]). The stability of this sequence is then assessed through a stationarity test, typically using a unit root test (such as the ADF test). For example, an ADF test outputs a significance level (for example, p-value = 0.02). If the p-value is less than 0.05, the sequence is considered stationary; conversely, if the p-value is greater than 0.05, it indicates that trend or seasonality still exists in the sequence, necessitating a second differencing step. Secondary differencing involves repeating the subtraction operation on the first-order differenced sequence (for example, if the first-order difference is [2, 3, -1], the second-order difference is [1, -4]), generating a new differenced sequence for which stationarity testing is repeated. This process is repeated until the statistical p-value is less than 0.05, and finally a time series that meets the stationary requirements is obtained for subsequent modeling and analysis.

[0116] It's worth noting that the process for calculating the autocorrelation and partial correlation function values ​​from a stationary time series, and determining the order of the autoregressive term and the moving average term, is as follows: First, when calculating the autocorrelation function (ACF), correlation analysis is performed between the time series and its own values ​​at different lags. For example, for lag 1, the correlation coefficient between the current value and the previous value is calculated; for lag 2, the correlation coefficient between the current value and the previous two values ​​is calculated, and so on. This process calculates the ratio of the covariance to the variance to obtain the autocorrelation coefficient. The ACF is then plotted to observe the decay trend of the correlation coefficient at different lags. If the ACF value rapidly approaches zero (i.e., truncated) after a certain lag, it corresponds to the order of the moving average term (MA); if the ACF value decays slowly (i.e., tailing), it corresponds to the order of the autoregressive term (AR).

[0117] It is worth noting that when calculating the partial correlation function (PACF), the influence of intermediate lags must be eliminated, focusing on the "net" correlation between a certain lag and the current time point. For example, when calculating the partial correlation of lag 2, the influence of lag 1 must be controlled first, and the independent contribution of lag 2 must be isolated through stepwise regression. Specifically, a regression model is established using the value at the current time point and the values ​​of lag 1 and lag 2, respectively, and the regression coefficient of lag 2 is taken as the partial correlation value. This process needs to be performed step by step to ensure that only the independent influence of the target lag is retained at each step. After plotting the PACF graph, if the partial correlation value is truncated after a certain lag order, it corresponds to the order of the autoregressive term (AR); if the partial correlation value is tailing, it corresponds to the order of the moving average term (MA).

[0118] It's worth noting that the AR and MA orders are determined by combining the characteristics of the ACF and PACF plots. For example, if the ACF is truncated after lag q and the PACF is tailing, the MA(q) model is selected; if the PACF is truncated after lag p and the ACF is tailing, the AR(p) model is selected; if both are tailing but decay rapidly, the ARMA(p,q) model is selected. Specifically, the order is preliminarily determined by observing the number of lag points in the ACF and PACF plots that are significantly non-zero (e.g., points outside the 95% confidence interval), and then further optimized by combining information criteria (e.g., AIC). For example, if the ACF is within the confidence interval after lag 2 and the PACF is truncated after lag 1, the AR(1) model is selected; if both the ACF and PACF are tailing but decay rapidly, the ARMA(1,1) model is selected. This process requires a comprehensive judgment based on statistical tests and actual data characteristics.

[0119] It is worth noting that the process of calculating temporal variation characteristics using maximum likelihood estimation (MLE) based on the determined order of the autoregressive term (p) and the order of the moving average term (q) is as follows: First, based on the stationary time series data, assuming the error term follows a normal distribution, a likelihood function is constructed. This function measures the probability of the observed data occurring under the model parameters. Subsequently, an iterative optimization algorithm (such as the steepest descent method) is used to find the parameter values ​​that maximize the likelihood function. Specifically, the model parameters include the autoregressive coefficient (corresponding to the weights of the p lagged terms), the moving average coefficient (corresponding to the weights of the q lagged error terms), and seasonal parameters (such as the order of seasonal differencing or the coefficients of the seasonal AR / MA term). For example, if the data exhibit annual seasonality (such as monthly data), a 12th-order seasonal differencing or seasonal autoregressive term (SAR) is introduced into the model, and its parameters are estimated in a similar manner. During the calculation process, the model gradually adjusts its parameters to minimize the residual error between the predicted and observed values, thereby maximizing the likelihood function. Finally, the output is the autoregressive coefficient (such as the contribution of the first-order lag is 0.8), the moving average coefficient (such as the weight of the first-order lag error term is -0.3), and the seasonal parameter (such as the coefficient of the 12th-order seasonal difference is 0.5). These parameters together describe the dynamic changes of the time series, including trend dependence, short-term fluctuation effects and seasonal patterns.

[0120] In step S14, spatiotemporal coupling is performed based on the element spatial distribution diagram and the time variation characteristics to obtain spatiotemporal prediction model parameters.

[0121] In one embodiment, a spatial covariant component is extracted based on the element spatial distribution diagram;

[0122] constructing a time covariant component according to the time variation characteristics;

[0123] Multiplying the spatial covariant component and the temporal covariant component to obtain a spatiotemporal covariant function;

[0124] Performing least squares optimization based on the soil element data set and the spatiotemporal covariance function to obtain an optimized covariance function;

[0125] The optimized covariance function is substituted into the Kriging equations to solve the Kriging coefficients as the parameters of the spatiotemporal prediction model.

[0126] In one embodiment, the process for extracting spatial covariance components based on the spatial distribution map of soil elements is as follows: First, based on the spatial distribution data of soil elements, the covariance or semivariogram is calculated between different spatial locations to describe the spatial correlation between data points. For example, if the closer the two locations are, the more similar their attribute values ​​are, then the spatial covariance component will be higher within short distances and gradually decay with increasing distance. If there are directional differences (such as anisotropy), the covariance characteristics in different directions must be calculated separately and adjusted to an isotropic model through geometric transformation. Next, the temporal covariance component is constructed based on the temporal variation characteristics: the covariance function of the time series is constructed using the known autoregressive coefficient (e.g., the contribution of the first lag is 0.8), the moving average coefficient (e.g., the weight of the first lag error term is -0.3), and the seasonal parameter (e.g., the coefficient of the 12th-order seasonal difference is 0.5). For example, if the data have annual periodicity, the temporal covariance component will contain a periodically decaying fluctuation pattern, describing the correlation between adjacent time points.

[0127] It is worth noting that the spatial and temporal covariance components are multiplied to generate a spatiotemporal covariance function, which incorporates the effects of spatial proximity and temporal continuity. For example, if a point is 50 meters away from a reference point and one month apart, its spatiotemporal covariance is the product of the corresponding spatial covariance value (e.g., 0.6) and the temporal covariance value (e.g., 0.4) (0.24). This function must be non-negatively definite (e.g., the covariance matrix must be semi-positive definite) to ensure stability in subsequent calculations. Next, the parameters of the spatiotemporal covariance function are optimized using the least squares method: the sum of squared residuals between the actual soil element observations and the model predictions is used as the objective function, and the parameters of the covariance function (e.g., spatial extent a, temporal decay b, and sill value σ²) are adjusted. For example, the initial parameters are set to a = 100 meters, b = 30 days, and σ² = 5. By iteratively adjusting a to 80 meters, b to 25 days, and σ² to 4.5, the sum of squared residuals is minimized.

[0128] It is worth noting that the optimized spatiotemporal covariance function is substituted into the kriging equations to obtain the kriging coefficients, which serve as parameters for the spatiotemporal prediction model. The specific steps include constructing a matrix containing the spatiotemporal covariance values ​​of all known points and the target point, adding unbiasedness constraints (such as the weights summing to 1), and forming a system of linear equations. The weight coefficients are then solved using matrix inversion or iterative methods. For example, if the spatiotemporal covariance values ​​of three known points are 0.8, 0.5, and 0.3, respectively, and the constraint is that the weights sum to 1, weights 0.4, 0.35, and 0.25 are assigned. Finally, the predicted value for the target point is calculated using a weighted sum (e.g., 0.4 × A + 0.35 × B + 0.25 × C). This process also outputs the uncertainty of the predicted value (e.g., the prediction variance), thus constructing a complete spatiotemporal prediction model.

[0129] In step S15, Bayesian optimization is performed based on the spatiotemporal prediction model parameters to obtain optimized prediction model parameters.

[0130] In one embodiment, historical monitoring data is obtained and a prior distribution is generated;

[0131] Define mean square error as the objective function and define the hyperparameter search space;

[0132] activating the spatiotemporal prediction model according to the spatiotemporal prediction model parameters;

[0133] generating candidate parameters according to the hyperparameter search space and the prior distribution;

[0134] Substituting the candidate parameters into the spatiotemporal prediction model to obtain an objective function value;

[0135] Incorporating the objective function value and the candidate parameters into the prior distribution and generating a next set of candidate parameters;

[0136] Repeat the iteration until the preset maximum number of iterations is reached or the change in the objective function value is less than the preset threshold, and the optimized prediction model parameters are obtained.

[0137] It is worth noting that the specific steps for generating a prior distribution are as follows: First, key statistical features are extracted from historical monitoring data, including calculation of the data's central tendency (e.g., arithmetic mean), dispersion (e.g., sample standard deviation), and distribution shape (e.g., skewness and kurtosis). Based on these statistical features, an appropriate distribution type is selected. For example, if the data is symmetrical and lacks extreme values, a normal distribution is used; if the data is right-skewed and has a lower bound (e.g., non-negative concentration), a lognormal or gamma distribution is used. For multi-parameter scenarios (e.g., model hyperparameters), a separate distribution must be designed for each parameter. Continuous parameters (e.g., learning rate) can be assigned a uniform distribution (covering a reasonable range, such as 0.001 to 0.1) or a truncated normal distribution (centered around a typical value with upper and lower bounds). For discrete parameters (e.g., the number of hidden layer nodes), a discrete probability mass function is defined (e.g., an equiprobable uniform distribution or a decreasing probability distribution, e.g., the probability increases with a smaller number of nodes). When historical data is insufficient or the distribution shape is unclear, a broad uniform distribution or a normal distribution with large variance is used as an uninformative prior to reduce the influence of subjective assumptions. Ultimately, the generation of the prior distribution is completed through empirical Bayesian methods (such as using historical data to fit the marginal likelihood function to estimate distribution parameters) or manual settings (directly specifying the distribution type and parameter values), ensuring that it can reflect the statistical laws of historical data and adapt to the needs of subsequent Bayesian optimization.

[0138] It is worth noting that the process of performing Bayesian optimization on the spatiotemporal prediction model parameters to obtain the optimized prediction model parameters is as follows: First, a prior distribution is generated based on historical monitoring data, that is, an assumed initial distribution of the model parameters (for example, if historical data indicates that the learning rate fluctuates between 0.001 and 0.1, the prior distribution is set to a uniform or normal distribution). Next, the objective function is defined as the mean squared error (MSE), which is the average squared difference between the predicted value and the actual observed value, and a hyperparameter search space is set (for example, the learning rate range is 0.001 to 0.1, the number of hidden layer nodes ranges from 10 to 100, and the batch size ranges from 16 to 128). After activating the spatiotemporal prediction model, the prior distribution and search space are used to generate candidate parameters. For example, the parameter combination of the objective function (such as a learning rate of 0.02, a number of hidden layer nodes of 50, and a batch size of 32) is selected through random sampling or based on an acquisition function (such as the expected improvement in EI).

[0139] It is worth noting that candidate parameters are substituted into the model for training and prediction, and the corresponding objective function value (e.g., MSE = 0.05) is calculated. This value, along with the candidate parameters, is then incorporated into the prior distribution to form a new posterior distribution (e.g., by fitting the uncertainty of the objective function using a Gaussian process). Subsequently, based on the updated posterior distribution and acquisition function, the next set of candidate parameters is generated (e.g., the learning rate is adjusted to 0.018 and the number of hidden layer nodes is adjusted to 60), and this process is repeated. Iterations terminate when a preset maximum number of iterations (e.g., 50) is reached or when the change in the objective function value falls below a preset threshold (e.g., the MSE change over three consecutive iterations is less than 0.001). Finally, the optimal parameter combination (e.g., a learning rate of 0.015, 45 hidden layer nodes, and a batch size of 64) is output as the optimized prediction model parameters. This process achieves efficient optimization by dynamically adjusting the parameter search direction, balancing exploration (trying new parameter regions) and exploitation (optimizing known effective parameters), avoiding blind search and ultimately finding the parameter configuration that minimizes prediction error.

[0140] In step S16, distribution prediction is performed based on the optimized prediction model parameters, the element spatial distribution diagram, and the time variation characteristics to obtain trace element prediction results, and visualize them.

[0141] In one embodiment, the optimized prediction model is activated according to the optimized prediction model parameters;

[0142] Predicting the uncollected area within the target area according to the optimized prediction model to obtain an unknown prediction value;

[0143] Performing Kriging correction based on the unknown predicted value in combination with the element spatial distribution map and the time variation characteristics to obtain a trace element prediction result;

[0144] A dynamic evolution graph is generated according to the trace element prediction results according to the timestamps and visualized.

[0145] It is worth noting that the process of activating the optimized prediction model parameters is as follows: First, the optimal parameters obtained through Bayesian optimization (such as a learning rate of 0.015, 45 hidden layer nodes, and a batch size of 64) are loaded into the spatiotemporal prediction model to activate the model's prediction function. Subsequently, the model integrates the element concentration data of known spatial points with time series characteristics to predict unsampled areas within the target area and infer the concentration values ​​of unknown points. For example, the model combines the spatial distribution pattern and temporal trend of historical monthly sampling data for a region and calculates the prediction results using a weighted spatiotemporal covariance function. Next, the prediction value is further optimized based on kriging correction: the spatial covariance component of the unsampled points is calculated using the element spatial distribution map, the temporal covariance component is constructed based on the temporal variation characteristics, the spatiotemporal covariance function is substituted into the kriging equation system, and the kriging coefficient is solved to weight the known point data. For example, if the spatiotemporal covariance values ​​of three known points are 0.8, 0.5, and 0.3, respectively, and the constraint weights sum to 1, weights 0.4, 0.35, and 0.25 are assigned. The final weighted sum is (0.4 × A + 0.35 × B + 0.25 × C), where A, B, and C are the trace element values ​​at the known points. This yields a revised trace element forecast, and the forecast variance is output to quantify uncertainty. Finally, a dynamic evolution graph is generated and visualized based on the revised forecast values ​​by timestamp. Forecast results at different times (e.g., monthly data from January 2023 to May 2025) are arranged along the timeline, and the spatial distribution changes are displayed using heat maps or contour maps. For example, using ArcGIS or Python's Plotly library, predicted values ​​can be mapped to a color gradient (red indicates high concentration, blue indicates low concentration), overlaid with a geographic basemap, and animated frame by frame to demonstrate dynamic evolution over time. Combined with timeline controls and interactive zoom capabilities, users can observe concentration trends in specific areas or time periods, ultimately creating intuitive spatiotemporal dynamic visualizations. The entire process combines parameter tuning with Bayesian optimization, the construction of spatiotemporal covariance functions, local optimization with kriging corrections, and dynamic visualization techniques to ensure comprehensive improvements in spatial accuracy, temporal continuity, and interpretability of prediction results.

[0146] In summary, the present invention discloses a method for dynamic monitoring of soil trace elements for land consolidation, which can improve the accuracy of dynamic monitoring of soil trace elements.

[0147] Reference Figure 2 The second embodiment of the present invention provides a soil trace element dynamic monitoring system for land consolidation, comprising:

[0148] Data acquisition module, used to obtain soil element data set;

[0149] A spatial distribution module is used to perform autocorrelation analysis on the soil element data set and perform spatial interpolation based on the obtained autocorrelation characteristics to obtain an element spatial distribution map;

[0150] A time series analysis module, configured to perform time series analysis based on the soil element dataset to obtain time variation characteristics;

[0151] A spatiotemporal coupling module, configured to perform spatiotemporal coupling according to the element spatial distribution diagram and the time variation characteristics to obtain spatiotemporal prediction model parameters;

[0152] A model optimization module, configured to perform Bayesian optimization based on the spatiotemporal prediction model parameters to obtain optimized prediction model parameters;

[0153] The result output module is used to perform distribution prediction based on the optimized prediction model parameters, the element spatial distribution map and the time variation characteristics, obtain the trace element prediction results, and visualize them.

[0154] It should be noted that the soil trace element dynamic monitoring system for land reclamation provided in an embodiment of the present invention is used to execute all the process steps of the soil trace element dynamic monitoring method for land reclamation in the above embodiment. The working principles and beneficial effects of the two correspond one to one, so they will not be repeated here.

[0155] An embodiment of the present invention further provides an electronic device. The electronic device includes: a processor, a memory, and a computer program stored in the memory and executable on the processor, such as a data acquisition program. When the processor executes the computer program, the steps of the above-mentioned embodiments of the method for dynamic monitoring of soil trace elements for land remediation are implemented, such as Figure 1 Alternatively, when the processor executes the computer program, the functions of the modules / units in the above-mentioned device embodiments are realized, such as the data acquisition module.

[0156] Exemplarily, the computer program may be divided into one or more modules / units, which are stored in the memory and executed by the processor to implement the present invention. The one or more modules / units may be a series of computer program instruction segments capable of implementing specific functions, and the instruction segments are used to describe the execution process of the computer program in the electronic device.

[0157] The electronic device may be a computing device such as a desktop computer, notebook, PDA, or smart tablet. The electronic device may include, but is not limited to, a processor and memory. Those skilled in the art will appreciate that the aforementioned components are merely examples of electronic devices and do not constitute a limitation of the electronic device. The electronic device may include more or fewer components than those described above, or a combination of certain components, or different components. For example, the electronic device may also include input / output devices, network access devices, buses, and the like.

[0158] The processor may be a central processing unit (CPU), other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA), other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. A general-purpose processor may be a microprocessor or any conventional processor. The processor is the control center of the electronic device and connects various parts of the entire electronic device using various interfaces and lines.

[0159] The memory can be used to store the computer programs and / or modules. The processor implements the various functions of the electronic device by running or executing the computer programs and / or modules stored in the memory and accessing the data stored in the memory. The memory may primarily include a program storage area and a data storage area. The program storage area may store an operating system and at least one application required for a function (such as a sound playback function or an image playback function); the data storage area may store data generated based on the use of the mobile phone (such as audio data, a phone book, etc.). Furthermore, the memory may include high-speed random access memory and non-volatile memory, such as a hard disk, internal memory, a plug-in hard disk, a smart media card (SMC), a secure digital (SD) card, a flash card, at least one disk storage device, a flash memory device, or other volatile solid-state storage device.

[0160] If the module / unit integrated into the electronic device is implemented as a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the present invention can implement all or part of the process steps in the above-mentioned method embodiments by using a computer program to instruct the relevant hardware. The computer program can be stored in a computer-readable storage medium. When executed by a processor, the computer program can implement the steps of each of the above-mentioned method embodiments. The computer program includes computer program code, which can be in source code form, object code form, executable file, or some intermediate form. The computer-readable medium can include: any entity or device capable of carrying the computer program code, recording medium, USB flash drive, mobile hard drive, magnetic disk, optical disk, computer memory, read-only memory (ROM), random access memory (RAM), electric carrier signal, telecommunication signal, and software distribution medium. It should be noted that the content of the computer-readable medium can be appropriately increased or decreased based on the requirements of legislation and patent practice in a jurisdiction. For example, in some jurisdictions, based on legislation and patent practice, computer-readable media does not include electric carrier signals and telecommunication signals.

[0161] It should be noted that the device embodiments described above are merely illustrative, wherein the units described as separate components may or may not be physically separated, and the components displayed as units may or may not be physical units, that is, they may be located in one place, or they may be distributed across multiple network units. Some or all of the modules may be selected according to actual needs to achieve the purpose of the present embodiment. In addition, in the drawings of the device embodiments provided by the present invention, the connection relationship between the modules indicates that there is a communication connection between them, which may be specifically implemented as one or more communication buses or signal lines. A person of ordinary skill in the art can understand and implement the present invention without inventive effort.

[0162] The specific embodiments described above further illustrate the objectives, technical solutions, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the scope of protection of the present invention. In particular, it should be noted that 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 for those skilled in the art.

Claims

1. A method for dynamic monitoring of soil trace elements for land consolidation, characterized in that: include: Obtain soil element dataset; Performing autocorrelation analysis on the soil element data set, and performing spatial interpolation based on the obtained autocorrelation characteristics to obtain an element spatial distribution map; Performing time series analysis on the soil element dataset to obtain temporal variation characteristics; Performing spatiotemporal coupling according to the element spatial distribution map and the time variation characteristics to obtain spatiotemporal prediction model parameters; Performing Bayesian optimization on the spatiotemporal prediction model parameters to obtain optimized prediction model parameters; Performing distribution prediction based on the optimized prediction model parameters, the element spatial distribution map, and the time variation characteristics to obtain trace element prediction results and visualize them; The step of performing autocorrelation analysis on the soil element dataset and performing spatial interpolation based on the obtained autocorrelation characteristics to obtain an element spatial distribution map includes: Calculating the Euclidean distance between each sampling point according to the soil element dataset and constructing a distance matrix; constructing a semivariogram based on the distance matrix; Calculating a function value according to a preset direction angle based on the semivariogram and drawing a rose diagram; Calculating anisotropy ratio according to the rose diagram; When the anisotropy ratio is greater than or equal to a preset anisotropy threshold, anisotropic Kriging interpolation is performed to obtain interpolation data; When the anisotropy ratio is less than the anisotropy threshold, ordinary Kriging interpolation is performed to obtain interpolation data; Performing spatial grid division according to the interpolation data and drawing an element spatial distribution map; The time series analysis performed on the soil element dataset to obtain time variation characteristics includes: sorting the soil element dataset in ascending order according to timestamps to obtain a time series; Performing first-order difference on the time series to obtain a first-order stationary series; Perform a stationary test on the first-order stationary sequence to obtain a stationary index; When the stability index is lower than a preset stability threshold, a stationary time series is obtained; When the stability index is higher than the stability threshold, performing secondary difference processing and then continuing the stability test until the stability index is lower than or equal to the stability threshold; Perform autocorrelation calculation according to the stationary time series to obtain an autocorrelation function value; Performing partial correlation calculation based on the stationary time series to obtain a partial correlation function value; Determining the autoregressive term order and the sliding average term order according to the autocorrelation function value and the partial correlation function value; Performing maximum likelihood estimation based on the autoregressive term order and the sliding average term order to obtain a time variation feature; The time variation characteristics include autoregressive coefficients, sliding mean coefficients and seasonal parameters; The performing spatiotemporal coupling according to the element spatial distribution diagram and the time variation characteristics to obtain spatiotemporal prediction model parameters includes: extracting a spatial covariant component according to the element spatial distribution map; constructing a time covariant component according to the time variation characteristics; Multiplying the spatial covariant component and the temporal covariant component to obtain a spatiotemporal covariant function; Performing least squares optimization based on the soil element data set and the spatiotemporal covariance function to obtain an optimized covariance function; Substituting the optimized covariance function into the Kriging equations to solve the Kriging coefficients as parameters of the spatiotemporal prediction model; The Bayesian optimization is performed on the spatiotemporal prediction model parameters to obtain optimized prediction model parameters, including: Obtain historical monitoring data and generate prior distribution; Define mean square error as the objective function and define the hyperparameter search space; activating the spatiotemporal prediction model according to the spatiotemporal prediction model parameters; generating candidate parameters according to the hyperparameter search space and the prior distribution; Substituting the candidate parameters into the spatiotemporal prediction model to obtain an objective function value; Incorporating the objective function value and the candidate parameters into the prior distribution and generating a next set of candidate parameters; Repeat the iteration until the preset maximum number of iterations is reached or the change in the objective function value is less than the preset threshold, and the optimized prediction model parameters are obtained.

2. The method for dynamic monitoring of soil trace elements for land consolidation according to claim 1, characterized in that: The performing distribution prediction based on the optimized prediction model parameters, the element spatial distribution diagram and the time variation characteristics to obtain the trace element prediction results and visualize them includes: activating the optimized prediction model according to the optimized prediction model parameters; Predicting the uncollected area within the target area according to the optimized prediction model to obtain an unknown prediction value; Performing Kriging correction based on the unknown predicted value in combination with the element spatial distribution map and the time variation characteristics to obtain a trace element prediction result; A dynamic evolution graph is generated according to the trace element prediction results according to the timestamps and visualized.

3. The method for dynamic monitoring of soil trace elements for land consolidation according to claim 1, characterized in that: Before obtaining the soil element dataset, the method further includes: Obtain soil texture distribution map; Dividing the regions according to the soil texture distribution map and calculating the coefficient of variation within different regions; When the coefficient of variation is greater than or equal to a preset coefficient of variation threshold, sampling is performed using a high-density sampling scheme; When the coefficient of variation is less than the coefficient of variation threshold, a low-density sampling scheme is adopted for sampling.

4. A soil trace element dynamic monitoring system for land consolidation, characterized in that: A method for dynamically monitoring soil trace elements for land remediation according to any one of claims 1 to 3, comprising: Data acquisition module, used to obtain soil element data set; A spatial distribution module is used to perform autocorrelation analysis on the soil element data set and perform spatial interpolation based on the obtained autocorrelation characteristics to obtain an element spatial distribution map; A time series analysis module, configured to perform time series analysis based on the soil element dataset to obtain time variation characteristics; A spatiotemporal coupling module, configured to perform spatiotemporal coupling according to the element spatial distribution diagram and the time variation characteristics to obtain spatiotemporal prediction model parameters; A model optimization module, configured to perform Bayesian optimization based on the spatiotemporal prediction model parameters to obtain optimized prediction model parameters; The result output module is used to perform distribution prediction based on the optimized prediction model parameters, the element spatial distribution map and the time variation characteristics, obtain the trace element prediction results, and visualize them.

5. An electronic device, characterized in that: The method comprises a processor, a memory, and a computer program stored in the memory and configured to be executed by the processor, wherein when the processor executes the computer program, the method for dynamic monitoring of soil trace elements for land remediation as described in any one of claims 1 to 3 is implemented.

6. A computer-readable storage medium, characterized in that The computer-readable storage medium includes a stored computer program, wherein when the computer program is running, the device where the computer-readable storage medium is located is controlled to execute the soil trace element dynamic monitoring method for land improvement according to any one of claims 1 to 3.

Images