A method for predicting soil moisture of multiple layers of farmland mulched soil

Abstract

The application discloses a mulched farmland multi-soil-layer soil moisture prediction method, and belongs to the technical field of soil moisture content prediction; the method comprises the following steps: acquiring multi-source data of a mulched farmland, wherein the multi-source data comprises soil moisture content data, meteorological data, soil property data, crop data, irrigation data and mulching management parameters; pre-processing the multi-source data, extracting input features and output target variables, and constructing a standardized modeling data set; taking the input features as input and the output target variables as output, training multiple machine learning prediction models, and screening a main prediction model and optimizing parameters thereof through cross-validation; the application realizes high-precision prediction of mulched farmland multi-soil-layer soil moisture and quantitative attribution of prediction sources by incorporating the mulching management parameters into the prediction model and combining soil layer explainability analysis, and provides technical support for mulched farmland water management and irrigation decision-making.

Description

Technical Field

[0001] This invention relates to the field of soil moisture prediction technology, specifically a method for predicting soil moisture in multi-soil-layered mulched farmland. Background Technology

[0002] Soil moisture is a key factor affecting crop growth, root water absorption, nutrient transport, and irrigation regime formulation. Accurately predicting the moisture dynamics of different soil layers is of great significance for precision irrigation and efficient water resource utilization.

[0003] Mulching affects the spatial and temporal distribution of soil moisture by inhibiting soil evaporation and regulating surface heat balance, essentially altering the water and heat exchange boundary conditions at the soil-atmosphere interface. However, the effect of mulching on soil moisture is not a single process. The mechanisms by which parameters such as mulch thickness, coverage, material, and color differentially influence the moisture content of shallow, middle, and deep soil layers remain unclear. Furthermore, the moisture effects of the same mulching measure vary significantly under different climatic, soil, and crop conditions. Therefore, accurate prediction of soil moisture across multiple soil layers under mulching conditions is a technical prerequisite for optimizing mulching schemes and irrigation regimes.

[0004] Existing soil moisture prediction technologies mainly include physical mechanism models, empirical statistical models, and machine learning models. Physical mechanism models (such as HYDRUS and SWAP) have a clear physical basis but their parameters are difficult to obtain; empirical statistical models are computationally simple but have limited generalization ability; machine learning models (such as support vector machines, random forests, and XGBoost) can handle multi-source data and nonlinear relationships, and some technologies further introduce interpretable analysis methods such as SHAP to interpret the model output.

[0005] However, the aforementioned existing technologies still have the following shortcomings:

[0006] First, there is a lack of adaptation to mulching scenarios and unified modeling for multiple soil layers. Existing methods do not treat mulching management parameters as independent input features, and cannot characterize the profile response characteristics of soil moisture under mulching conditions; at the same time, most of them are single-layer or surface layer predictions, lacking unified modeling for moisture in multiple soil layers, and lacking means to standardize soil layer data from different sources into a unified depth system.

[0007] Secondly, there is a lack of in-depth interpretability for mulching scenarios. Existing SHAP analysis provides a global interpretation and cannot reveal the differences in the contribution of the same variable in different soil layers and the law of its variation with depth; at the same time, there is a lack of quantitative analysis on the interaction effects between mulching management parameters, making it difficult to guide the refined design of mulching schemes.

[0008] Third, there is a lack of diagnostic analysis of the sources of prediction and validation of generalization ability. Existing methods lack a systematic diagnosis of the sources of model prediction ability and cannot quantitatively distinguish whether the prediction ability originates from temporal memory, water input, soil properties, or mulch regulation; moreover, they rely on cross-validation and lack external validation using independent field observation data and quantitative diagnosis of the reasons for the decline in generalization ability. Summary of the Invention

[0009] The purpose of this invention is to provide a method for predicting soil moisture in multi-soil-layered mulched farmland, in order to address the shortcomings mentioned in the background art.

[0010] To achieve the above objectives, the present invention provides the following technical solution: a method for predicting soil moisture in multi-layered mulched farmland, comprising the following steps:

[0011] Acquire multi-source data of mulched farmland, including soil moisture content data, meteorological data, soil property data, crop data, irrigation data, and mulching management parameters;

[0012] The multi-source data is preprocessed to extract input features and output target variables, and a standardized modeling dataset is constructed.

[0013] Using the input features as input and the output target variable as output, train multiple machine learning prediction models, select the master prediction model through cross-validation, and optimize its parameters;

[0014] The input features of the sample to be predicted are input into the optimized master prediction model to obtain the predicted soil moisture content values ​​for different soil layers.

[0015] Based on the optimized master prediction model, perform SHAP contribution analysis and / or variable group ablation diagnosis of soil layers, and output the variable contribution diagnosis results.

[0016] Preferably, the preprocessing includes:

[0017] When the soil layer division intervals from different data sources are inconsistent with the target standard depth interval, the original soil layer interval is defined by the upper and lower boundaries of the soil layer. The overlap thickness between the original soil layer interval and the target standard depth interval is calculated. The proportion of the overlap thickness to the thickness of the original soil layer interval or the thickness of the target standard depth interval is used as a weighting coefficient to perform a weighted calculation on the soil moisture content data within the original soil layer interval, thereby obtaining the soil moisture content data corresponding to the target standard depth interval.

[0018] Assign the soil moisture content data corresponding to the target standard depth range to the output target variable;

[0019] The different soil layers include shallow, intermediate, and deep layers; the shallow layer is a soil layer whose midpoint depth is less than or equal to a first preset threshold, the intermediate layer is a soil layer whose midpoint depth is greater than the first preset threshold and less than or equal to a second preset threshold, and the deep layer is a soil layer whose midpoint depth is greater than the second preset threshold; the midpoint depth of the soil layer is the arithmetic mean of the upper boundary and the lower boundary of the soil layer.

[0020] Preferably, the SHAP contribution analysis of the soil layers includes:

[0021] Based on the depth of the midpoint of the soil layer, the samples were divided into shallow sample subsets, middle sample subsets, and deep sample subsets.

[0022] The TreeSHAP algorithm is used to calculate the SHAP contribution value of each input feature in each sample subset to the prediction result.

[0023] The mean of the absolute values ​​of the SHAP contribution values ​​of all samples in each sample subset is calculated and used as the contribution intensity of the input feature in the soil layer; the mean of the SHAP contribution values ​​of all samples in each sample subset is directly calculated and used as the contribution direction of the input feature in the soil layer.

[0024] For each sample subset, sort the input features from largest to smallest contribution intensity to generate shallow contribution ranking, middle contribution ranking, and deep contribution ranking respectively; compare the changes in the position of the same input feature in different soil layers, and output the input features whose position rises above a preset position threshold or whose position falls below a preset position threshold as the key driving variable or non-sensitive variable of that soil layer.

[0025] Plot a line graph or scatter plot of the SHAP contribution value of each input feature as a function of depth, with the depth of the midpoint of the soil layer as the x-axis and the SHAP contribution value as the y-axis. Output the line graph or scatter plot as the variable contribution depth response curve.

[0026] Preferably, the variable group ablation diagnosis includes:

[0027] The input features are divided into time-series feature group, water input variable group, soil parameter variable group, mulch regulation variable group, vertical migration variable group, and climate / crop background variable group;

[0028] Remove all features corresponding to any variable group from the input features, keep the optimized hyperparameters of the main prediction model unchanged, retrain the model on the remaining features using the same hyperparameters of the main prediction model to obtain an updated model, and calculate the prediction accuracy of the updated model on the same test set.

[0029] The difference between the prediction accuracy of the main prediction model before removal and the prediction accuracy of the updated model after removal is used as the quantitative value of the contribution of the variable group to the model's prediction ability.

[0030] Preferably, the optimized master prediction model is further validated externally using independent field observation data.

[0031] Obtain independent field observation data, which comes from regions, years, or crop types that did not participate in the training of the master prediction model;

[0032] Extract the input features from the independent field observation data and input them into the optimized master prediction model to obtain the predicted soil moisture content at each observation time.

[0033] Calculate at least one of the following error indicators between the predicted and measured values ​​at each observation time: coefficient of determination R², root mean square error RMSE, and mean absolute error MAE.

[0034] Using the observation time as the horizontal axis and soil moisture content as the vertical axis, the rising and falling directions of the predicted and measured values ​​are compared moment by moment. The proportion of observation times with consistent rising and falling directions to the total number of observation times is counted as the trend direction consistency rate. The Pearson correlation coefficient and Nash efficiency coefficient (NSE) between the predicted and measured values ​​are calculated, and the trend direction consistency rate, the Pearson correlation coefficient, and the Nash efficiency coefficient are used together as the temporal consistency evaluation index. When the trend direction consistency rate is greater than a preset first threshold and the Nash efficiency coefficient is greater than a preset second threshold, the optimized master prediction model is determined to have passed the temporal consistency verification.

[0035] Preferably, the mulching management parameters include mulch film thickness, mulch film coverage, mulch film material, and mulch film color; the mulching management parameters are listed alongside the meteorological data, soil property data, crop data, and irrigation data, and are encoded together as components of the input features after preprocessing.

[0036] Preferably, the SHAP contribution analysis of soil layers further includes the SHAP contribution analysis of the interaction effects of mulch management parameters:

[0037] Select at least two film management parameters from the input features, normalize the numerical encoded values ​​of each film management parameter involved in the interaction, scale them to the same numerical range, and multiply the normalized value of the first film management parameter with the normalized value of the second film management parameter to construct the interaction features.

[0038] The interaction features are added to the input features, the model training is re-executed, and the TreeSHAP algorithm is used to calculate the SHAP contribution value of the interaction features to the prediction result.

[0039] The direction of the contribution of the interaction between the at least two mulch management parameters to the predicted soil moisture content is determined based on the sign of the SHAP contribution value of the interaction feature, and the contribution intensity is determined based on the absolute value of the SHAP contribution value of the interaction feature.

[0040] Preferably, scenario simulation and comparison of coating solutions are also included:

[0041] Under the condition of fixed meteorological data, soil property data, crop data and irrigation data, multiple mulching scheme combinations are constructed. Each mulching scheme combination includes different values ​​for mulch film thickness, mulch film coverage rate, mulch film material coding value and / or mulch film color coding value.

[0042] Each mulching scheme is combined with the fixed meteorological data, soil property data, crop data, and irrigation data to form the input features corresponding to each mulching scheme.

[0043] The input features corresponding to each mulching scheme are input into the optimized master prediction model to obtain the predicted soil moisture content of different soil layers under each mulching scheme.

[0044] For each mulching scheme, calculate the deviation score between the predicted soil moisture content of a specified soil layer and the preset target range, as well as the difference score between the shallow and deep soil moisture content. Calculate the weighted sum of the deviation score and the difference score to obtain the comprehensive score of the scheme. Sort all mulching schemes from high to low according to the comprehensive score and output the ranking results.

[0045] Preferably, it also includes the quantification of prediction uncertainty:

[0046] When performing the prediction of the optimized master prediction model, the quantile loss function of quantile regression forest or XGBoost is used, and at least two quantile levels are set to calculate the predicted soil moisture content of each soil layer at each quantile level; the at least two quantile levels include a first quantile and a second quantile, and the first quantile is less than the second quantile.

[0047] The interval between the predicted value corresponding to the first quantile and the predicted value corresponding to the second quantile is used as the confidence interval for the predicted value of soil moisture content in the soil layer.

[0048] Output the predicted soil moisture content for each soil layer and its corresponding confidence interval.

[0049] Preferably, the external validation further includes domain difference analysis between the training data and the external validation data:

[0050] Calculate the mean vector and covariance matrix of each input feature in the training data, and the mean vector and covariance matrix of the corresponding input features in the independent field observation data, respectively.

[0051] For each input feature, the difference between the mean of the independent field observation data and the mean of the training data is calculated as a first-order domain offset.

[0052] Calculate the difference matrix between the covariance matrix of the independent field observation data and the covariance matrix of the training data, and extract the diagonal elements corresponding to each feature in the difference matrix as second-order domain offsets;

[0053] The input features whose first-order domain offset is greater than a first preset threshold and the input features whose second-order domain offset is greater than a second preset threshold are jointly output as key domain offset variables that affect the generalization ability of the model.

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

[0055] 1. This invention encodes mulch film thickness, coverage rate, material, and color as independent mulch film management parameters alongside meteorological, soil, crop, and irrigation data as input features. Simultaneously, during the preprocessing stage, it performs overlapping thickness weighted calculations on soil layer intervals from different sources, uniformly mapping them to the target standard depth interval, thus constructing a three-layer prediction architecture of shallow, middle, and deep layers. The inclusion of mulch film management parameters enables the model to capture the changes in evaporation boundaries caused by mulch film, and the multi-layer soil architecture allows the model to distinguish the differences in moisture response at different depths.

[0056] 2. This invention uses the midpoint depth of the soil layer as the dividing criterion, dividing the sample into three subsets: shallow, middle, and deep. The SHAP contribution value of each input feature in each subset is calculated, generating a ranking of contributions for each soil layer, and a response curve showing the contribution value changing with depth is plotted. Simultaneously, the input features are divided into time-series feature groups, moisture input variable groups, soil parameter variable groups, mulch regulation variable groups, vertical migration variable groups, and climate / crop background variable groups. After removing each group, the model is retrained while keeping the hyperparameters unchanged. The difference in accuracy before and after removal is used as the quantified contribution of that group to the predictive ability. This allows for layer-specific SHAP to reveal the transfer patterns of dominant driving factors at different soil depths; variable group ablation quantifies the independent contribution of each variable group to the prediction accuracy.

[0057] 3. This invention selects at least two mulching management parameters from the input features, constructs interaction features after normalization, and calculates their SHAP contribution values ​​to quantitatively determine the synergistic or antagonistic effects between mulching parameters. Based on this, multiple mulching scheme combinations are constructed under the premise of fixing other conditions. The models are then input into each scheme to calculate the deviation score between the predicted value of the specified soil layer and the preset target interval, as well as the difference score between the shallow and deep layers. After weighted summation, a comprehensive score ranking is output. The synergistic mechanism between mulching parameters is revealed through interaction effect analysis, and this mechanism is used for scheme comparison in scenario simulation.

[0058] 4. This invention uses independent field observation data for external validation, and uses the trend direction consistency rate and Nash efficiency coefficient as time series consistency evaluation indicators, and sets a dual threshold joint judgment condition; further, it calculates the difference between the mean vector and the difference between the covariance matrix of the training data and the external validation data, extracts the first-order and second-order domain offsets of each feature, and identifies input features whose domain offsets exceed the preset threshold; external validation confirms the effectiveness of the model on independent data, and domain difference analysis points out the key domain offset variables that lead to the decay of generalization ability. Attached Figure Description

[0059] Figure 1 This is a flowchart of a method for predicting soil moisture in multi-soil-layered mulched farmland according to the present invention. Detailed Implementation

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

[0061] The embodiments of the present invention will be described below with reference to the accompanying drawings. It should be understood that the embodiments described are for illustrative purposes only and are not intended to limit the scope of protection of the present invention.

[0062] like Figure 1 As shown in the figure, the method for predicting soil moisture in multi-soil-layer farmland under plastic film mulching provided by this invention includes the following specific steps:

[0063] Step S1: Obtain multi-source data

[0064] Acquire multi-source data for mulched farmland. The multi-source data includes soil moisture content data, meteorological data, soil property data, crop data, irrigation data, and mulching management parameters.

[0065] Soil moisture content data includes soil moisture content at the target time, soil moisture content at the previous time, and the sampling interval. The previous time refers to the observation time determined according to the preset sampling interval before the target time, used to characterize the temporal continuity of soil moisture. The sampling interval is used to characterize the time difference between two adjacent observations. Depending on the data source, the preset sampling interval is one of 1 hour, 6 hours, 12 hours, 24 hours, or a longer interval.

[0066] Meteorological data includes at least one of precipitation, average temperature, and drought index. Precipitation is used to characterize external water supply, average temperature is used to characterize evaporation demand, and drought index is used to characterize the regional water supply and demand balance. The drought index is at least one of the following: Standardized Precipitation Evapotranspiration Index (SPEI), Palmer Drought Index (PDSI), or a drought index calculated from the ratio of precipitation to potential evapotranspiration.

[0067] Soil property data include soil organic matter content, soil bulk density, pH value, and soil texture. These parameters collectively characterize the soil's water retention capacity and infiltration performance.

[0068] Crop data includes at least one of crop type, crop coefficient, and days after sowing. The crop coefficient adopts the FAO56 crop coefficient to unify the water consumption characteristics of different crops and growth stages; the days after sowing is used to unify the time expression of data from different sources.

[0069] Irrigation data includes irrigation volume, which characterizes anthropogenic water input.

[0070] The mulching management parameters include mulch film thickness, mulch film coverage, mulch film material, and mulch film color. Mulch film thickness is measured in millimeters or micrometers, ranging from 0.005 mm to 0.030 mm; mulch film coverage is the proportion of mulch film-covered area to the total area of ​​the corresponding farmland observation unit, ranging from 0 to 1; mulch film material is selected from at least one of polyethylene film, polylactic acid biodegradable film, and starch-based biodegradable film; mulch film color is selected from at least one of transparent, black, white, and silver-gray.

[0071] The data used for model training comes from publicly published literature and field observation data. Each sample from different sources includes the study location, observation year, crop type, mulching treatment method, observation time, upper soil boundary, lower soil boundary, soil moisture at the target time, and meteorological, irrigation, and mulching management information corresponding to that observation time.

[0072] The aforementioned data come from different sources and have varying formats. Directly using them for model training can lead to problems such as conflicting variable names, inconsistent units, inability of categorical variables to participate in numerical calculations, and inconsistent soil layer classifications across different literature, resulting in a lack of unified modeling. Therefore, subsequent step S2, data preprocessing, is required to convert these multi-source data into a standardized dataset with a unified format and clear physical meaning, providing usable input-output pairs for model training.

[0073] Step S2: Data Preprocessing

[0074] The multi-source data obtained in step S1 is preprocessed to extract input features and output target variables, and a standardized modeling dataset is constructed.

[0075] Preprocessing includes the following operations: unifying fields with the same meaning from different sources into the same variable name and converting all variables to a unified unit; for continuous variables with a missing proportion below a preset threshold, imputation is performed using the mean, median, or linear interpolation, and samples with a missing proportion exceeding the preset threshold are removed; outliers are identified using the interquartile range (IQR) method, where IQR = Q3 - Q1, and observations below Q1 - 1.5 × IQR or above Q3 + 1.5 × IQR are marked as statistical outliers, and verified in conjunction with the physical range of the variables. Observations that are verified to be true extreme responses are retained or marked as outliers, while those that are recorded as errors are corrected or removed; one-hot encoding is performed on categorical variables, converting each category into a binary feature; for example, when the color of the mulch film includes transparent, black, and white, three binary features are generated respectively, and when the sample mulch film color is black, the above three features are set to 0, 1, and 0 respectively; mulch film material, soil texture, and crop type are all encoded in the same way.

[0076] Based on the above conventional pretreatment, this step also standardizes the soil depth to meet the needs of predicting multiple soil layers in mulched farmland.

[0077] Before standardization, soil moisture content data from different sources are uniformly converted into volumetric water content, with units of cm³ / cm³ or %. When the original data is mass water content, it is converted into volumetric water content based on soil bulk density and water density.

[0078] When the soil layer division intervals from different data sources are inconsistent with the target standard depth interval, the original soil layer interval is defined by the upper boundary UBD and the lower boundary LBD of the soil layer. The overlap thickness between the original soil layer interval and the target standard depth interval is calculated. The proportion of the overlap thickness to the thickness of the original soil layer interval or the thickness of the target standard depth interval is used as a weighting coefficient to perform weighted calculation on the soil moisture content data in the original soil layer interval, so as to obtain the soil moisture content data corresponding to the target standard depth interval.

[0079] In one embodiment, the depth of the midpoint of the soil layer is calculated using the following formula:

[0080]

[0081] in, The depth at the midpoint of the soil layer is represented by UBD, the upper boundary of the soil layer is represented by UBD, and the lower boundary of the soil layer is represented by LBD.

[0082] For example, if the original soil layer ranges from 0 to 30 cm and the target standard depth ranges from 0 to 20 cm (shallow layer), with a 20 cm overlap, the soil moisture content data from the 0 to 20 cm depth range is directly used as the target standard depth range data, with a weighting coefficient of 1.0 (the proportion of the overlap to the 20 cm target standard depth range). If the original soil layer ranges from 0 to 30 cm and the target standard depth ranges from 20 to 60 cm (middle layer), with a 10 cm overlap, the soil moisture content data from the 20 to 30 cm depth range is weighted and assigned to the target standard depth range data, with a weighting coefficient of 0.25 (the proportion of the overlap to the 40 cm target standard depth range).

[0083] After the soil depth standardization is completed, the soil moisture content data corresponding to the target standard depth range are assigned to the output target variable. Through the above standardization process, data from different literature, different experiments, or different monitoring depths can be unified under the same depth system, providing a data foundation for subsequent unified modeling of multiple soil layers.

[0084] Finally, the preprocessed data are uniformly encoded into input features X and output target variable Y. Input features X is a vector containing 19 features, including: soil moisture content at the previous time step, sampling interval, rainfall, average temperature, drought index, soil organic matter content, soil bulk density, pH value, soil texture code value, crop type code value, crop coefficient, number of days after sowing, irrigation amount, mulch film thickness, mulch film coverage rate, mulch film material code value, mulch film color code value, upper soil boundary value, and lower soil boundary value. Output target variable Y is the soil moisture content value of the corresponding soil layer at the target time.

[0085] For each continuous variable in the input feature X, Z-score standardization is used, and the calculation formula is as follows: , where x is the original value of the continuous variable, and μ and σ are the mean and standard deviation of the variable in the training set, respectively. The standardized parameters are calculated and stored only on the training set and applied to the validation set, test set, and subsequent samples to be predicted, thus avoiding data leakage.

[0086] At this point, the various raw data obtained in step S1, after the above preprocessing, have been unified into a standardized dataset with consistent format, clear physical meaning, and direct applicability for model training. This dataset consists of input features X and output target variable Y. All subsequent model training, comparison, selection, and optimization are conducted based on this dataset, ensuring data consistency and reproducibility in the modeling process.

[0087] Step S3: Train multiple machine learning prediction models

[0088] Based on the standardized modeling dataset constructed in step S2, multiple machine learning prediction models are established. Each model is trained using input feature X as input and output target variable Y as output. In one implementation, the input-output relationship of the above prediction models can be uniformly expressed as: ;

[0089] Where X represents the input feature vector, and Y represents the predicted soil moisture value at the target time. Indicates the main prediction model, Indicates model parameters.

[0090] The multiple machine learning prediction models include at least two of the following: multiple linear regression, support vector regression, random forest, and XGBoost. All models are trained and evaluated using the same data partitioning and the same K-fold cross-validation strategy to ensure fairness in model comparison.

[0091] Specifically, the standardized modeling dataset is divided into a training set and a test set using a fixed random seed of random_state=42. The training set is used for model fitting and hyperparameter optimization, while the test set is used to evaluate the generalization performance of each model. K-fold cross-validation is employed, dividing the training data into K subsets (K = 5). One subset is selected sequentially as the validation set, and the remaining K-1 subsets are used as the training set. This process is repeated K times to obtain the cross-validation prediction values ​​for each sample.

[0092] The initial parameter settings for each model are as follows: the multiple linear regression model uses least squares regression with an intercept term; the support vector regression model uses the radial basis kernel function; the random forest model sets the number of trees n_estimators=100 and the tree depth max_depth=None; the XGBoost model sets the number of trees n_estimators=100, the learning rate learning_rate=0.3, and the tree depth max_depth=6.

[0093] Using the coefficient of determination Root mean square error (RMSE) and mean absolute error (MAE) are used as evaluation indicators for prediction accuracy. In one implementation, these evaluation indicators are calculated using the following formula:

[0094]

[0095]

[0096]

[0097] in, This represents the measured soil moisture content of the i-th sample. This represents the predicted soil moisture content of the i-th sample. RMSE represents the average measured soil moisture content, and n represents the sample size; the dimensions of RMSE and MAE are consistent with those of soil moisture content.

[0098] Calculate the cross-validation results of each model. RMSE and MAE are used as the arithmetic mean of the cross-validation performance metrics for each model.

[0099] By following the steps above, under the same data partitioning and evaluation strategy, the prediction accuracy of different models on the same standardized dataset can be obtained, providing a comparable basis for the subsequent selection of the main prediction model. Compared to a single model, this step can objectively identify the model structure most suitable for predicting soil moisture in mulched farmland, avoiding the loss of prediction accuracy due to inappropriate model selection.

[0100] Step S4: Select the main prediction model and optimize the parameters

[0101] Based on the cross-validation results of step S3, a master predictor model is selected from multiple machine learning prediction models. In a preferred embodiment, the master predictor model is the XGBoost model, which performs best in cross-validation.

[0102] Hyperparameter optimization of the master prediction model is performed to further improve prediction accuracy. A Bayesian optimization method is used, with the goal of minimizing cross-validation error. The specific steps are as follows:

[0103] Select several initial parameter combinations, such as 10, within the preset hyperparameter space, and calculate the cross-validation error for each combination. Build a surrogate model based on the evaluated parameter combinations, using a Gaussian process or tree-structured Parzen estimation. Select the next set of parameters to be evaluated based on the acquisition function, such as the desired improvement in EI. Repeat the iteration until the preset number of iterations, such as 50, or the error convergence condition is reached, and output the optimized hyperparameter combination.

[0104] The hyperparameter search space of the XGBoost model is as follows: the number of trees n_estimators ranges from 50 to 500, the tree depth max_depth ranges from 3 to 15, the learning rate ranges from 0.01 to 0.3, the sample sampling ratio subsample ranges from 0.6 to 1.0, the feature sampling ratio colsample_bytree ranges from 0.6 to 1.0, the minimum child weight min_child_weight ranges from 1 to 10, the L1 regularization coefficient reg_alpha ranges from 0 to 10, the L2 regularization coefficient reg_lambda ranges from 0 to 10, and the split loss threshold gamma ranges from 0 to 5.

[0105] Through the above steps, the optimized master prediction model and its optimal hyperparameter combination are obtained, denoted as . This optimization process allows the model to maintain high accuracy while avoiding overfitting, providing a reliable model foundation for subsequent multi-soil layer prediction and interpretability analysis.

[0106] Step S5: Multi-layer prediction and uncertainty quantification

[0107] Input features of the sample to be predicted Substitute into the optimized master prediction model ;

[0108] in Based on the optimized hyperparameter combination from step S4, the predicted soil moisture content of different soil layers in the mulched farmland was calculated. .

[0109] The different soil layers include shallow, intermediate, and deep layers. The shallow layer is the soil layer with a midpoint depth less than or equal to a first preset threshold; the intermediate layer is the soil layer with a midpoint depth greater than the first preset threshold and less than or equal to a second preset threshold; and the deep layer is the soil layer with a midpoint depth greater than the second preset threshold.

[0110] In one implementation, the first preset threshold is 20cm, and the second preset threshold is 60cm; the depth of the midpoint of the soil layer is calculated using the formula in step S2. Obtained through calculation.

[0111] The input feature format of the sample to be predicted is consistent with that of the training phase. After prediction output, it is destandardized according to the mean and standard deviation saved during the training phase, and the output unit is volumetric water content in cm³ / cm³ or %.

[0112] Simultaneously with outputting the predicted values, the prediction uncertainty is quantified. The XGBoost quantile loss function is used, setting at least two quantile levels, for example, the first quantile is 0.1 and the second quantile is 0.9. The predicted soil moisture content for each soil layer is calculated at each quantile level. The mathematical form of the quantile loss function is:

[0113]

[0114] in, Quantile level, These are measured values. These are predicted values.

[0115] By minimizing the aforementioned loss function, the conditional quantile prediction value at this quantile level can be obtained. The interval between the predicted value corresponding to the first quantile and the predicted value corresponding to the second quantile is used as the confidence interval for the predicted soil moisture content of this soil layer.

[0116] Output the predicted soil moisture content for each soil layer and its corresponding confidence interval.

[0117] The predicted values ​​themselves provide the optimal estimate of the moisture content of each soil layer, while the confidence interval quantifies the uncertainty of this estimate—the width of the interval reflects the reliability of the prediction. For irrigation decisions, a wider confidence interval suggests that decision-makers need to incorporate more information or adopt a conservative strategy; a narrower confidence interval allows the prediction results to be directly used to formulate irrigation plans. This output format provides more complete input information for the scenario simulation and scheme comparison in subsequent step S8.

[0118] Step S6: SHAP Contribution Analysis of Soil Layers

[0119] Based on the optimized master prediction model in step S4, a soil layer SHAP contribution analysis was performed, specifically including the following steps:

[0120] Based on the depth of the midpoint of the soil layer, the samples were divided into shallow sample subsets, middle sample subsets, and deep sample subsets.

[0121] The TreeSHAP algorithm is used to calculate the SHAP contribution value of each input feature. This algorithm is specifically designed for tree models. By traversing each decision tree in the model, it assigns the prediction results to each input feature based on the Shapley value in game theory, and calculates the SHAP contribution value of each input feature in each sample subset to the prediction result.

[0122] For each sample subset, perform the following operations: Calculate the mean of the absolute values ​​of the SHAP contribution values ​​of all samples within the subset, and use this mean as the contribution intensity of the input feature in that soil layer; directly calculate the mean of the SHAP contribution values ​​of all samples within the subset, and use this mean as the contribution direction of the input feature in that soil layer. A positive contribution direction indicates that the input feature increases the predicted soil moisture content, while a negative contribution direction indicates that it decreases the predicted soil moisture content.

[0123] Furthermore, each sample subset is sorted from largest to smallest contribution intensity of each input feature, generating shallow contribution ranking, mid-layer contribution ranking, and deep contribution ranking respectively. The change in the ranking of the same input feature in different soil layers is compared—input features whose ranking rises above a preset ranking threshold are marked as key driving variables of that soil layer, and input features whose ranking falls below a preset ranking threshold are marked as non-sensitive variables of that soil layer.

[0124] The preset position threshold is determined based on the total number of input features. In one embodiment, when the total number of input features is 19, the preset position threshold is set to 3.

[0125] Plot a line graph or scatter plot of the SHAP contribution value of each input feature as a function of depth, with the depth of the midpoint of the soil layer as the x-axis and the SHAP contribution value as the y-axis. Output this as the variable contribution depth response curve.

[0126] The above steps provide a clear visual representation of the contribution of each input feature to soil moisture prediction as soil depth changes. A continuously decreasing curve for a particular feature indicates that its influence weakens with increasing depth; a peak in the curve within a certain depth range indicates that this range represents the specific soil layer where the feature plays a dominant role. This information can be directly used to guide the development of stratified management strategies.

[0127] Step S7: SHAP contribution analysis of the interaction effect of coating management parameters

[0128] Based on step S6, a SHAP contribution analysis of the interaction effect of the coating management parameters is further performed.

[0129] Select at least two film management parameters from the input features. Normalize the numerical encoded values ​​of each film management parameter involved in the interaction, scaling them to the same numerical range, such as [0,1], to eliminate scale differences between different parameters caused by different physical units, such as the film thickness in mm and the coverage rate as a dimensionless ratio of 0 to 1.

[0130] The interaction feature is constructed by multiplying the normalized values ​​of the first film management parameter and the second film management parameter.

[0131] Interactive features are added to the input features, the model is retrained, and the TreeSHAP algorithm is used to calculate the SHAP contribution value of the interactive features to the prediction results. The sign of the SHAP contribution value of the interactive features determines the direction of the interaction between at least two mulch management parameters on the predicted soil moisture content—positive values ​​indicate a synergistic effect, meaning the combination of the two parameters enhances each other's water retention effect; negative values ​​indicate an antagonistic effect, meaning the combination of the two parameters weakens each other's effect. The contribution strength is determined by the absolute value of the SHAP contribution value of the interactive features.

[0132] Through the above steps, we can quantitatively answer questions such as whether there is a synergistic effect between mulch film thickness and mulch film color. If the SHAP contribution value of the interaction feature between film thickness and color is positive and has a large absolute value, it indicates that the combined effect of "thick film + specific color" is better than the simple addition of their independent contributions; this analysis provides a quantitative basis for the construction and optimization of mulch film scheme combinations in the subsequent step S8.

[0133] Step S8: Simulation and Comparison of Lamination Solution Scenario

[0134] After completing model training and interpretability analysis, scenario simulation and comparison of coating schemes are conducted.

[0135] Under the condition of fixed meteorological data, soil property data, crop data and irrigation data, multiple mulching scheme combinations are constructed.

[0136] Each mulching scheme combination includes different values ​​for mulch film thickness, mulch film coverage, mulch film material coding, and / or mulch film color coding; for example, Scheme 1 is "black film + 0.010mm + 95% coverage", Scheme 2 is "transparent film + 0.010mm + 95% coverage", and Scheme 3 is "black film + 0.015mm + 95% coverage". Each mulching scheme combination is combined with fixed meteorological data, soil property data, crop data, and irrigation data to form the input features corresponding to each mulching scheme.

[0137] The input features corresponding to each coating scheme are respectively input into the optimized master prediction model. The predicted soil moisture content of each soil layer in the shallow, middle and deep layers under each mulching scheme was obtained.

[0138] For each mulching scheme, calculate the deviation score between the predicted soil moisture content of the specified soil layer and the preset target range, as well as the difference score between the shallow and deep soil moisture content.

[0139] The deviation score reflects whether the scheme can maintain the moisture of the target soil layer within the ideal range, while the difference score reflects whether the scheme will cause excessive differences in moisture between soil layers. The deviation score and the difference score are weighted and summed to obtain the comprehensive score of the scheme. All mulching schemes are ranked from high to low according to the comprehensive score, and the ranking results are output.

[0140] Through the steps described above, the synergistic or antagonistic mechanisms of the mulching parameters identified in step S7 are directly applied to scheme comparison—the scheme with the highest score is the mulching scheme with the best overall performance under the given soil, weather, and crop conditions. Users can quickly screen out the optimal mulching scheme without conducting field control trials, providing a direct quantitative decision-making basis for mulching scheme optimization.

[0141] Step S9: External Verification

[0142] The optimized master prediction model is externally validated using independent field observation data. This independent field observation data comes from regions, years, or crop types that were not involved in the training of the master prediction model; that is, the data source is completely independent of the training data in step S1 and is not involved in the model training and hyperparameter optimization in steps S3 and S4.

[0143] The specific verification process is as follows: Obtain independent field observation data, extract its input features and input them into the optimized master prediction model to obtain the predicted soil moisture content values ​​at each observation time; calculate the coefficient of determination between the predicted values ​​and the measured values ​​at each observation time. At least one of the following error indicators: root mean square error (RMSE) and mean absolute error (MAE), wherein the calculation formula for the above indicator is the same as that in step S3.

[0144] Based on this, a temporal consistency evaluation is conducted: with the observation time as the horizontal axis and soil moisture content as the vertical axis, the rising and falling directions of the predicted and measured values ​​are compared moment by moment; if the direction of change of the predicted value at the current moment relative to the predicted value at the previous moment is consistent with the direction of change of the measured value, it is counted as one instance of directional consistency; the proportion of observation times with consistent rising and falling directions to the total number of observation times is counted as the trend direction consistency rate; the Pearson correlation coefficient and Nash efficiency coefficient (NSE) between the predicted and measured values ​​are calculated.

[0145] The formula for calculating the Nash efficiency coefficient (NSE) is:

[0146]

[0147] in, Let be the measured value at the t-th observation time. Let be the predicted value at the t-th observation time. is the mean of the measured values, and T is the total number of observation times.

[0148] The trend direction consistency rate, Pearson correlation coefficient, and Nash efficiency coefficient are used as time series consistency evaluation indicators. When the trend direction consistency rate is greater than a preset first threshold and the Nash efficiency coefficient is greater than a preset second threshold, the optimized master prediction model is deemed to have passed the time series consistency verification. In one implementation, the preset first threshold is set to 70%, and the preset second threshold is set to 0.5.

[0149] The above verification steps objectively confirm the model's effectiveness and ability to capture time-series trends on independent data. This verification method differs from the cross-validation in step S3, which tests the model's stability within the training data; this step tests the model's generalization ability on entirely new data.

[0150] If the verification fails, further analysis is needed to determine the cause of the insufficient generalization ability, which is completed in step S10.

[0151] Step S10: Domain Difference Analysis

[0152] Based on the external validation in step S9, a domain difference analysis is performed between the training data and the external validation data.

[0153] Calculate the mean vector and covariance matrix of each input feature in the training data, and the mean vector and covariance matrix of the corresponding input features in the independent field observation data.

[0154] For each input feature, the difference between the mean of the independent field observation data and the mean of the training data is calculated as a first-order domain offset, which reflects the degree of offset of a single feature in the data center location.

[0155] Calculate the difference matrix between the covariance matrix of independent field observation data and the covariance matrix of training data. Extract the diagonal elements corresponding to each feature in the difference matrix as second-order domain offsets. The second-order domain offsets reflect the degree of shift of each feature in the data distribution structure, including the shrinkage or expansion of feature variance and the change in the correlation structure between features. This information cannot be captured by the first-order moments, i.e., the mean.

[0156] Input features with a first-order domain offset greater than a first preset threshold and input features with a second-order domain offset greater than a second preset threshold are used together as key domain offset variables that affect the generalization ability of the model.

[0157] Through the above steps, the external validation results of step S9 are expanded from a qualitative conclusion of whether the validation is passed to a quantitative diagnosis of which variable distribution differences lead to the decrease in accuracy. When the model is applied to a new region, the prediction accuracy decreases, and the output list directly points to the key variables that cause the generalization decay, such as "excessive differences in rainfall distribution" or "soil organic matter range exceeds the training set", providing a clear direction for local correction or resampling of the model.

[0158] Thus, a complete technical loop has been formed, from data acquisition, model training, multi-soil layer prediction, interpretable diagnosis to generalization verification and diagnosis.

[0159] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the invention can be implemented in other specific forms without departing from its spirit or essential characteristics. Therefore, the embodiments should be considered in all respects as exemplary and non-limiting, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of equivalents of the claims are intended to be included within the present invention. No reference numerals in the claims should be construed as limiting the scope of the claims.

Claims

1. A method for predicting soil moisture in multi-layered mulched farmland, characterized in that, Includes the following steps: Acquire multi-source data of mulched farmland, including soil moisture content data, meteorological data, soil property data, crop data, irrigation data, and mulching management parameters; The multi-source data is preprocessed to extract input features and output target variables, and a standardized modeling dataset is constructed. Using the input features as input and the output target variable as output, train multiple machine learning prediction models, select the master prediction model through cross-validation, and optimize its parameters. The input features of the sample to be predicted are input into the optimized master prediction model to obtain the predicted soil moisture content values ​​for different soil layers. Based on the optimized master prediction model, perform SHAP contribution analysis and / or variable group ablation diagnosis of soil layers, and output the variable contribution diagnosis results.

2. The method for predicting soil moisture in multi-soil-layered mulched farmland according to claim 1, characterized in that, The preprocessing includes: When the soil layer division intervals from different data sources are inconsistent with the target standard depth interval, the original soil layer interval is defined by the upper and lower boundaries of the soil layer. The overlap thickness between the original soil layer interval and the target standard depth interval is calculated. The proportion of the overlap thickness to the thickness of the original soil layer interval or the thickness of the target standard depth interval is used as a weighting coefficient to perform a weighted calculation on the soil moisture content data within the original soil layer interval, thereby obtaining the soil moisture content data corresponding to the target standard depth interval. Assign the soil moisture content data corresponding to the target standard depth range to the output target variable; The different soil layers include shallow, intermediate, and deep layers; the shallow layer is a soil layer whose midpoint depth is less than or equal to a first preset threshold, the intermediate layer is a soil layer whose midpoint depth is greater than the first preset threshold and less than or equal to a second preset threshold, and the deep layer is a soil layer whose midpoint depth is greater than the second preset threshold; the midpoint depth of the soil layer is the arithmetic mean of the upper boundary and the lower boundary of the soil layer.

3. The method for predicting soil moisture in multi-layered mulched farmland according to claim 1, characterized in that, The SHAP contribution analysis of the soil layers includes: Based on the depth of the midpoint of the soil layer, the samples were divided into shallow sample subsets, middle sample subsets, and deep sample subsets. The TreeSHAP algorithm is used to calculate the SHAP contribution value of each input feature in each sample subset to the prediction result. The mean of the absolute values ​​of the SHAP contribution values ​​of all samples in each sample subset is calculated and used as the contribution intensity of the input feature in the soil layer; the mean of the SHAP contribution values ​​of all samples in each sample subset is directly calculated and used as the contribution direction of the input feature in the soil layer. For each sample subset, sort the input features from largest to smallest contribution intensity to generate shallow contribution ranking, middle contribution ranking, and deep contribution ranking respectively; compare the changes in the position of the same input feature in different soil layers, and output the input features whose position rises above a preset position threshold or whose position falls below a preset position threshold as the key driving variable or non-sensitive variable of that soil layer. Plot a line graph or scatter plot of the SHAP contribution value of each input feature as a function of depth, with the depth of the midpoint of the soil layer as the x-axis and the SHAP contribution value as the y-axis. Output the line graph or scatter plot as the variable contribution depth response curve.

4. The method for predicting soil moisture in multi-layered mulched farmland according to claim 1, characterized in that, The variable group ablation diagnosis includes: The input features are divided into time-series feature group, water input variable group, soil parameter variable group, mulch regulation variable group, vertical migration variable group, and climate / crop background variable group; Remove all features corresponding to any variable group from the input features, keep the optimized hyperparameters of the main prediction model unchanged, retrain the model on the remaining features using the same hyperparameters of the main prediction model to obtain an updated model, and calculate the prediction accuracy of the updated model on the same test set. The difference between the prediction accuracy of the main prediction model before removal and the prediction accuracy of the updated model after removal is used as the quantitative value of the contribution of the variable group to the model's prediction ability.

5. The method for predicting soil moisture in multi-layered mulched farmland according to claim 1, characterized in that, This also includes external validation of the optimized master prediction model using independent field observation data: Obtain independent field observation data, which comes from regions, years, or crop types that did not participate in the training of the master prediction model; Extract the input features from the independent field observation data and input them into the optimized master prediction model to obtain the predicted soil moisture content at each observation time. Calculate at least one of the following error indicators between the predicted and measured values ​​at each observation time: coefficient of determination R², root mean square error RMSE, and mean absolute error MAE. Using the observation time as the horizontal axis and soil moisture content as the vertical axis, the rising and falling directions of the predicted value and the measured value are compared moment by moment. The proportion of observation times with consistent upward and downward directions to the total number of observation times is used as the trend direction consistency rate; the Pearson correlation coefficient and Nash efficiency coefficient (NSE) between the predicted value and the measured value are calculated, and the trend direction consistency rate, the Pearson correlation coefficient, and the Nash efficiency coefficient are used together as the time series consistency evaluation index. When the trend direction consistency rate is greater than a preset first threshold and the Nash efficiency coefficient is greater than a preset second threshold, the optimized master prediction model is determined to have passed the time series consistency verification.

6. A method for predicting soil moisture in multi-layered mulched farmland according to any one of claims 1 to 3, characterized in that, The mulching management parameters include mulch film thickness, mulch film coverage, mulch film material, and mulch film color; these mulching management parameters, along with the meteorological data, soil property data, crop data, and irrigation data, are preprocessed and encoded together as components of the input features.

7. The method for predicting soil moisture in multi-soil-layered mulched farmland according to claim 3, characterized in that, The SHAP contribution analysis for soil layers also includes the SHAP contribution analysis of the interaction effects of mulch management parameters: Select at least two film management parameters from the input features, normalize the numerical encoded values ​​of each film management parameter involved in the interaction, scale them to the same numerical range, and multiply the normalized value of the first film management parameter with the normalized value of the second film management parameter to construct the interaction features. The interaction features are added to the input features, the model training is re-executed, and the TreeSHAP algorithm is used to calculate the SHAP contribution value of the interaction features to the prediction result. The direction of the contribution of the interaction between the at least two mulch management parameters to the predicted soil moisture content is determined based on the sign of the SHAP contribution value of the interaction feature, and the contribution intensity is determined based on the absolute value of the SHAP contribution value of the interaction feature.

8. The method for predicting soil moisture in multi-soil-layered mulched farmland according to claim 6, characterized in that, It also includes scenario simulation and comparison of coating solutions: Under the condition of fixed meteorological data, soil property data, crop data and irrigation data, multiple mulching scheme combinations are constructed. Each mulching scheme combination includes different values ​​for mulch film thickness, mulch film coverage rate, mulch film material coding value and / or mulch film color coding value. Each mulching scheme is combined with the fixed meteorological data, soil property data, crop data, and irrigation data to form the input features corresponding to each mulching scheme. The input features corresponding to each mulching scheme are input into the optimized master prediction model to obtain the predicted soil moisture content of different soil layers under each mulching scheme. For each mulching scheme, calculate the deviation score between the predicted soil moisture content of a specified soil layer and the preset target range, as well as the difference score between the shallow and deep soil moisture content; then, sum the deviation score and the difference score by weight to obtain the comprehensive score of the scheme. All coating schemes are sorted from highest to lowest according to the comprehensive score, and the sorting results are output.

9. The method for predicting soil moisture in multi-soil-layered mulched farmland according to claim 1, characterized in that, It also includes the quantification of prediction uncertainty: When performing the prediction of the optimized master prediction model, the quantile loss function of quantile regression forest or XGBoost is used, and at least two quantile levels are set to calculate the predicted soil moisture content of each soil layer at each quantile level; the at least two quantile levels include a first quantile and a second quantile, and the first quantile is less than the second quantile. The interval between the predicted value corresponding to the first quantile and the predicted value corresponding to the second quantile is used as the confidence interval for the predicted value of soil moisture content in the soil layer. Output the predicted soil moisture content for each soil layer and its corresponding confidence interval.

10. A method for predicting soil moisture in multi-layered mulched farmland according to claim 5, characterized in that, The external validation also includes domain difference analysis between the training data and the external validation data: Calculate the mean vector and covariance matrix of each input feature in the training data, and the mean vector and covariance matrix of the corresponding input features in the independent field observation data, respectively. For each input feature, the difference between the mean of the independent field observation data and the mean of the training data is calculated as a first-order domain offset. Calculate the difference matrix between the covariance matrix of the independent field observation data and the covariance matrix of the training data, and extract the diagonal elements corresponding to each feature in the difference matrix as second-order domain offsets; The input features whose first-order domain offset is greater than a first preset threshold and the input features whose second-order domain offset is greater than a second preset threshold are jointly output as key domain offset variables that affect the generalization ability of the model.

Images